Working with catalog-based forecasts

This example shows some basic interactions with data-based forecasts. We will load in a forecast stored in the CSEP data format, and compute the expected rates on a 0.1° x 0.1° grid covering the state of California. We will plot the expected rates in the spatial cells.

  1. Define forecast properties (time horizon, spatial region, etc).

  2. Compute the expected rates in space and magnitude bins

  3. Plot expected rates in the spatial cells

Load required libraries

Most of the core functionality can be imported from the top-level csep package. Utilities are available from the csep.utils subpackage.

import numpy

import csep
from csep.core import regions
from csep.utils import datasets

Load data forecast

PyCSEP contains some basic forecasts that can be used to test of the functionality of the package. This forecast has already been filtered to the California RELM region.

Define spatial and magnitude regions

Before we can compute the bin-wise rates we need to define a spatial region and a set of magnitude bin edges. The magnitude bin edges # are the lower bound (inclusive) except for the last bin, which is treated as extending to infinity. We can bind these # to the forecast object. This can also be done by passing them as keyword arguments into csep.load_catalog_forecast().

# Magnitude bins properties
min_mw = 4.95
max_mw = 8.95
dmw = 0.1

# Create space and magnitude regions
magnitudes = regions.magnitude_bins(min_mw, max_mw, dmw)
region = regions.california_relm_region()

# Bind region information to the forecast (this will be used for binning of the catalogs)
forecast.region = regions.create_space_magnitude_region(region, magnitudes)

Compute spatial event counts

The csep.core.forecasts.CatalogForecast provides a method to compute the expected number of events in spatial cells. This requires a region with magnitude information.


Processed 1 catalogs in 0.0020678043365478516 seconds
Processed 2 catalogs in 0.003952503204345703 seconds
Processed 3 catalogs in 0.0055806636810302734 seconds
Processed 4 catalogs in 0.007014274597167969 seconds
Processed 5 catalogs in 0.00838613510131836 seconds
Processed 6 catalogs in 0.00983428955078125 seconds
Processed 7 catalogs in 0.011170148849487305 seconds
Processed 8 catalogs in 0.012708902359008789 seconds
Processed 9 catalogs in 0.015000104904174805 seconds
Processed 10 catalogs in 0.016472816467285156 seconds
Processed 20 catalogs in 0.02996969223022461 seconds
Processed 30 catalogs in 0.04437565803527832 seconds
Processed 40 catalogs in 0.05838727951049805 seconds
Processed 50 catalogs in 0.07305192947387695 seconds
Processed 60 catalogs in 0.08679389953613281 seconds
Processed 70 catalogs in 0.10055685043334961 seconds
Processed 80 catalogs in 0.11487603187561035 seconds
Processed 90 catalogs in 0.1297016143798828 seconds
Processed 100 catalogs in 0.14351439476013184 seconds
Processed 200 catalogs in 0.2807955741882324 seconds
Processed 300 catalogs in 0.42336130142211914 seconds
Processed 400 catalogs in 0.5637788772583008 seconds
Processed 500 catalogs in 0.7492990493774414 seconds
Processed 600 catalogs in 0.8895461559295654 seconds
Processed 700 catalogs in 1.0328381061553955 seconds
Processed 800 catalogs in 1.2236864566802979 seconds
Processed 900 catalogs in 1.363717794418335 seconds
Processed 1000 catalogs in 1.5060994625091553 seconds
Processed 2000 catalogs in 3.139495372772217 seconds
Processed 3000 catalogs in 4.692880392074585 seconds
Processed 4000 catalogs in 6.291493654251099 seconds
Processed 5000 catalogs in 7.870847463607788 seconds
Processed 6000 catalogs in 9.423059701919556 seconds
Processed 7000 catalogs in 11.010499000549316 seconds
Processed 8000 catalogs in 12.49993085861206 seconds
Processed 9000 catalogs in 14.128604650497437 seconds
Processed 10000 catalogs in 15.664738655090332 seconds

Plot expected event counts

We can plot the expected event counts the same way that we plot a csep.core.forecasts.GriddedForecast

ax = forecast.expected_rates.plot(plot_args={'clim': [-3.5, 0]}, show=True)

The images holes in the image are due to under-sampling from the forecast.

Quick sanity check

The forecasts were filtered to the spatial region so all events should be binned. We loop through each data in the forecast and count the number of events and compare that with the expected rates. The expected rate is an average in each space-magnitude bin, so we have to multiply this value by the number of catalogs in the forecast.

Total running time of the script: ( 0 minutes 17.969 seconds)

Gallery generated by Sphinx-Gallery