Repository: SCECcode/pycsep (language: python)

Issue:
initial commit



Discussion / hints:

Fix patch:
diff --git a/csep/core/forecasts.py b/csep/core/forecasts.py
index 6fac74b3..323e1299 100644
--- a/csep/core/forecasts.py
+++ b/csep/core/forecasts.py
@@ -639,17 +639,15 @@ def min_magnitude(self):
 
     def spatial_counts(self, cartesian=False):
         """ Returns the expected spatial counts from forecast """
-        if self.expected_rates is not None:
-            return self.expected_rates.spatial_counts(cartesian=cartesian)
-        else:
-            return None
+        if self.expected_rates is None:
+            self.get_expected_rates()
+        return self.expected_rates.spatial_counts(cartesian=cartesian)
 
     def magnitude_counts(self):
         """ Returns expected magnitude counts from forecast """
-        if self.expected_rates is not None:
-            return self.expected_rates.magnitude_counts()
-        else:
-            return None
+        if self.expected_rates is None:
+            self.get_expected_rates()
+        return self.expected_rates.magnitude_counts()
 
     def get_event_counts(self):
         """ Returns a numpy array containing the number of event counts for each catalog.
diff --git a/csep/core/regions.py b/csep/core/regions.py
index bc26a9c2..26ccbd66 100644
--- a/csep/core/regions.py
+++ b/csep/core/regions.py
@@ -749,7 +749,20 @@ def tight_bbox(self):
         unique_poly = np.append(unique_poly, [unique_poly[0, :]], axis=0)
         return unique_poly
 
-def geographical_area_from_bounds(lon1,lat1,lon2,lat2):
+    def get_cell_area(self):
+        """ Compute the area of each polygon in sq. kilometers.
+
+            Returns:
+                out (numpy.array): numpy array containing cell area in km^2
+        """
+        area = numpy.zeros(self.num_nodes)
+        for idx, origin in enumerate(self.origins()):
+            top_right = origin + self.dh
+            area[idx] = geographical_area_from_bounds(origin[0], origin[1], top_right[0], top_right[1])
+        return area
+
+
+def geographical_area_from_bounds(lon1, lat1, lon2, lat2):
     """
     Computes area of spatial cell identified by origin coordinate and top right cooridnate.
     The functions computes area only for square/rectangle bounding box by based on spherical earth assumption.
diff --git a/csep/utils/comcat.py b/csep/utils/comcat.py
index 19c27305..da0d2558 100644
--- a/csep/utils/comcat.py
+++ b/csep/utils/comcat.py
@@ -21,11 +21,10 @@
 # PyCSEP imports
 from csep.utils.time_utils import HistoricTime
 
-# constants
 # url template for counting events
 HOST = 'earthquake.usgs.gov'
 SEARCH_TEMPLATE = 'https://[HOST]/fdsnws/event/1/query?format=geojson'
-TIMEOUT = 60  # how long do we wait for a url to return?
+TIMEOUT = 120  # how long do we wait for a url to return?
 TIMEFMT = '%Y-%m-%dT%H:%M:%S'
 WAITSECS = 3  # number of seconds to wait after failing download before trying again
 SEARCH_LIMIT = 20000  # maximum number of events ComCat will return in one search
diff --git a/csep/utils/constants.py b/csep/utils/constants.py
index 13f401b9..adbe780b 100644
--- a/csep/utils/constants.py
+++ b/csep/utils/constants.py
@@ -10,6 +10,8 @@
 
 MW_5_EQS_PER_YEAR = 10
 
+earth_radius_km = 6371.
+
 # Magnitude Bins
 min_mw = 2.5
 max_mw = 8.95
diff --git a/csep/utils/plots.py b/csep/utils/plots.py
index f4308202..6afcfc16 100644
--- a/csep/utils/plots.py
+++ b/csep/utils/plots.py
@@ -17,6 +17,7 @@
 from csep.utils.constants import SECONDS_PER_DAY, CSEP_MW_BINS
 from csep.utils.calc import bin1d_vec
 from csep.utils.time_utils import datetime_to_utc_epoch
