.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/quadtree_gridded_forecast_evaluation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_quadtree_gridded_forecast_evaluation.py: .. _quadtree_gridded-forecast-evaluation: Quadtree Grid-based Forecast Evaluation ======================================= This example demonstrates how to create a quadtree based single resolution-grid and multi-resolution grid. Multi-resolution grid is created using earthquake catalog, in which seismic density determines the size of a grid cell. In creating a multi-resolution grid we select a threshold (:math:`N_{max}`) as a maximum number of earthquake in each cell. In single-resolution grid, we simply select a zoom-level (L) that determines a single resolution grid. The number of cells in single-resolution grid are equal to :math:`4^L`. The zoom-level L=11 leads to 4.2 million cells, nearest to 0.1x0.1 grid. We use these grids to create and evaluate a time-independent forecast. Grid-based forecasts assume the variability of the forecasts is Poissonian. Therefore, poisson-based evaluations should be used to evaluate grid-based forecasts defined using quadtree regions. Overview: 1. Define spatial grids - Multi-resolution grid - Single-resolution grid 2. Load forecasts - Multi-resolution forecast - Single-resolution forecast 3. Load evaluation catalog 4. Apply Poissonian evaluations for both grid-based forecasts 5. Visualize evaluation results .. GENERATED FROM PYTHON SOURCE LINES 31-36 Load required libraries ----------------------- Most of the core functionality can be imported from the top-level :mod:`csep` package. Utilities are available from the :mod:`csep.utils` subpackage. .. GENERATED FROM PYTHON SOURCE LINES 36-45 .. code-block:: Python import numpy import pandas from csep.core import poisson_evaluations as poisson from csep.utils import time_utils, plots from csep.core.regions import QuadtreeGrid2D from csep.core.forecasts import GriddedForecast from csep.utils.time_utils import decimal_year_to_utc_epoch from csep.core.catalogs import CSEPCatalog .. GENERATED FROM PYTHON SOURCE LINES 46-52 Load Training Catalog for Multi-resolution grid ---------------------------------------------- We define a multi-resolution quadtree using earthquake catalog. We load a training catalog in CSEP and use that catalog to create a multi-resolution grid. Sometimes, we do not the catalog in exact format as requried by pyCSEP. So we can read a catalog using Pandas and convert it into the format accepable by PyCSEP. Then we instantiate an object of class CSEPCatalog by calling function :func:`csep.core.regions.CSEPCatalog.from_dataframe` .. GENERATED FROM PYTHON SOURCE LINES 52-71 .. code-block:: Python dfcat = pandas.read_csv('cat_train_2013.csv') column_name_mapper = { 'lon': 'longitude', 'lat': 'latitude', 'mag': 'magnitude', 'index': 'id' } # maps the column names to the dtype expected by the catalog class dfcat = dfcat.reset_index().rename(columns=column_name_mapper) # create the origin_times from decimal years dfcat['origin_time'] = dfcat.apply(lambda row: decimal_year_to_utc_epoch(row.year), axis=1) # create catalog from dataframe catalog_train = CSEPCatalog.from_dataframe(dfcat) print(catalog_train) .. rst-class:: sphx-glr-script-out .. code-block:: none Name: None Start Date: 1976-01-01 00:00:00+00:00 End Date: 2013-01-01 00:00:00+00:00 Latitude: (-77.16000366, 87.01999664) Longitude: (-180.0, 180.0) Min Mw: 5.150024414 Max Mw: 9.08350563 Event Count: 28465 .. GENERATED FROM PYTHON SOURCE LINES 72-77 Define Multi-resolution Gridded Region ------------------------------------------------ Now use define a threshold for maximum number of earthquake allowed per cell, i.e. Nmax and call :func:`csep.core.regions.QuadtreeGrid_from_catalog` to create a multi-resolution grid. For simplicity we assume only single magnitude bin, i.e. all the earthquakes greater than and equal to 5.95 .. GENERATED FROM PYTHON SOURCE LINES 77-84 .. code-block:: Python mbins = numpy.array([5.95]) Nmax = 25 r_multi = QuadtreeGrid2D.from_catalog(catalog_train, Nmax, magnitudes=mbins) print('Number of cells in Multi-resolution grid :', r_multi.num_nodes) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of cells in Multi-resolution grid : 3502 .. GENERATED FROM PYTHON SOURCE LINES 85-90 Define Single-resolution Gridded Region ---------------------------------------- Here as an example we define a single resolution grid at zoom-level L=6. For this purpose we call :func:`csep.core.regions.QuadtreeGrid2D_from_single_resolution` to create a single resolution grid. .. GENERATED FROM PYTHON SOURCE LINES 90-99 .. code-block:: Python # For simplicity of example, we assume only single magnitude bin, # i.e. all the earthquakes greater than and equal to 5.95 mbins = numpy.array([5.95]) r_single = QuadtreeGrid2D.from_single_resolution(6, magnitudes=mbins) print('Number of cells in Single-Resolution grid :', r_single.num_nodes) .. rst-class:: sphx-glr-script-out .. code-block:: none Number of cells in Single-Resolution grid : 4096 .. GENERATED FROM PYTHON SOURCE LINES 100-105 Load forecast of multi-resolution grid --------------------------------- An example time-independent forecast had been created for this grid and provided the example forecast data set along with the main repository. We load the time-independent global forecast which has time horizon of 1 year. The filepath is relative to the root directory of the package. You can specify any file location for your forecasts. .. GENERATED FROM PYTHON SOURCE LINES 105-119 .. code-block:: Python forecast_data = numpy.loadtxt('example_rate_zoom=EQ10L11.csv') #Reshape forecast as Nx1 array forecast_data = forecast_data.reshape(-1,1) forecast_multi_grid = GriddedForecast(data = forecast_data, region = r_multi, magnitudes = mbins, name = 'Example Multi-res Forecast') #The loaded forecast is for 1 year. The test catalog we will use to evaluate is for 6 years. So we can rescale the forecast. print(f"expected event count before scaling: {forecast_multi_grid.event_count}") forecast_multi_grid.scale(6) print(f"expected event count after scaling: {forecast_multi_grid.event_count}") .. rst-class:: sphx-glr-script-out .. code-block:: none expected event count before scaling: 116.18568954606255 expected event count after scaling: 697.1141372763753 .. GENERATED FROM PYTHON SOURCE LINES 120-124 Load forecast of single-resolution grid ------------------------------------- We have already created a time-independent global forecast with time horizon of 1 year and provided with the reporsitory. The filepath is relative to the root directory of the package. You can specify any file location for your forecasts. .. GENERATED FROM PYTHON SOURCE LINES 124-139 .. code-block:: Python forecast_data = numpy.loadtxt('example_rate_zoom=6.csv') #Reshape forecast as Nx1 array forecast_data = forecast_data.reshape(-1,1) forecast_single_grid = GriddedForecast(data = forecast_data, region = r_single, magnitudes = mbins, name = 'Example Single-res Forecast') # The loaded forecast is for 1 year. The test catalog we will use is for 6 years. So we can rescale the forecast. print(f"expected event count before scaling: {forecast_single_grid.event_count}") forecast_single_grid.scale(6) print(f"expected event count after scaling: {forecast_single_grid.event_count}") .. rst-class:: sphx-glr-script-out .. code-block:: none expected event count before scaling: 116.18568954606256 expected event count after scaling: 697.1141372763753 .. GENERATED FROM PYTHON SOURCE LINES 140-145 Load evaluation catalog ----------------------- We have a test catalog stored here. We can read the test catalog as a pandas frame and convert it into a format that is acceptable to PyCSEP Then we instantiate an object of catalog .. GENERATED FROM PYTHON SOURCE LINES 145-163 .. code-block:: Python dfcat = pandas.read_csv('cat_test.csv') column_name_mapper = { 'lon': 'longitude', 'lat': 'latitude', 'mag': 'magnitude' } # maps the column names to the dtype expected by the catalog class dfcat = dfcat.reset_index().rename(columns=column_name_mapper) # create the origin_times from decimal years dfcat['origin_time'] = dfcat.apply(lambda row: decimal_year_to_utc_epoch(row.year), axis=1) # create catalog from dataframe catalog = CSEPCatalog.from_dataframe(dfcat) print(catalog) .. rst-class:: sphx-glr-script-out .. code-block:: none Name: None Start Date: 2014-01-01 00:00:00+00:00 End Date: 2019-01-01 00:00:00+00:00 Latitude: (-63.26, 74.39) Longitude: (-179.23, 179.66) Min Mw: 5.95047692260089 Max Mw: 8.27271203001144 Event Count: 651 .. GENERATED FROM PYTHON SOURCE LINES 164-173 Compute Poisson spatial test and Number test ------------------------------------------------------ Simply call the :func:`csep.core.poisson_evaluations.spatial_test` and :func:`csep.core.poisson_evaluations.number_test` functions to evaluate the forecast using the specified evaluation catalog. The spatial test requires simulating from the Poisson forecast to provide uncertainty. The verbose option prints the status of the simulations to the standard output. Note: But before we use evaluation catalog, we need to link gridded region with observed catalog. Since we have two different grids here, so we do it separately for both grids. .. GENERATED FROM PYTHON SOURCE LINES 173-186 .. code-block:: Python #For Multi-resolution grid, linking region to catalog. catalog.region = forecast_multi_grid.region spatial_test_multi_res_result = poisson.spatial_test(forecast_multi_grid, catalog) number_test_multi_res_result = poisson.number_test(forecast_multi_grid, catalog) #For Single-resolution grid, linking region to catalog. catalog.region = forecast_single_grid.region spatial_test_single_res_result = poisson.spatial_test(forecast_single_grid, catalog) number_test_single_res_result = poisson.number_test(forecast_single_grid, catalog) .. GENERATED FROM PYTHON SOURCE LINES 187-192 Plot spatial test results ------------------------- We provide the function :func:`csep.utils.plotting.plot_poisson_consistency_test` to visualize the evaluation results from consistency tests. .. GENERATED FROM PYTHON SOURCE LINES 192-201 .. code-block:: Python stest_result = [spatial_test_single_res_result, spatial_test_multi_res_result] ax_spatial = plots.plot_poisson_consistency_test(stest_result, plot_args={'xlabel': 'Spatial likelihood'}) ntest_result = [number_test_single_res_result, number_test_multi_res_result] ax_number = plots.plot_poisson_consistency_test(ntest_result, plot_args={'xlabel': 'Number of Earthquakes'}) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /tutorials/images/sphx_glr_quadtree_gridded_forecast_evaluation_001.png :alt: Poisson S-Test :srcset: /tutorials/images/sphx_glr_quadtree_gridded_forecast_evaluation_001.png :class: sphx-glr-multi-img * .. image-sg:: /tutorials/images/sphx_glr_quadtree_gridded_forecast_evaluation_002.png :alt: Poisson N-Test :srcset: /tutorials/images/sphx_glr_quadtree_gridded_forecast_evaluation_002.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.797 seconds) .. _sphx_glr_download_tutorials_quadtree_gridded_forecast_evaluation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: quadtree_gridded_forecast_evaluation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: quadtree_gridded_forecast_evaluation.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_