build: added rasterio to conda requirements

docs: Added "plot composition" section to  Plotting tutorial. Added plot_magnitude_vs_time and plot_magnitude_histogram examples to Catalog-forecast tutorial. Added reflinks to tutorial headers. Added load_evaluation_result to API ref.
ft: Added default plot_region=True to forecast.plot(). Added log_scale option to plot_magnitude_histogram for y axis, and also the mean is plotted now instead of median. Remove tight_layout from plot_catalog() option.
fix: added transform argument when plotting the region in plot_gridded_dataset.

Author: pabloitu

csep/core/catalogs.py

@@ -857,7 +857,7 @@ def plot(self, ax=None, show=False, extent=None, set_global=False, **kwargs):
 
         # no mutable function arguments
         plot_args_default = {
-             'basemap': kwargs.pop('basemap', None) or 'ESRI_terrain',
+             'basemap': kwargs.pop('basemap', 'ESRI_terrain') if ax is None else None
         }
 
         # Plot the region border (if it exists) by default
@@ -870,10 +870,11 @@ def plot(self, ax=None, show=False, extent=None, set_global=False, **kwargs):
 
         plot_args = kwargs.get('plot_args', {})
         plot_args_default.update(plot_args)
+        plot_args_default.update(kwargs)
 
         # this call requires internet connection and basemap
         ax = plot_catalog(self, ax=ax, show=show, extent=extent,
-                          set_global=set_global, **plot_args_default, **kwargs)
+                          set_global=set_global, **plot_args_default)
         return ax
 
     def plot_magnitude_versus_time(self, ax=None, show=False, **kwargs):

csep/core/forecasts.py

