Merge pull request #263 from pabloitu/260-plot-refactoring

Refactored Plot Module, updating plotting functions and changed API to `csep.plots`

Author: pabloitu

.coveragerc

@@ -5,7 +5,6 @@ omit =
     tests/*
     docs/*
     examples/*
-    csep/utils/plots.py
     csep/utils/constants.py
     csep/utils/datasets.py
     csep/utils/documents.py

csep/__init__.py

@@ -3,23 +3,17 @@
 import time
 
 from csep._version import __version__
-
-from csep.core import forecasts
+from csep.core import catalog_evaluations
 from csep.core import catalogs
+from csep.core import forecasts
 from csep.core import poisson_evaluations
-from csep.core import catalog_evaluations
 from csep.core import regions
+from csep.core.exceptions import CSEPCatalogException
+from csep.core.forecasts import GriddedForecast, CatalogForecast
 from csep.core.repositories import (
     load_json,
     write_json
 )
-
-from csep.core.exceptions import CSEPCatalogException
-
-from csep.utils import datasets
-from csep.utils import readers
-
-from csep.core.forecasts import GriddedForecast, CatalogForecast
 from csep.models import (
     EvaluationResult,
     CatalogNumberTestResult,
@@ -28,7 +22,8 @@
     CatalogPseudolikelihoodTestResult,
     CalibrationTestResult
 )
-
+from csep.utils import datasets
+from csep.utils import readers
 from csep.utils.time_utils import (
     utc_now_datetime,
     strptime_to_utc_datetime,
@@ -333,7 +328,7 @@ def query_gns(start_time, end_time,  min_magnitude=2.950,
         verbose (bool): print catalog summary statistics
 
     Returns:
-        :class:`csep.core.catalogs.CSEPCatalog
+        :class:`csep.core.catalogs.CSEPCatalog`
     """
 
     # Timezone should be in UTC
@@ -366,7 +361,7 @@ def query_gcmt(start_time, end_time, min_magnitude=5.0,
                max_depth=None,
                catalog_id=None,
                min_latitude=None, max_latitude=None,
-               min_longitude=None, max_longitude=None):
+               min_longitude=None, max_longitude=None, verbose=True):
 
     eventlist = readers._query_gcmt(start_time=start_time,
                                      end_time=end_time,
@@ -381,6 +376,17 @@ def query_gcmt(start_time, end_time, min_magnitude=5.0,
                                    name='gCMT',
                                    catalog_id=catalog_id,
                                    date_accessed=utc_now_datetime())
+
+    if verbose:
+        print("Downloaded catalog from GCMT with following parameters")
+        print("Start Date: {}\nEnd Date: {}".format(str(catalog.start_time),
+                                                    str(catalog.end_time)))
+        print("Min Latitude: {} and Max Latitude: {}".format(catalog.min_latitude,
+                                                             catalog.max_latitude))
+        print("Min Longitude: {} and Max Longitude: {}".format(catalog.min_longitude,
+                                                               catalog.max_longitude))
+        print("Min Magnitude: {}".format(catalog.min_magnitude))
+        print(f"Found {catalog.event_count} events in the gns catalog.")
     return catalog
 
 

csep/artifacts/ExampleBasemaps/basemap_california.tif

@@ -65,7 +65,7 @@ def negative_binomial_number_test(gridded_forecast, observed_catalog, variance):
     delta1, delta2 = _nbd_number_test_ndarray(fore_cnt, obs_cnt, variance, epsilon=epsilon)
 
     # store results
-    result.test_distribution = ('negative_binomial', fore_cnt)
+    result.test_distribution = ('negative_binomial', fore_cnt, variance)
     result.name = 'NBD N-Test'
     result.observed_statistic = obs_cnt
     result.quantile = (delta1, delta2)
@@ -126,7 +126,7 @@ def _simulate_catalog(sim_cells, sampling_weights, sim_fore, random_numbers=None
 
 
 def _binary_likelihood_test(forecast_data, observed_data, num_simulations=1000, random_numbers=None, 
-                              seed=None, use_observed_counts=True, verbose=True, normalize_likelihood=False):
+                              seed=None, use_observed_counts=True, verbose=False, normalize_likelihood=False):
     """  Computes binary conditional-likelihood test from CSEP using an efficient simulation based approach.
     
     Args:

csep/core/binomial_evaluations.py

@@ -71,7 +71,7 @@ def _simulate_catalog(sim_cells, sampling_weights, random_numbers=None):
 
 
 def _brier_score_test(forecast_data, observed_data, num_simulations=1000,
-                      random_numbers=None, seed=None, verbose=True):
+                      random_numbers=None, seed=None, verbose=False):
     """  Computes the Brier consistency test conditional on the total observed
     number of events
 

csep/core/brier_evaluations.py

@@ -23,7 +23,7 @@
     from csep.core.forecasts import CatalogForecast
 
 
-def number_test(forecast, observed_catalog, verbose=True):
+def number_test(forecast, observed_catalog, verbose=False):
     """ Performs the number test on a catalog-based forecast.
 
     The number test builds an empirical distribution of the event counts for each data. By default, this
@@ -62,7 +62,7 @@ def number_test(forecast, observed_catalog, verbose=True):
     return result
 
 
-def spatial_test(forecast, observed_catalog, verbose=True):
+def spatial_test(forecast, observed_catalog, verbose=False):
     """ Performs spatial test for catalog-based forecasts.
 
 
@@ -136,7 +136,7 @@ def spatial_test(forecast, observed_catalog, verbose=True):
         delta_1, delta_2 = get_quantiles(test_distribution_spatial_1d, obs_lh_norm)
 
     result = CatalogSpatialTestResult(test_distribution=test_distribution_spatial_1d,
-                                      name='S-Test',
+                                      name='Catalog S-Test',
                                       observed_statistic=obs_lh_norm,
                                       quantile=(delta_1, delta_2),
                                       status=message,
@@ -148,7 +148,7 @@ def spatial_test(forecast, observed_catalog, verbose=True):
     return result
 
 
-def magnitude_test(forecast, observed_catalog, verbose=True):
+def magnitude_test(forecast, observed_catalog, verbose=False):
     """ Performs magnitude test for catalog-based forecasts """
     test_distribution = []
 
@@ -160,7 +160,7 @@ def magnitude_test(forecast, observed_catalog, verbose=True):
         print("Cannot perform magnitude test when observed event count is zero.")
         # prepare result
         result = CatalogMagnitudeTestResult(test_distribution=test_distribution,
-                                            name='M-Test',
+                                            name='Catalog M-Test',
                                             observed_statistic=None,
                                             quantile=(None, None),
                                             status='not-valid',
@@ -212,7 +212,7 @@ def magnitude_test(forecast, observed_catalog, verbose=True):
 
     # prepare result
     result = CatalogMagnitudeTestResult(test_distribution=test_distribution,
-                              name='M-Test',
+                              name='Catalog M-Test',
                               observed_statistic=obs_d_statistic,
                               quantile=(delta_1, delta_2),
                               status='normal',
@@ -224,7 +224,7 @@ def magnitude_test(forecast, observed_catalog, verbose=True):
     return result
 
 
-def pseudolikelihood_test(forecast, observed_catalog, verbose=True):
+def pseudolikelihood_test(forecast, observed_catalog, verbose=False):
     """ Performs the spatial pseudolikelihood test for catalog forecasts.
 
     Performs the spatial pseudolikelihood test as described by Savran et al., 2020. The tests uses a pseudolikelihood
@@ -307,7 +307,7 @@ def pseudolikelihood_test(forecast, observed_catalog, verbose=True):
     # prepare evaluation result
     result = CatalogPseudolikelihoodTestResult(
         test_distribution=test_distribution_1d,
-        name='PL-Test',
+        name='Catalog PL-Test',
         observed_statistic=obs_plh,
         quantile=(delta_1, delta_2),
         status=message,

csep/core/catalog_evaluations.py

@@ -23,7 +23,7 @@
 from csep.utils.log import LoggingMixin
 from csep.utils.readers import csep_ascii
 from csep.utils.file import get_file_extension
-from csep.utils.plots import plot_catalog
+from csep.plots import plot_catalog, plot_magnitude_versus_time
 
 
 class AbstractBaseCatalog(LoggingMixin):
@@ -840,48 +840,59 @@ def b_positive(self):
         """ Implements the b-positive indicator from Nicholas van der Elst """
         pass
 
-    def plot(self, ax=None, show=False, extent=None, set_global=False, plot_args=None):
-        """ Plot catalog according to plate-carree projection
+    def plot(self, ax=None, show=False, extent=None, set_global=False, **kwargs):
+        """ Plots the catalog epicenters.
+
+        See :func:`csep.utils.plots.plot_catalog` for a description of keyword arguments.
 
         Args:
-            ax (`matplotlib.pyplot.axes`): Previous axes onto which catalog can be drawn
-            show (bool): if true, show the figure. this call is blocking.
+            ax (matplotlib.pyplot.axes): Previous axes onto which catalog can be drawn
+            show (bool): If True, shows the figure.
             extent (list): Force an extent [lon_min, lon_max, lat_min, lat_max]
-            plot_args (optional/dict): dictionary containing plotting arguments for making figures
+            set_global (bool): Whether to plot using a global projection
+            **kwargs (dict): Keyword arguments passed to
+                :func:`csep.utils.plots.plot_catalog`
 
         Returns:
             axes: matplotlib.Axes.axes
         """
 
         # no mutable function arguments
-        plot_args_default = {
-             'basemap': 'ESRI_terrain',
-             'markersize': 2,
-             'markercolor': 'red',
-             'alpha': 0.3,
-             'mag_scale': 7,
-             'legend': True,
-             'grid_labels': True,
-             'legend_loc': 3,
-             'figsize': (8, 8),
-             'title': self.name,
-             'mag_ticks': False
+        plot_args = {
+             'basemap': kwargs.pop('basemap', 'ESRI_terrain') if ax is None else None
         }
 
         # Plot the region border (if it exists) by default
         try:
             # This will throw error if catalog does not have region
             _ = self.region.num_nodes
-            plot_args_default['region_border'] = True
+            plot_args['region_border'] = True
         except AttributeError:
             pass
 
-        plot_args = plot_args or {}
-        plot_args_default.update(plot_args)
+        plot_args.update(kwargs.get('plot_args', {}))
+        plot_args.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=plot_args_default)
+        ax = plot_catalog(self, ax=ax, show=show, extent=extent,
+                          set_global=set_global, **plot_args)
+        return ax
+
+    def plot_magnitude_versus_time(self, ax=None, show=False, **kwargs):
+        """ Plot the magnitude-time series of a catalog. See
+        https://docs.cseptesting.org/reference/generated/csep.utils.plots.plot_magnitude_versus_time.html for
+        a description of keyword arguments.
+
+        Args:
+            ax (`matplotlib.pyplot.axes`): Previous axes onto which catalog can be drawn
+            show (bool): if true, show the figure. this call is blocking.
+
+        Returns:
+            axes: matplotlib.Axes.axes
+        """
+
+        ax = plot_magnitude_versus_time(self, ax=ax, show=show, **kwargs)
+
         return ax
 
 

csep/core/catalogs.py

@@ -13,7 +13,7 @@
 from csep.utils.time_utils import decimal_year, datetime_to_utc_epoch
 from csep.core.catalogs import AbstractBaseCatalog
 from csep.utils.constants import SECONDS_PER_ASTRONOMICAL_YEAR
-from csep.utils.plots import plot_spatial_dataset
+from csep.plots import plot_gridded_dataset
 
 
 # idea: should this be a SpatialDataSet and the class below SpaceMagnitudeDataSet, bc of functions like
@@ -432,17 +432,27 @@ def load_ascii(cls, ascii_fname, start_date=None, end_date=None, name=None, swap
         gds = cls(start_date, end_date, magnitudes=mws, name=name, region=region, data=rates)
         return gds
 
-    def plot(self, ax=None, show=False, log=True, extent=None, set_global=False, plot_args=None):
-        """ Plot gridded forecast according to plate-carree projection
+    def plot(self, ax=None, show=False, log=True, extent=None, set_global=False, plot_args=None,
+             **kwargs):
+        """ Plot the spatial rate of the forecast
+
+        See :func:`csep.utils.plots.plot_gridded_dataset` for a detailed description of the
+        keyword arguments.
 
         Args:
-            show (bool): if true, show the figure. this call is blocking.
-            plot_args (optional/dict): dictionary containing plotting arguments for making figures
+            ax (`matplotlib.pyplot.axes`): Previous axes onto which catalog can be drawn
+            show (bool): If True, shows the figure.
+            log (bool): If True, plots the base-10 logarithm of the spatial rates
+            extent (list): Force an extent [lon_min, lon_max, lat_min, lat_max]
+            set_global (bool): Whether to plot using a global projection
+            **kwargs (dict): Keyword arguments passed to
+                :func:`csep.utils.plots.plot_gridded_dataset`
 
         Returns:
             axes: matplotlib.Axes.axes
         """
-        # no mutable function arguments
+
+
         if self.start_time is None or self.end_time is None:
             time = 'forecast period'
         else:
@@ -451,19 +461,24 @@ def plot(self, ax=None, show=False, log=True, extent=None, set_global=False, plo
             time = f'{round(end-start,3)} years'
 
         plot_args = plot_args or {}
-        plot_args.setdefault('figsize', (10, 10))
-        plot_args.setdefault('title', self.name)
-
+        plot_args.update({
+             '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
         if log:
             plot_args.setdefault('clabel', f'log10 M{self.min_magnitude}+ rate per cell per {time}')
             with numpy.errstate(divide='ignore'):
-                ax = plot_spatial_dataset(numpy.log10(self.spatial_counts(cartesian=True)), self.region, ax=ax,
-                                          show=show, extent=extent, set_global=set_global, plot_args=plot_args)
+                ax = plot_gridded_dataset(numpy.log10(self.spatial_counts(cartesian=True)), self.region, ax=ax,
+                                          show=show, extent=extent, set_global=set_global,
+                                          **plot_args)
         else:
             plot_args.setdefault('clabel', f'M{self.min_magnitude}+ rate per cell per {time}')
-            ax = plot_spatial_dataset(self.spatial_counts(cartesian=True), self.region, ax=ax,show=show, extent=extent,
-                                      set_global=set_global, plot_args=plot_args)
+            ax = plot_gridded_dataset(self.spatial_counts(cartesian=True), self.region, ax=ax, show=show, extent=extent,
+                                      set_global=set_global, **plot_args)
         return ax
 
 
@@ -654,7 +669,7 @@ def magnitude_counts(self):
             self.get_expected_rates()
         return self.expected_rates.magnitude_counts()
 
-    def get_event_counts(self, verbose=True):
+    def get_event_counts(self, verbose=False):
         """ Returns a numpy array containing the number of event counts for each catalog.
 
             Note: This function can take a while to compute if called without already iterating through a forecast that
@@ -715,7 +730,7 @@ def get_expected_rates(self, verbose=False):
                                                   magnitudes=self.magnitudes, name=self.name)
             return self.expected_rates
 
-    def plot(self, plot_args = None, verbose=True, **kwargs):
+    def plot(self, plot_args = None, verbose=False, **kwargs):
         plot_args = plot_args or {}
         if self.expected_rates is None:
             self.get_expected_rates(verbose=verbose)

csep/core/forecasts.py

@@ -47,7 +47,7 @@ def paired_t_test(forecast, benchmark_forecast, observed_catalog,
     result.name = 'Paired T-Test'
     result.test_distribution = (out['ig_lower'], out['ig_upper'])
     result.observed_statistic = out['information_gain']
-    result.quantile = (out['t_statistic'], out['t_critical'])
+    result.quantile = (out['t_statistic'], out['t_critical'], alpha)
     result.sim_name = (forecast.name, benchmark_forecast.name)
     result.obs_name = observed_catalog.name
     result.status = 'normal'
@@ -601,7 +601,7 @@ def _simulate_catalog(num_events, sampling_weights, sim_fore,
 
 def _poisson_likelihood_test(forecast_data, observed_data,
                              num_simulations=1000, random_numbers=None,
-                             seed=None, use_observed_counts=True, verbose=True,
+                             seed=None, use_observed_counts=True, verbose=False,
                              normalize_likelihood=False):
     """
 	Computes the likelihood-test from CSEP using an efficient simulation based approach.

csep/core/poisson_evaluations.py

@@ -658,7 +658,11 @@ def get_cartesian(self, data):
         """Returns 2d ndrray representation of the data set, corresponding to the bounding box.
 
         Args:
-            data:
+            data: An array of values, whose indices corresponds to each cell of the region
+
+        Returns:
+            A 2D array of values, corresponding to the original data projected onto the region.
+            Values outside the region polygon are represented as np.nan
         """
         assert len(data) == len(self.polygons)
         results = numpy.zeros(self.bbox_mask.shape[:2])