+from csep.core import regions
 
 
 """
@@ -1727,6 +1728,146 @@ def plot_calibration_test(evaluation_result, axes=None, plot_args=None, show=Fal
 
     return ax
 
+def plot_ROC(forecast, catalog, axes=None, plot_uniform=True, savepdf=True, savepng=True, show=True,
+                    plot_args=None):
+    """
+    Plot Receiver operating characteristic (ROC) Curves based on forecast and test catalog.
+
+    The ROC is computed following this procedure:
+        (1) Obtain spatial rates from GriddedForecast
+        (2) Rank the rates in descending order (highest rates first).
+        (3) Sort forecasted rates by ordering found in (2), and normalize rates so the cumulative sum equals unity.
+        (4) Obtain binned spatial rates from observed catalog
+        (5) Sort gridded observed rates by ordering found in (2), and normalize so the cumulative sum equals unity.
+        (6) Compute spatial bin areas, sort by ordering found in (2), and normalize so the cumulative sum equals unity.
+        (7) Plot ordered and cumulated rates (forecast and catalog) against ordered and cumulated bin areas.
+
+    Note that: 
+        (1) The testing catalog and forecast should have exactly the same time-window (duration)
+        (2) Forecasts should be defined over the same region
+        (3) If calling this function multiple times, update the color in plot_args
+
+    Args:
+        forecast (:class: `csep.forecast.GriddedForecast`):
+        catalog (:class:`AbstractBaseCatalog`): evaluation catalog
+        axes (:class:`matplotlib.pyplot.ax`): Previously defined ax object
+        savepdf (str):    output option of pdf file
+        savepng (str):    output option of png file
+        plot_uniform (bool): if true, include uniform forecast on plot
+
+        Optional plotting arguments:
+            * figsize: (:class:`list`/:class:`tuple`) - default: [9, 8]
+            * forecast_linecolor: (:class:`str`) - default: 'black'
+            * title_fontsize: (:class:`float`) Fontsize of the plot title - default: 18
+            * forecast_linecolor: (:class:`str`) - default: 'black'
+            * forecast_linestyle: (:class:`str`) - default: '-'
+            * observed_linecolor: (:class:`str`) - default: 'blue'
+            * observed_linestyle: (:class:`str`) - default: '-'
+            * forecast_label: (:class:`str`) - default: Observed (Forecast)
+            * legend_fontsize: (:class:`float`) Fontsize of the plot title - default: 16
+            * legend_loc: (:class:`str`) - default: 'upper left'
+            * title_fontsize: (:class:`float`) Fontsize of the plot title - default: 18
+            * label_fontsize: (:class:`float`) Fontsize of the plot title - default: 14
+            * title: (:class:`str`) - default: 'ROC Curve'
+            * filename: (:class:`str`) - default: roc_curve.
+
+    Returns:
+        :class:`matplotlib.pyplot.ax` object
+
+    Raises:
+        TypeError: throws error if CatalogForecast-like object is provided
+        RuntimeError: throws error if Catalog and Forecast do not have the same region
+
+        Written by Han Bao, UCLA, March 2021. Modified June 2021.
+    """
+    if not catalog.region == forecast.region:
+        raise RuntimeError("catalog region and forecast region must be identical.")
+
+    # Parse plotting arguments
+    plot_args = plot_args or {}
+    figsize = plot_args.get('figsize', (9, 8))
+    forecast_linecolor = plot_args.get('forecast_linecolor', 'black')
+    forecast_linestyle = plot_args.get('forecast_linestyle', '-')
+    observed_linecolor = plot_args.get('observed_linecolor', 'blue')
+    observed_linestyle = plot_args.get('observed_linestyle', '-')
+    legend_fontsize = plot_args.get('legend_fontsize', 16)
+    legend_loc = plot_args.get('legend_loc', 'upper left')
+    title_fontsize = plot_args.get('title_fontsize', 18)
+    label_fontsize = plot_args.get('label_fontsize', 14)
+    filename = plot_args.get('filename', 'roc_figure')
+    title = plot_args.get('title', 'ROC Curve')
+
+    # Plot catalog ordered by forecast rates
+    name = forecast.name
+    if not name:
+        name = ''
+    else:
+        name = f'({name})'
+
+    forecast_label = plot_args.get('forecast_label', f'Forecast {name}')
+    observed_label = plot_args.get('observed_label', f'Observed {name}')
+
+    # Initialize figure
+    if axes is not None:
+        ax = axes
+    else:
+        fig, ax = pyplot.subplots(figsize=figsize)
+
+    # This part could be vectorized if optimizations become necessary
+    # Initialize array to store cell area in km^2
+    area_km2 = catalog.region.get_cell_area()
+    obs_counts = catalog.spatial_counts()
+
+    # Obtain rates (or counts) aggregated in spatial cells
+    # If CatalogForecast, needs expected rates. Might take some time to compute.
+    rate = forecast.spatial_counts()
+
+    # Get index of rates (descending sort)
+    I = numpy.argsort(rate)
+    I = numpy.flip(I)
+
+    # Order forecast and cells rates by highest rate cells first
+    fore_norm_sorted = numpy.cumsum(rate[I]) / numpy.sum(rate)
+    area_norm_sorted = numpy.cumsum(area_km2[I]) / numpy.sum(area_km2)
+
+    # Compute normalized and sorted rates of observations
+    obs_norm_sorted = numpy.cumsum(obs_counts[I]) / numpy.sum(obs_counts)
+
+    # Plot uniform forecast
+    if plot_uniform:
+        ax.plot(area_norm_sorted, area_norm_sorted, 'k--', label='Uniform')
+
+    # Plot sorted and normalized forecast (descending order)
+    ax.plot(area_norm_sorted, fore_norm_sorted,
+            label=forecast_label,
+            color=forecast_linecolor,
+            linestyle = forecast_linestyle)
+
+    # Plot cell-wise rates of observed catalog ordered by forecast rates (descending order)
+    ax.step(area_norm_sorted, obs_norm_sorted,
+            label=observed_label,
+            color=observed_linecolor,
+            linestyle = observed_linestyle)
+
+    # Plotting arguments
+    ax.set_ylabel("True Positive Rate", fontsize=label_fontsize)
+    ax.set_xlabel('False Positive Rate (Normalized Area)', fontsize=label_fontsize)
+    ax.set_xscale('log')
+    ax.legend(loc=legend_loc, shadow=True, fontsize=legend_fontsize)
+    ax.set_title(title, fontsize=title_fontsize)
+
+    if filename:
+        if savepdf:
+            outFile = "{}.pdf".format(filename)
+            pyplot.savefig(outFile, format='pdf')
+        if savepng:
+            outFile = "{}.png".format(filename)
+            pyplot.savefig(outFile,format='png')
+
+    if show:
+        pyplot.show()
+    return ax
+    
 def add_labels_for_publication(figure, style='bssa', labelsize=16):
     """ Adds publication labels too the outside of a figure. """
     all_axes = figure.get_axes()
@@ -1739,4 +1880,4 @@ def add_labels_for_publication(figure, style='bssa', labelsize=16):
         if style == 'bssa':
             ax.annotate(f'({annot})', (0.025, 1.025), xycoords='axes fraction', fontsize=labelsize)
 
-    return
\ No newline at end of file
+    return
diff --git a/csep/utils/stats.py b/csep/utils/stats.py
index 3e135ceb..3550b0de 100644
--- a/csep/utils/stats.py
+++ b/csep/utils/stats.py
@@ -1,6 +1,8 @@
 import numpy
 import scipy.stats
 import scipy.special
+# PyCSEP imports
+from csep.core import regions
 
 def sup_dist(cdf1, cdf2):
     """