@@ -453,9 +453,10 @@ def plot(self, ax=None, show=False, log=True, extent=None, set_global=False, plo
 
         plot_args = plot_args or {}
         plot_args.update({
-             'basemap': kwargs.pop('basemap', None) or 'ESRI_terrain',
+             'basemap': kwargs.pop('basemap', 'ESRI_terrain') if ax is None else None,
              'title': kwargs.pop('title', None) or self.name,
              'figsize': kwargs.pop('figsize', None) or (8, 8),
+             'plot_region': True
         })
         plot_args.update(**kwargs)
         # this call requires internet connection and basemap

csep/utils/plots.py

@@ -393,6 +393,7 @@ def plot_magnitude_histogram(
     observation: "CSEPCatalog",
     magnitude_bins: Optional[Union[List[float], numpy.ndarray]] = None,
     percentile: int = 95,
+    log_scale: bool = True,
     ax: Optional["matplotlib.axes.Axes"] = None,
     show: bool = False,
     **kwargs,
@@ -411,6 +412,8 @@ def plot_magnitude_histogram(
             Defaults to `None`.
         percentile (int, optional): The percentile used for uncertainty intervals. Defaults to
             `95`.
+        log_scale (bool, optional): Whether to plot the y-axis in logarithmic scale. Defaults to
+            True.
         ax (matplotlib.axes.Axes, optional): The axes object to draw the plot on. If `None`, a
             new figure and axes are created. Defaults to `None`.
         show (bool, optional): Whether to display the plot immediately. Defaults to `False`.
@@ -482,6 +485,7 @@ def get_histogram_synthetic_cat(x, mags, normed=True):
     bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2
 
     # Compute statistics for the forecast histograms
+    forecast_mean = numpy.mean(forecast_hist, axis=0)
     forecast_median = numpy.median(forecast_hist, axis=0)
     forecast_low = numpy.percentile(forecast_hist, (100 - percentile) / 2.0, axis=0)
     forecast_high = numpy.percentile(forecast_hist, 100 - (100 - percentile) / 2.0, axis=0)
@@ -490,7 +494,12 @@ def get_histogram_synthetic_cat(x, mags, normed=True):
     forecast_err_upper = forecast_high - forecast_median
 
     # Plot observed histogram
-    ax.semilogy(
+    if log_scale:
+        plot_func = ax.semilogy
+    else:
+        plot_func = ax.plot
+
+    plot_func(
         bin_centers,
         obs_hist,
         color=plot_args["color"],
@@ -504,11 +513,11 @@ def get_histogram_synthetic_cat(x, mags, normed=True):
     # Plot forecast histograms as bar plot with error bars
     ax.plot(
         bin_centers,
-        forecast_median,
+        forecast_mean,
         ".",
         markersize=plot_args["markersize"],
         color="darkred",
-        label="Forecast Median",
+        label="Forecast Mean",
     )
     ax.errorbar(
         bin_centers,
@@ -815,9 +824,6 @@ def plot_catalog(
 
     ax.set_title(plot_args["title"], fontsize=plot_args["title_fontsize"], y=1.06)
 
-    if plot_args["tight_layout"]:
-        ax.figure.tight_layout()
-
     if show:
         pyplot.show()
 
@@ -934,7 +940,8 @@ def plot_gridded_dataset(
     if plot_region and not set_global:
         try:
             pts = region.tight_bbox()
-            ax.plot(pts[:, 0], pts[:, 1], lw=1, color=plot_args["region_color"])
+            ax.plot(pts[:, 0], pts[:, 1], lw=1, color=plot_args["region_color"],
+                    transform=ccrs.PlateCarree())
         except AttributeError:
             pass
 

docs/reference/api_reference.rst

@@ -7,20 +7,21 @@ This contains a reference document to the PyCSEP API.
 
 .. :currentmodule:: csep
 
-Loading catalogs and forecasts
-------------------------------
+Loading catalogs, forecasts and results
+---------------------------------------
 
 .. autosummary::
    :toctree: generated
 
-   load_stochastic_event_sets
    load_catalog
    query_comcat
    query_bsi
    query_gns
    query_gcmt
    load_gridded_forecast
    load_catalog_forecast
+   load_stochastic_event_sets
+   load_evaluation_result
 
 Catalogs
 --------

examples/tutorials/catalog_filtering.py

@@ -1,7 +1,7 @@
 """
 .. _tutorial-catalog-filtering:
 
-Catalogs operations
+Catalogs Operations
 ===================
 
 This example demonstrates how to perform standard operations on a catalog. This example requires an internet

examples/tutorials/catalog_forecast_evaluation.py

@@ -1,5 +1,6 @@
 """
-.. _catalog-forecast-evaluation:
+
+.. _tutorial-catalog-forecast-evaluation:
 
 Catalog-based Forecast Evaluation
 =================================
@@ -24,7 +25,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
@@ -71,12 +72,13 @@
 # More fine-grain control and optimizations can be achieved by creating a :class:`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)
+                                      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}']
+forecast.filters = [f'origin_time >= {forecast.start_epoch}',
+                    f'origin_time < {forecast.end_epoch}']
 
 ####################################################################################################################################
 # Obtain evaluation catalog from ComCat
@@ -99,6 +101,22 @@
 # Plot the catalog
 comcat_catalog.plot(show=True)
 
+####################################################################################################################################
+# Exploratory visualizations
+# --------------------------
+#
+# We can visualize differences between the distribution of event counts of the observations against the forecast:
+
+plots.plot_cumulative_events_versus_time(forecast, comcat_catalog, time_axis='days', show=True)
+
+
+####################################################################################################################################
+#
+# and their magnitude histograms
+#
+
+plots.plot_magnitude_histogram(forecast, comcat_catalog, log_scale=False, show=True)
+
 ####################################################################################################################################
 # Perform number test
 # -------------------
@@ -113,4 +131,4 @@
 #
 # We can create a simple visualization of the number test from the evaluation result class.
 
-ax = number_test_result.plot(show=True)
+ax = number_test_result.plot(show=True)

examples/tutorials/gridded_forecast_evaluation.py

@@ -28,17 +28,14 @@
 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
-
 ####################################################################################################################################
 # 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.
-from csep.utils.stats import get_Kagan_I1_score
+# because they can be rescaled to any arbitrary time period.
+
 
 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')
@@ -103,8 +100,8 @@
 # consistency tests.
 
 ax = plots.plot_consistency_test(spatial_test_result,
-                                 xlabel='Spatial likelihood')
-plt.show()
+                                 xlabel='Spatial likelihood',
+                                 show=True)
 
 
 ####################################################################################################################################
@@ -182,5 +179,6 @@
 # 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).
 
+from csep.utils.stats import get_Kagan_I1_score
 I_1 = get_Kagan_I1_score(forecast, catalog)
 print("I_1score is: ", I_1)