Plot Refactoring

.github/workflows/build-test.yml

@@ -34,10 +34,12 @@ jobs:
       - name: Install pyCSEP
         run: |
           pip install --no-deps -e .
+          pip install rasterio
           python -c "import csep; print('Version: ', csep.__version__)"
 
       - name: Test with pytest
         run: |
+
           pip install vcrpy==4.3.1 pytest pytest-cov
           pytest --cov=./ --cov-config=.coveragerc
 

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

@@ -55,6 +55,7 @@ def number_test(forecast, observed_catalog, verbose=True):
                                      obs_name=observed_catalog.name)
     return result
 
+
 def spatial_test(forecast, observed_catalog, verbose=True):
     """ Performs spatial test for catalog-based forecasts.
 
@@ -129,7 +130,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,
@@ -140,6 +141,7 @@ def spatial_test(forecast, observed_catalog, verbose=True):
 
     return result
 
+
 def magnitude_test(forecast, observed_catalog, verbose=True):
     """ Performs magnitude test for catalog-based forecasts """
     test_distribution = []
@@ -152,7 +154,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',
@@ -215,6 +217,7 @@ def magnitude_test(forecast, observed_catalog, verbose=True):
 
     return result
 
+
 def pseudolikelihood_test(forecast, observed_catalog, verbose=True):
     """ Performs the spatial pseudolikelihood test for catalog forecasts.
 
@@ -310,6 +313,7 @@ def pseudolikelihood_test(forecast, observed_catalog, verbose=True):
 
     return result
 
+
 def calibration_test(evaluation_results, delta_1=False):
     """ Perform the calibration test by computing a Kilmogorov-Smirnov test of the observed quantiles against a uniform
     distribution.

csep/core/catalogs.py

@@ -829,7 +829,8 @@ 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):
+    def plot(self, ax=None, show=False, extent=None, set_global=False, plot_args=None,
+             **kwargs):
         """ Plot catalog according to plate-carree projection
 
         Args:
@@ -844,17 +845,7 @@ def plot(self, ax=None, show=False, extent=None, set_global=False, plot_args=Non
 
         # 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
+             'basemap': kwargs.pop('basemap', None) or 'ESRI_terrain',
         }
 
         # Plot the region border (if it exists) by default
@@ -869,8 +860,8 @@ def plot(self, ax=None, show=False, extent=None, set_global=False, plot_args=Non
         plot_args_default.update(plot_args)
 
         # 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_default, **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
@@ -427,7 +427,8 @@ 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):
+    def plot(self, ax=None, show=False, log=True, extent=None, set_global=False, plot_args=None,
+             **kwargs):
         """ Plot gridded forecast according to plate-carree projection
 
         Args:
@@ -446,19 +447,23 @@ 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', None) or 'ESRI_terrain',
+             'title': kwargs.pop('title', None) or self.name,
+             'figsize': kwargs.pop('figsize', None) or (8, 8),
+        })
+        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

@@ -674,7 +674,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_distribution_test(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_distribution_test(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_distribution_test(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_distribution_test(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)]