changes to catalogs and added defaults for plot functions. closes #99 and closes #95. (#100)

- added defaults for plot catalog to make the basic plot look prettier
- added defaults for forecast plots
- forecast.plot() properly returns the figure axes
- added function to compute tight_bbox around region
- added option region_border to plot tight bbox around regions in catalog plots
- added option region_border to plot tight boox around region in forecast plots
- added start of plots.rst for more in-depth documentation on plotting
- fixes the catalog evaluation example. closes #99

Author: wsavran

csep/core/catalogs.py

@@ -845,14 +845,33 @@ 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': [4.0, 5.0, 6.0, 7.0]
+        }
+        # 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
+        except AttributeError:
+            pass
 
         plot_args = plot_args or {}
-        plot_args.setdefault('figsize', (10, 10))
-        plot_args.setdefault('title', self.name)
+        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)
+                          set_global=set_global, plot_args=plot_args_default)
         return ax
 
 

csep/core/forecasts.py

@@ -673,10 +673,20 @@ def get_expected_rates(self, verbose=False, return_skipped=False):
             else:
                 return self.expected_rates
 
-    def plot(self, **kwargs):
+    def plot(self, plot_args = None, **kwargs):
+        plot_args = plot_args or {}
         if self.expected_rates is None:
             self.get_expected_rates()
-        self.expected_rates.plot(**kwargs)
+        args_dict = {'title': self.name,
+                     'grid_labels': True,
+                     'grid': True,
+                     'borders': True,
+                     'feature_lw': 0.5,
+                     'basemap': 'ESRI_terrain',
+                     }
+        args_dict.update(plot_args)
+        ax = self.expected_rates.plot(**kwargs, plot_args=args_dict)
+        return ax
 
     def get_dataframe(self):
         """Return a single dataframe with all of the events from all of the catalogs."""

csep/core/regions.py

@@ -11,7 +11,7 @@
 import pyproj
 
 # PyCSEP imports
-from csep.utils.calc import bin1d_vec, cleaner_range
+from csep.utils.calc import bin1d_vec, cleaner_range, first_nonnan, last_nonnan
 from csep.utils.scaling_relationships import WellsAndCoppersmith
 
 def california_relm_collection_region(dh_scale=1, magnitudes=None, name="relm-california-collection"):
@@ -730,3 +730,26 @@ def _build_bitmask_vec(self):
 
         return a, xs, ys
 
+    def tight_bbox(self):
+        # creates tight bounding box around the region, probably a faster way to do this.
+        ny, nx = self.idx_map.shape
+        asc = []
+        desc = []
+        for j in range(ny):
+            row = self.idx_map[j, :]
+            argmin = first_nonnan(row)
+            argmax = last_nonnan(row)
+            # points are stored clockwise
+            poly_min = self.polygons[int(row[argmin])].points
+            asc.insert(0, poly_min[0])
+            asc.insert(0, poly_min[1])
+            poly_max = self.polygons[int(row[argmax])].points
+            desc.append(poly_max[3])
+            desc.append(poly_max[2])
+        # close the loop
+        poly = np.array(asc + desc)
+        sorted_idx = np.sort(np.unique(poly, return_index=True, axis=0)[1], kind='stable')
+        unique_poly = poly[sorted_idx]
+        unique_poly = np.append(unique_poly, [unique_poly[0, :]], axis=0)
+        return unique_poly
+

csep/utils/calc.py

@@ -213,4 +213,13 @@ def cleaner_range(start, end, h):
     start = numpy.floor(const * start)
     end = numpy.floor(const * end)
     d = const * h
-    return numpy.arange(start, end + d / 2, d) / const
+    return numpy.arange(start, end + d / 2, d) / const
+
+def first_nonnan(arr, axis=0, invalid_val=-1):
+    mask = arr==arr
+    return numpy.where(mask.any(axis=axis), mask.argmax(axis=axis), invalid_val)
+
+def last_nonnan(arr, axis=0, invalid_val=-1):
+    mask = arr==arr
+    val = arr.shape[axis] - numpy.flip(mask, axis=axis).argmax(axis=axis) - 1
+    return numpy.where(mask.any(axis=axis), val, invalid_val)

csep/utils/plots.py

