.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/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_gridded_forecast_evaluation.py: .. _grid-forecast-evaluation: Grid-based Forecast Evaluation ============================== This example demonstrates how to evaluate a grid-based and 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. Overview: 1. Define forecast properties (time horizon, spatial region, etc). 2. Obtain evaluation catalog 3. Apply Poissonian evaluations for grid-based forecasts 4. Store evaluation results using JSON format 5. Visualize evaluation results .. GENERATED FROM PYTHON SOURCE LINES 21-26 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 26-34 .. code-block:: Python import csep from csep.core import poisson_evaluations as poisson from csep.utils import datasets, time_utils, plots # Needed to show plots from the terminal import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 35-41 Define forecast properties -------------------------- We choose a :ref:`time-independent-forecast` to show how to evaluate a grid-based earthquake forecast using PyCSEP. Note, the start and end date should be chosen based on the creation of the forecast. This is important for time-independent forecasts because they can be rescale to any arbitrary time period. .. GENERATED FROM PYTHON SOURCE LINES 41-46 .. code-block:: Python from csep.utils.stats import get_Kagan_I1_score start_date = time_utils.strptime_to_utc_datetime('2006-11-12 00:00:00.0') end_date = time_utils.strptime_to_utc_datetime('2011-11-12 00:00:00.0') .. GENERATED FROM PYTHON SOURCE LINES 47-52 Load forecast ------------- For this example, we provide the example forecast data set along with the main repository. 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 52-58 .. code-block:: Python forecast = csep.load_gridded_forecast(datasets.helmstetter_aftershock_fname, start_date=start_date, end_date=end_date, name='helmstetter_aftershock') .. GENERATED FROM PYTHON SOURCE LINES 59-65 Load evaluation catalog ----------------------- We will download the evaluation catalog from ComCat (this step requires an internet connection). We can use the ComCat API to filter the catalog in both time and magnitude. See the catalog filtering example, for more information on how to filter the catalog in space and time manually. .. GENERATED FROM PYTHON SOURCE LINES 65-70 .. code-block:: Python print("Querying comcat catalog") catalog = csep.query_comcat(forecast.start_time, forecast.end_time, min_magnitude=forecast.min_magnitude) print(catalog) .. rst-class:: sphx-glr-script-out .. code-block:: none Querying comcat catalog Fetched ComCat catalog in 5.163585662841797 seconds. Downloaded catalog from ComCat with following parameters Start Date: 2007-02-26 12:19:54.530000+00:00 End Date: 2011-02-18 17:47:35.770000+00:00 Min Latitude: 31.9788333 and Max Latitude: 41.1444 Min Longitude: -125.0161667 and Max Longitude: -114.8398 Min Magnitude: 4.96 Found 34 events in the ComCat catalog. Name: None Start Date: 2007-02-26 12:19:54.530000+00:00 End Date: 2011-02-18 17:47:35.770000+00:00 Latitude: (31.9788333, 41.1444) Longitude: (-125.0161667, -114.8398) Min Mw: 4.96 Max Mw: 7.2 Event Count: 34 .. GENERATED FROM PYTHON SOURCE LINES 71-75 Filter evaluation catalog in space ---------------------------------- We need to remove events in the evaluation catalog outside the valid region specified by the forecast. .. GENERATED FROM PYTHON SOURCE LINES 75-79 .. code-block:: Python catalog = catalog.filter_spatial(forecast.region) print(catalog) .. rst-class:: sphx-glr-script-out .. code-block:: none Name: None Start Date: 2007-02-26 12:19:54.530000+00:00 End Date: 2011-02-18 17:47:35.770000+00:00 Latitude: (31.9788333, 41.1155) Longitude: (-125.0161667, -115.0481667) Min Mw: 4.96 Max Mw: 7.2 Event Count: 32 .. GENERATED FROM PYTHON SOURCE LINES 80-86 Compute Poisson spatial test ---------------------------- Simply call the :func:`csep.core.poisson_evaluations.spatial_test` function 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. .. GENERATED FROM PYTHON SOURCE LINES 86-89 .. code-block:: Python spatial_test_result = poisson.spatial_test(forecast, catalog) .. GENERATED FROM PYTHON SOURCE LINES 90-95 Store evaluation results ------------------------ PyCSEP provides easy ways of storing objects to a JSON format using :func:`csep.write_json`. The evaluations can be read back into the program for plotting using :func:`csep.load_evaluation_result`. .. GENERATED FROM PYTHON SOURCE LINES 95-98 .. code-block:: Python csep.write_json(spatial_test_result, 'example_spatial_test.json') .. GENERATED FROM PYTHON SOURCE LINES 99-104 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 104-109 .. code-block:: Python ax = plots.plot_poisson_consistency_test(spatial_test_result, plot_args={'xlabel': 'Spatial likelihood'}) plt.show() .. image-sg:: /tutorials/images/sphx_glr_gridded_forecast_evaluation_001.png :alt: Poisson S-Test :srcset: /tutorials/images/sphx_glr_gridded_forecast_evaluation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 110-123 Plot ROC Curves --------------- We can also plot the Receiver operating characteristic (ROC) Curves based on forecast and testing-catalog. In the figure below, True Positive Rate is the normalized cumulative forecast rate, after sorting cells in decreasing order of rate. The “False Positive Rate” is the normalized cumulative area. The dashed line is the ROC curve for a uniform forecast, meaning the likelihood for an earthquake to occur at any position is the same. The further the ROC curve of a forecast is to the uniform forecast, the specific the forecast is. When comparing the forecast ROC curve against a catalog, one can evaluate if the forecast is more or less specific (or smooth) at different level or seismic rate. Note: This figure just shows an example of plotting an ROC curve with a catalog forecast. If "linear=True" the diagram is represented using a linear x-axis. If "linear=False" the diagram is represented using a logarithmic x-axis. .. GENERATED FROM PYTHON SOURCE LINES 123-131 .. code-block:: Python print("Plotting concentration ROC curve") _= plots.plot_concentration_ROC_diagram(forecast, catalog, linear=True) .. image-sg:: /tutorials/images/sphx_glr_gridded_forecast_evaluation_002.png :alt: Concentration ROC Curve :srcset: /tutorials/images/sphx_glr_gridded_forecast_evaluation_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Plotting concentration ROC curve .. GENERATED FROM PYTHON SOURCE LINES 132-134 Plot ROC and Molchan curves using the alarm-based approach ---------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 134-155 .. code-block:: Python #In this script, we generate ROC diagrams and Molchan diagrams using the alarm-based approach to evaluate the predictive #performance of models. This method exploits contingency table analysis to evaluate the predictive capabilities of #forecasting models. By analysing the contingency table data, we determine the ROC curve and Molchan trajectory and #estimate the Area Skill Score to assess the accuracy and reliability of the prediction models. The generated graphs #visually represent the prediction performance. # Note: If "linear=True" the diagram is represented using a linear x-axis. # If "linear=False" the diagram is represented using a logarithmic x-axis. print("Plotting ROC curve from the contingency table") # Set linear True to obtain a linear x-axis, False to obtain a logical x-axis. _ = plots.plot_ROC_diagram(forecast, catalog, linear=True) print("Plotting Molchan curve from the contingency table and the Area Skill Score") # Set linear True to obtain a linear x-axis, False to obtain a logical x-axis. _ = plots.plot_Molchan_diagram(forecast, catalog, linear=True) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /tutorials/images/sphx_glr_gridded_forecast_evaluation_003.png :alt: ROC Curve from contingency table :srcset: /tutorials/images/sphx_glr_gridded_forecast_evaluation_003.png :class: sphx-glr-multi-img * .. image-sg:: /tutorials/images/sphx_glr_gridded_forecast_evaluation_004.png :alt: gridded forecast evaluation :srcset: /tutorials/images/sphx_glr_gridded_forecast_evaluation_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Plotting ROC curve from the contingency table Plotting Molchan curve from the contingency table and the Area Skill Score .. GENERATED FROM PYTHON SOURCE LINES 156-161 Calculate Kagan's I_1 score --------------------------- We can also get the Kagan's I_1 score for a gridded forecast (see Kagan, YanY. [2009] Testing long-term earthquake forecasts: likelihood methods and error diagrams, Geophys. J. Int., v.177, pages 532-542). .. GENERATED FROM PYTHON SOURCE LINES 161-164 .. code-block:: Python I_1 = get_Kagan_I1_score(forecast, catalog) print("I_1score is: ", I_1) .. rst-class:: sphx-glr-script-out .. code-block:: none I_1score is: [2.31435371] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 19.388 seconds) .. _sphx_glr_download_tutorials_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: gridded_forecast_evaluation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: gridded_forecast_evaluation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: gridded_forecast_evaluation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_