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.002006053924560547 seconds
Processed 2 catalogs in 0.00373077392578125 seconds
Processed 3 catalogs in 0.0051877498626708984 seconds
Processed 4 catalogs in 0.006409168243408203 seconds
Processed 5 catalogs in 0.007552623748779297 seconds
Processed 6 catalogs in 0.008825540542602539 seconds
Processed 7 catalogs in 0.00997614860534668 seconds
Processed 8 catalogs in 0.011281728744506836 seconds
Processed 9 catalogs in 0.013590335845947266 seconds
Processed 10 catalogs in 0.01488947868347168 seconds
Processed 20 catalogs in 0.02682781219482422 seconds
Processed 30 catalogs in 0.03980064392089844 seconds
Processed 40 catalogs in 0.05248427391052246 seconds
Processed 50 catalogs in 0.06592440605163574 seconds
Processed 60 catalogs in 0.07828855514526367 seconds
Processed 70 catalogs in 0.0904853343963623 seconds
Processed 80 catalogs in 0.10318279266357422 seconds
Processed 90 catalogs in 0.11621451377868652 seconds
Processed 100 catalogs in 0.12835907936096191 seconds
Processed 200 catalogs in 0.2501952648162842 seconds
Processed 300 catalogs in 0.37784457206726074 seconds
Processed 400 catalogs in 0.5035185813903809 seconds
Processed 500 catalogs in 0.679222583770752 seconds
Processed 600 catalogs in 0.8030292987823486 seconds
Processed 700 catalogs in 0.9286062717437744 seconds
Processed 800 catalogs in 1.1094391345977783 seconds
Processed 900 catalogs in 1.2325737476348877 seconds
Processed 1000 catalogs in 1.357952356338501 seconds
Processed 2000 catalogs in 2.8450138568878174 seconds
Processed 3000 catalogs in 4.259924411773682 seconds
Processed 4000 catalogs in 5.715805530548096 seconds
Processed 5000 catalogs in 7.1464104652404785 seconds
Processed 6000 catalogs in 8.570996284484863 seconds
Processed 7000 catalogs in 10.042028665542603 seconds
Processed 8000 catalogs in 11.41623330116272 seconds
Processed 9000 catalogs in 12.9328773021698 seconds
Processed 10000 catalogs in 14.343766212463379 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]})


/opt/hostedtoolcache/Python/3.7.9/x64/lib/python3.7/site-packages/cartopy/mpl/ UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
  warnings.warn('The .ylabels_right attribute is deprecated. Please '

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 22.564 seconds)

Gallery generated by Sphinx-Gallery