@@ -652,11 +652,6 @@ def plot_catalog(catalog, ax=None, show=False, extent=None, set_global=False, pl
 
     """
     # Get spatial information for plotting
-    bbox = catalog.get_bbox()
-    if extent is None:
-        dh = (bbox[1] - bbox[0])/20.
-        dv = (bbox[3] - bbox[2]) / 20.
-        extent = [bbox[0] - dh, bbox[1]+dh, bbox[2] -dv, bbox[3] + dv]
 
     # Retrieve plot arguments
     plot_args = plot_args or {}
@@ -674,6 +669,7 @@ def plot_catalog(catalog, ax=None, show=False, extent=None, set_global=False, pl
     legend_loc = plot_args.get('legend_loc', 1)
     mag_ticks = plot_args.get('mag_ticks', False)
     labelspacing = plot_args.get('labelspacing', 1)
+    region_border = plot_args.get('region_border', True)
     # cartopy properties
     projection = plot_args.get('projection', ccrs.PlateCarree(central_longitude=0.0))
     grid = plot_args.get('grid', True)
@@ -685,8 +681,19 @@ def plot_catalog(catalog, ax=None, show=False, extent=None, set_global=False, pl
     linecolor = plot_args.get('linecolor', 'black')
 
 
-    # Instantiage GeoAxes object
+    bbox = catalog.get_bbox()
+    if plot_args['region_border']:
+        try:
+            bbox = catalog.region.get_bbox()
+        except AttributeError:
+            pass
 
+    if extent is None:
+        dh = (bbox[1] - bbox[0])/20.
+        dv = (bbox[3] - bbox[2]) / 20.
+        extent = [bbox[0] - dh, bbox[1]+dh, bbox[2] -dv, bbox[3] + dv]
+
+    # Instantiage GeoAxes object
     if ax is None:
         fig = pyplot.figure(figsize=figsize)
         ax = fig.add_subplot(111, projection=projection)
@@ -727,6 +734,13 @@ def size_map(markersize, values, scale):
                   loc=legend_loc, title=r"Magnitudes",title_fontsize=16,
                   labelspacing=labelspacing, handletextpad=5, framealpha=False)
 
+    if region_border:
+        try:
+            pts = catalog.region.tight_bbox()
+            ax.plot(pts[:, 0], pts[:, 1], lw=1, color='black')
+        except AttributeError:
+            print("unable to get tight bbox")
+
     # Gridline options
     if grid:
         gl = ax.gridlines(draw_labels=grid_labels, alpha=0.5)
@@ -801,6 +815,7 @@ def plot_spatial_dataset(gridded, region, ax=None, show=False, extent=None, set_
     borders = plot_args.get('borders', False)
     linewidth = plot_args.get('linewidth', True)
     linecolor = plot_args.get('linecolor', 'black')
+    region_border = plot_args.get('region_border', True)
     # color bar properties
     cmap = plot_args.get('cmap', None)
     clabel = plot_args.get('clabel', '')
@@ -857,6 +872,10 @@ def plot_spatial_dataset(gridded, region, ax=None, show=False, extent=None, set_
         gl.xformatter = LONGITUDE_FORMATTER
         gl.yformatter = LATITUDE_FORMATTER
 
+    if region_border:
+        pts = region.tight_bbox()
+        ax.plot(pts[:,0], pts[:,1], lw=1, color='black')
+
     # matplotlib figure options
     ax.set_title(title, y=1.06)
     if filename is not None:

docs/concepts/plots.rst

@@ -0,0 +1,27 @@
+.. _plots-reference:
+
+#####
+Plots
+#####
+
+PyCSEP provides several functions to produce commonly used plots, such as an earthquake forecast or the evaluation catalog
+or perhaps a combination of the two.
+
+.. contents:: Table of Contents
+    :local:
+    :depth: 2
+
+************
+Introduction
+************
+
+
+
+**************
+Plot arguments
+**************
+
+***************
+Available plots
+***************
+

examples/tutorials/catalog_forecast_evaluation.py

@@ -33,8 +33,8 @@
 # Forecasts should define a time horizon in which they are valid. The choice is flexible for catalog-based forecasts, because
 # the catalogs can be filtered to accommodate multiple end-times. Conceptually, these should be separate forecasts.
 
-start_time = time_utils.strptime_to_utc_datetime("1992-06-28 11:57:34.14")
-end_time = time_utils.strptime_to_utc_datetime("1992-07-28 11:57:34.14")
+start_time = time_utils.strptime_to_utc_datetime("1992-06-28 11:57:35.0")
+end_time = time_utils.strptime_to_utc_datetime("1992-07-28 11:57:35.0")
 
 ####################################################################################################################################
 # Define spatial and magnitude regions
@@ -72,7 +72,8 @@
 
 forecast = csep.load_catalog_forecast(datasets.ucerf3_ascii_format_landers_fname,
                                       start_time = start_time, end_time = end_time,
-                                      region = space_magnitude_region)
+                                      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}']
@@ -92,7 +93,12 @@
 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)
+
+# Plot the catalog
+comcat_catalog.plot()
 
 ####################################################################################################################################
 # Perform number test