Catalog-based Forecast Evaluation

This example shows how to evaluate a catalog-based forecasting using the Number test. This test is the simplest of the evaluations.

Overview:

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

2. Access catalog from ComCat

3. Filter catalog to be consistent with the forecast properties

4. Apply catalog-based number test to catalog

5. Visualize results for catalog-based forecast

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 csep
from csep.core import regions, catalog_evaluations
from csep.utils import datasets, time_utils

Define start and end times of forecast

Forecasts should define a time horizon in which they are valid. The choice is flexible for catalog-based forecasts, because the catalogs can be filtered to accommodate multiple end-times. Conceptually, these should be separate forecasts.

start_time = time_utils.strptime_to_utc_datetime("1992-06-28 11:57:35.0")
end_time = time_utils.strptime_to_utc_datetime("1992-07-28 11:57:35.0")

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. The forecast is already filtered in space and magnitude
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)
space_magnitude_region = regions.create_space_magnitude_region(region, magnitudes)

Load catalog forecast

To reduce the file size of this example, we’ve already filtered the catalogs to the appropriate magnitudes and spatial locations. The original forecast was computed for 1 year following the start date, so we still need to filter the catalog in time. We can do this by passing a list of filtering arguments to the forecast or updating the class.

By default, the forecast loads catalogs on-demand, so the filters are applied as the catalog loads. On-demand means that until we loop over the forecast in some capacity, none of the catalogs are actually loaded.

More fine-grain control and optimizations can be achieved by creating a csep.core.forecasts.CatalogForecast directly.

forecast = csep.load_catalog_forecast(datasets.ucerf3_ascii_format_landers_fname,
                                      start_time = start_time, end_time = end_time,
                                      region = space_magnitude_region,
                                      apply_filters = True)

# Assign filters to forecast
forecast.filters = [f'origin_time >= {forecast.start_epoch}', f'origin_time < {forecast.end_epoch}']

Obtain evaluation catalog from ComCat

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.

We need to filter the ComCat catalog to be consistent with the forecast. This can be done either through the ComCat API or using catalog filtering strings. Here we’ll use the ComCat API to make the data access quicker for this example. We still need to filter the observed catalog in space though.

# Obtain Comcat catalog and filter to region.
comcat_catalog = csep.query_comcat(start_time, end_time, min_magnitude=forecast.min_magnitude)

# Filter observed catalog using the same region as the forecast
comcat_catalog = comcat_catalog.filter_spatial(forecast.region)
print(comcat_catalog)

# Plot the catalog
comcat_catalog.plot()

Fetched ComCat catalog in 0.5070505142211914 seconds.

Downloaded catalog from ComCat with following parameters
Start Date: 1992-06-28 12:00:45+00:00
End Date: 1992-07-24 18:14:36.250000+00:00
Min Latitude: 33.901 and Max Latitude: 36.705
Min Longitude: -118.067 and Max Longitude: -116.285
Min Magnitude: 4.95
Found 19 events in the ComCat catalog.

        Name: None

        Start Date: 1992-06-28 12:00:45+00:00
        End Date: 1992-07-24 18:14:36.250000+00:00

        Latitude: (33.901, 36.705)
        Longitude: (-118.067, -116.285)

        Min Mw: 4.95
        Max Mw: 6.3

        Event Count: 19


<GeoAxes: >

Perform number test

We can perform the Number test on the catalog based forecast using the observed catalog we obtained from Comcat.

number_test_result = catalog_evaluations.number_test(forecast, comcat_catalog)

Processed 1 catalogs in 0.001165628433227539 seconds
Processed 2 catalogs in 0.0015769004821777344 seconds
Processed 3 catalogs in 0.001963376998901367 seconds
Processed 4 catalogs in 0.0022356510162353516 seconds
Processed 5 catalogs in 0.0024819374084472656 seconds
Processed 6 catalogs in 0.002842426300048828 seconds
Processed 7 catalogs in 0.0031113624572753906 seconds
Processed 8 catalogs in 0.003477334976196289 seconds
Processed 9 catalogs in 0.004408359527587891 seconds
Processed 10 catalogs in 0.0047495365142822266 seconds
Processed 20 catalogs in 0.008414506912231445 seconds
Processed 30 catalogs in 0.012235164642333984 seconds
Processed 40 catalogs in 0.015888452529907227 seconds
Processed 50 catalogs in 0.019929885864257812 seconds
Processed 60 catalogs in 0.023308753967285156 seconds
Processed 70 catalogs in 0.026613235473632812 seconds
Processed 80 catalogs in 0.030157089233398438 seconds
Processed 90 catalogs in 0.03392767906188965 seconds
Processed 100 catalogs in 0.03716707229614258 seconds
Processed 200 catalogs in 0.06937623023986816 seconds
Processed 300 catalogs in 0.10452413558959961 seconds
Processed 400 catalogs in 0.1388390064239502 seconds
Processed 500 catalogs in 0.20105242729187012 seconds
Processed 600 catalogs in 0.2340714931488037 seconds
Processed 700 catalogs in 0.2680518627166748 seconds
Processed 800 catalogs in 0.3326606750488281 seconds
Processed 900 catalogs in 0.3651440143585205 seconds
Processed 1000 catalogs in 0.3985569477081299 seconds
Processed 2000 catalogs in 0.8680195808410645 seconds
Processed 3000 catalogs in 1.2929298877716064 seconds
Processed 4000 catalogs in 1.7397739887237549 seconds
Processed 5000 catalogs in 2.174591302871704 seconds
Processed 6000 catalogs in 2.610468864440918 seconds
Processed 7000 catalogs in 3.0657706260681152 seconds
Processed 8000 catalogs in 3.4687514305114746 seconds
Processed 9000 catalogs in 3.9487407207489014 seconds
Processed 10000 catalogs in 4.3683555126190186 seconds

Plot number test result

We can create a simple visualization of the number test from the evaluation result class.

ax = number_test_result.plot(show=True)