Plot Refactoring

.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

@@ -333,7 +333,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

csep/core/binomial_evaluations.py

@@ -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)

csep/core/catalog_evaluations.py

@@ -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,
@@ -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',
@@ -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/catalogs.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.utils.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/forecasts.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.utils.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
 
 

csep/core/poisson_evaluations.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'

csep/core/regions.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])

csep/models.py

@@ -143,7 +143,7 @@ def plot(self, show=False, plot_args=None):
                               'bins': bins}
         # looks funny, but will update the defaults with the user defined arguments
         plot_args_defaults.update(plot_args)
-        ax = plots.plot_number_test(self, show=show, plot_args=plot_args)
+        ax = plots.plot_test_distribution(self, show=show)
         return ax
 
 
@@ -160,7 +160,7 @@ def plot(self, show=False, plot_args=None):
                               'bins': 'auto'}
         # looks funny, but will update the defaults with the user defined arguments
         plot_args_defaults.update(plot_args)
-        ax = plots.plot_likelihood_test(self, show=show, plot_args=plot_args)
+        ax = plots.plot_test_distribution(self, show=show)
         return ax
 
 
@@ -175,7 +175,7 @@ def plot(self, show=False, plot_args=None):
                               'title': 'Magnitude Test',
                               'bins': 'auto'}
         plot_args_defaults.update(plot_args)
-        ax = plots.plot_magnitude_test(self, show=show, plot_args=plot_args)
+        ax = plots.plot_test_distribution(self, show=show)
         return ax
 
 
@@ -194,7 +194,7 @@ def plot(self, show=False, plot_args=None):
         }
         # looks funny, but will update the defaults with the user defined arguments
         plot_args_defaults.update(plot_args)
-        ax = plots.plot_spatial_test(self, show=show, plot_args=plot_args)
+        ax = plots.plot_test_distribution(self, show=show)
         return ax
 
 
@@ -211,7 +211,7 @@ def plot(self, show=False, axes=None, plot_args=None):
             'title': self.name
         }
         plot_args_defaults.update(plot_args)
-        ax = plots.plot_calibration_test(self, show=show, axes=axes, plot_args=plot_args)
+        ax = plots.plot_calibration_test(self, show=show, ax=axes, plot_args=plot_args)
         return ax
 
 
@@ -361,4 +361,4 @@ def from_great_circle_radius(cls, centroid, radius, num_points=10):
         # get new lons and lats
         endlon, endlat, backaz = geod.fwd(center_lons, center_lats, azim, radius)
         # class method
-        return cls(np.column_stack([endlon, endlat]))
+        return cls(np.column_stack([endlon, endlat]))

csep/utils/datasets.py

@@ -5,6 +5,7 @@
 _gridded_forecast_root = os.path.join(_root_dir, 'artifacts', 'ExampleForecasts', 'GriddedForecasts')
 _catalog_forecast_root = os.path.join(_root_dir, 'artifacts', 'ExampleForecasts', 'CatalogForecasts')
 _observed_catalog_root = os.path.join(_root_dir, 'artifacts', 'ObservedCatalogs')
+_basemap_root = os.path.join(_root_dir, 'artifacts', 'ExampleBasemaps')
 _polygon_region_root = os.path.join(_root_dir, 'artifacts', 'Regions', 'Polygons')
 _l_test_example_root = os.path.join(_root_dir, 'artifacts', 'ExampleResults', 'Poisson_l-test')
 
@@ -35,6 +36,10 @@
     _polygon_region_root, 'Italy', 'ItalyCollectionPolygon.txt'
 )
 
+# basemap file name
+basemap_california = os.path.join(_basemap_root, 'basemap_california.tif')
+
+
 # Example test results
 l_test_examples = [os.path.join(_l_test_example_root, j) for j in os.listdir(_l_test_example_root)]