@@ -205,3 +207,65 @@ def poisson_inverse_cdf(random_matrix, lam):
         sample from the poisson distribution
     """
     return scipy.stats.poisson.ppf(random_matrix, lam)
+
+def get_Kagan_I1_score(forecasts, catalog):
+    """
+    A program for scoring (I_1) earthquake-forecast grids by the methods of:
+    Kagan, Yan Y. [2009] Testing long-term earthquake forecasts: likelihood methods
+                         and error diagrams, Geophys. J. Int., v. 177, pages 532-542.
+
+    Some advantages of these methods are that they:
+        - are insensitive to the grid used to cover the Earth;
+        - are insensitive to changes in the overall seismicity rate;
+        - do not modify the locations or magnitudes of test earthquakes;
+        - do not require simulation of virtual catalogs;
+        - return relative quality measures, not just "pass" or "fail;" and
+        - indicate relative specificity of forecasts as well as relative success.
+    
+    Written by Han Bao, UCLA, March 2021. Modified June 2021
+    
+    Note that: 
+        (1) The testing catalog and forecast should have exactly the same time-window (duration)
+        (2) Forecasts and catalogs have identical regions
+
+    Args:
+        forecasts:  csep.forecast or a list of csep.forecast (one catalog to test against different forecasts)
+        catalog:    csep.catalog 
+        
+    Returns:
+       I_1 (numpy.array): containing I1 for each forecast in inputs
+
+    """
+    ### Determine if input 'forecasts' is a list of csep.forecasts or a single csep.forecasts
+    try:
+        n_forecast = len(forecasts) # the input forecasts is a list of csep.forecast
+    except:  
+        n_forecast = 1             # the input forecasts is a single csep.forecast
+        forecasts = [forecasts]
+
+    # Sanity checks can go here
+    for forecast in forecasts:
+        if forecast.region != catalog.region:
+            raise RuntimeError("Catalog and forecasts must have identical regions.")
+
+    # Initialize array
+    I_1 = numpy.zeros(n_forecast, dtype=numpy.float64)
+
+    # Compute cell areas
+    area_km2 = catalog.region.get_cell_area()
+    total_area = numpy.sum(area_km2)
+    
+    for j, forecast in enumerate(forecasts):
+        # Eq per cell per duration in forecast; note, if called on a CatalogForecast this could require computed expeted rates
+        rate = forecast.spatial_counts()
+        # Get Rate Density and uniform_forecast of the Forecast
+        rate_den = rate / area_km2
+        uniform_forecast = numpy.sum(rate) / total_area
+        # Compute I_1 score
+        n_event = catalog.event_count
+        counts = catalog.spatial_counts()
+        non_zero_idx = numpy.argwhere(rate_den > 0)
+        non_zero_idx = non_zero_idx[:,0]
+        I_1[j] = numpy.dot(counts[non_zero_idx], numpy.log2(rate_den[non_zero_idx] / uniform_forecast)) / n_event
+    
+    return I_1
diff --git a/examples/tutorials/catalog_forecast_evaluation.py b/examples/tutorials/catalog_forecast_evaluation.py
index 6fcee589..c8f022a8 100644
--- a/examples/tutorials/catalog_forecast_evaluation.py
+++ b/examples/tutorials/catalog_forecast_evaluation.py
@@ -24,7 +24,7 @@
 
 import csep
 from csep.core import regions, catalog_evaluations
-from csep.utils import datasets, time_utils
+from csep.utils import datasets, time_utils, plots
 
 ####################################################################################################################################
 # Define start and end times of forecast
@@ -86,14 +86,13 @@
 # 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
+# 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)
 
@@ -115,3 +114,20 @@
 # We can create a simple visualization of the number test from the evaluation result class.
 
 ax = number_test_result.plot()
+
+####################################################################################################################################
+# Plot ROC Curves
+# -----------------------
+#
+# We can also plot the Receiver operating characteristic (ROC) Curves based on forecast and testing-catalog.
+# In the figure below, False Positive Rate is the normalized cumulative forecast rate, after sorting cells in decreasing order of rate.
+# The "True Positive Rate" is the normalized cumulative area. The dashed line is the ROC curver of 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 an 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.
+
+print("Plotting ROC curve")
+plots.plot_ROC(forecast, comcat_catalog)
\ No newline at end of file
diff --git a/examples/tutorials/gridded_forecast_evaluation.py b/examples/tutorials/gridded_forecast_evaluation.py
index daab0e46..5b3bf647 100644
--- a/examples/tutorials/gridded_forecast_evaluation.py
+++ b/examples/tutorials/gridded_forecast_evaluation.py
@@ -27,6 +27,7 @@
 import csep
 from csep.core import poisson_evaluations as poisson
 from csep.utils import datasets, time_utils, plots
+from csep.utils.stats import get_Kagan_I1_score
 
 ####################################################################################################################################
 # Define forecast properties
@@ -59,8 +60,8 @@
 # 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.
 
-catalog = csep.query_comcat(forecast.start_time, forecast.end_time,
-                            min_magnitude=forecast.min_magnitude)
+print("Querying comcat catalog")
+catalog = csep.query_comcat(forecast.start_time, forecast.end_time, min_magnitude=forecast.min_magnitude)
 print(catalog)
 
 ####################################################################################################################################
@@ -100,3 +101,21 @@
 
 ax = plots.plot_poisson_consistency_test(spatial_test_result,
                                         plot_args={'xlabel': 'Spatial likelihood'})
+
+####################################################################################################################################
+# Calculate Kagan's I_1 score 
+# -----------------------
+#
+# We can also get the Kagan's I_1 score (see Kagan, Yan Y. [2009] Testing long-term earthquake forecasts: likelihood methods and error diagrams, Geophys. J. Int., v. 177, pages 532-542.
+
+I_1 = get_Kagan_I1_score(forecast, catalog)
+print("I_1 score is: ", I_1)
+
+####################################################################################################################################
+# Plot ROC Curves
+# -----------------------
+#
+# We can also plot the Receiver operating characteristic (ROC) Curves based on forecast and testing-catalog.
+# In the figure below, False Positive Rate is the normalized cumulative forecast rate, after sorting cells in decreasing order of rate. The "True Positive Rate" is the normalized cumulative area. The dashed line is the ROC curver of 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 an catalog, one can evaluate if the forecast is more or less specific (or smooth) at different level or seismic rate.
+
+plots.plot_ROC(forecast, catalog)


Test patch:
diff --git a/tests/test_regions.py b/tests/test_regions.py
index 5c9ad44c..ac5921ec 100644
--- a/tests/test_regions.py
+++ b/tests/test_regions.py
@@ -43,17 +43,19 @@ def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
     def setUp(self):
+        self.r = california_relm_region()
         self.from_dat = numpy.loadtxt(get_california_region_fname())
         self.num_nodes = len(self.from_dat)
 
     def test_node_count(self):
         """ Ensures the node counts are consistent between the two files. """
-        r = california_relm_region()
-        self.assertEqual(self.num_nodes, r.num_nodes)
+        self.assertEqual(self.num_nodes, self.r.num_nodes)
 
     def test_origins(self):
         """ Compares XML file against the simple .dat file containing the region. """
-        r = california_relm_region()
         # they dont have to be in the same order, but they need
-        numpy.testing.assert_array_equal(r.midpoints().sort(), self.from_dat.sort())
+        numpy.testing.assert_array_equal(self.r.midpoints().sort(), self.from_dat.sort())
 
+    def test_eq_oper(self):
+        r = california_relm_region()
+        assert self.r == r
\ No newline at end of file