v0.6.1 release branch

csep/core/regions.py

@@ -517,6 +517,9 @@ def __init__(self, polygons, dh, name='cartesian2d', mask=None):
         # index values of polygons array into the 2d cartesian grid, based on the midpoint.
         self.xs = xs
         self.ys = ys
+        # Bounds [origin, top_right]
+        orgs = self.origins()
+        self.bounds = numpy.column_stack((orgs, orgs + dh))
 
     def __eq__(self, other):
         return self.to_dict() == other.to_dict()
@@ -772,13 +775,16 @@ def geographical_area_from_bounds(lon1, lat1, lon2, lat2):
     Returns:
         Area of cell in Km2
     """
-    earth_radius_km = 6371.
-    R2 = earth_radius_km ** 2
-    rad_per_deg = numpy.pi / 180.0e0
+    if lon1 == lon2 or lat1 == lat2:
+        return 0
+    else:
+        earth_radius_km = 6371.
+        R2 = earth_radius_km ** 2
+        rad_per_deg = numpy.pi / 180.0e0
 
-    strip_area_steradian = 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat1) * rad_per_deg)) \
+        strip_area_steradian = 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat1) * rad_per_deg)) \
                            - 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat2) * rad_per_deg))
-    area_km2 = strip_area_steradian * R2 / (360.0 / (lon2 - lon1))
+        area_km2 = strip_area_steradian * R2 / (360.0 / (lon2 - lon1))
     return area_km2
 
 def quadtree_grid_bounds(quadk):

csep/utils/plots.py

@@ -6,6 +6,7 @@
 
 import scipy.stats
 import matplotlib
+import matplotlib.lines
 from matplotlib import cm
 from matplotlib.collections import PatchCollection
 import matplotlib.pyplot as pyplot
@@ -1512,6 +1513,7 @@ def plot_comparison_test(results_t, results_w=None, axes=None, plot_args=None):
 
     return ax
 
+
 def plot_poisson_consistency_test(eval_results, normalize=False, one_sided_lower=False, axes=None, plot_args=None, show=False):
     """ Plots results from CSEP1 tests following the CSEP1 convention.
 
@@ -1883,3 +1885,249 @@ def add_labels_for_publication(figure, style='bssa', labelsize=16):
             ax.annotate(f'({annot})', (0.025, 1.025), xycoords='axes fraction', fontsize=labelsize)
 
     return
+
+
+def plot_consistency_test(eval_results, normalize=False, one_sided_lower=True, plot_args=None, variance=None):
+    """ Plots results from CSEP1 tests following the CSEP1 convention.
+
+    Note: All of the evaluations should be from the same type of evaluation, otherwise the results will not be
+          comparable on the same figure.
+
+    Args:
+        eval_results (list): Contains the tests results :class:`csep.core.evaluations.EvaluationResult` (see note above)
+        normalize (bool): select this if the forecast likelihood should be normalized by the observed likelihood. useful
+                          for plotting simulation based simulation tests.
+        one_sided_lower (bool): select this if the plot should be for a one sided test
+        plot_args(dict): optional argument containing a dictionary of plotting arguments, with keys as strings and items as described below
+
+    Optional plotting arguments:
+        * figsize: (:class:`list`/:class:`tuple`) - default: [6.4, 4.8]
+        * title: (:class:`str`) - default: name of the first evaluation result type
+        * title_fontsize: (:class:`float`) Fontsize of the plot title - default: 10
+        * xlabel: (:class:`str`) - default: 'X'
+        * xlabel_fontsize: (:class:`float`) - default: 10
+        * xticks_fontsize: (:class:`float`) - default: 10
+        * ylabel_fontsize: (:class:`float`) - default: 10
+        * color: (:class:`float`/:class:`None`) If None, sets it to red/green according to :func:`_get_marker_style` - default: 'black'
+        * linewidth: (:class:`float`) - default: 1.5
+        * capsize: (:class:`float`) - default: 4
+        * hbars:  (:class:`bool`)  Flag to draw horizontal bars for each model - default: True
+        * tight_layout: (:class:`bool`) Set matplotlib.figure.tight_layout to remove excess blank space in the plot - default: True
+
+    Returns:
+        ax (:class:`matplotlib.pyplot.axes` object)
+    """
+
+
+    try:
+        results = list(eval_results)
+    except TypeError:
+        results = [eval_results]
+    results.reverse()
+    # Parse plot arguments. More can be added here
+    if plot_args is None:
+        plot_args = {}
+    figsize= plot_args.get('figsize', (7,8))
+    xlabel = plot_args.get('xlabel', 'X')
+    xlabel_fontsize = plot_args.get('xlabel_fontsize', None)
+    xticks_fontsize = plot_args.get('xticks_fontsize', None)
+    ylabel_fontsize = plot_args.get('ylabel_fontsize', None)
+    color = plot_args.get('color', 'black')
+    linewidth = plot_args.get('linewidth', None)
+    capsize = plot_args.get('capsize', 4)
+    hbars = plot_args.get('hbars', True)
+    tight_layout = plot_args.get('tight_layout', True)
+    percentile = plot_args.get('percentile', 95)
+
+    fig, ax = pyplot.subplots(figsize=figsize)
+    xlims = []
+
+    for index, res in enumerate(results):
+        # handle analytical distributions first, they are all in the form ['name', parameters].
+        if res.test_distribution[0] == 'poisson':
+            plow = scipy.stats.poisson.ppf((1 - percentile/100.)/2., res.test_distribution[1])
+            phigh = scipy.stats.poisson.ppf(1 - (1 - percentile/100.)/2., res.test_distribution[1])
+            observed_statistic = res.observed_statistic
+
+        elif res.test_distribution[0] == 'negative_binomial':
+            var = variance
+            observed_statistic = res.observed_statistic
+            mean = res.test_distribution[1]
+            upsilon = 1.0 - ((var - mean) / var)
+            tau = (mean**2 /(var - mean))
+            phigh = scipy.stats.nbinom.ppf((1 - percentile/100.)/2., tau, upsilon)
+            plow = scipy.stats.nbinom.ppf(1 - (1 - percentile/100.)/2., tau, upsilon)
+
+        # empirical distributions
+        else:
+            if normalize:
+                test_distribution = numpy.array(res.test_distribution) - res.observed_statistic
+                observed_statistic = 0
+            else:
+                test_distribution = numpy.array(res.test_distribution)
+                observed_statistic = res.observed_statistic
+            # compute distribution depending on type of test
+            if one_sided_lower:
+                plow = numpy.percentile(test_distribution, 5)
+                phigh = numpy.percentile(test_distribution, 100)
+            else:
+                plow = numpy.percentile(test_distribution, 2.5)
+                phigh = numpy.percentile(test_distribution, 97.5)
+
+        if not numpy.isinf(observed_statistic): # Check if test result does not diverges
+            low = observed_statistic - plow
+            high = phigh - observed_statistic
+            ax.errorbar(observed_statistic, index, xerr=numpy.array([[low, high]]).T,
+                        fmt=_get_marker_style(observed_statistic, (plow, phigh), one_sided_lower),
+                        capsize=4, linewidth=linewidth, ecolor=color, markersize = 10, zorder=1)
+            # determine the limits to use
+            xlims.append((plow, phigh, observed_statistic))
+            # we want to only extent the distribution where it falls outside of it in the acceptable tail
+            if one_sided_lower:
+                if observed_statistic >= plow and phigh < observed_statistic:
+                    # draw dashed line to infinity
+                    xt = numpy.linspace(phigh, 99999, 100)
+                    yt = numpy.ones(100) * index
+                    ax.plot(xt, yt, linestyle='--', linewidth=linewidth, color=color)
+
+        else:
+            print('Observed statistic diverges for forecast %s, index %i.'
+                  ' Check for zero-valued bins within the forecast'% (res.sim_name, index))
+            ax.barh(index, 99999, left=-10000, height=1, color=['red'], alpha=0.5)
+
+
+    try:
+        ax.set_xlim(*_get_axis_limits(xlims))
+    except ValueError:
+        raise ValueError('All EvaluationResults have infinite observed_statistics')
+    ax.set_yticks(numpy.arange(len(results)))
+    ax.set_yticklabels([res.sim_name for res in results], fontsize=14)
+    ax.set_ylim([-0.5, len(results)-0.5])
+    if hbars:
+        yTickPos = ax.get_yticks()
+        if len(yTickPos) >= 2:
+            ax.barh(yTickPos, numpy.array([99999] * len(yTickPos)), left=-10000,
+                    height=(yTickPos[1] - yTickPos[0]), color=['w', 'gray'], alpha=0.2, zorder=0)
+    ax.set_xlabel(xlabel, fontsize=14)
+    ax.tick_params(axis='x', labelsize=13)
+    if tight_layout:
+        ax.figure.tight_layout()
+        fig.tight_layout()
+    return ax
+
+
+def plot_pvalues_and_intervals(test_results, ax, var=None):
+    """ Plots p-values and intervals for a list of Poisson or NBD test results
+
+    Args:
+        test_results (list): list of EvaluationResults for N-test. All tests should use the same distribution
+                             (ie Poisson or NBD).
+        ax (matplotlib.axes.Axes.axis): axes to use for plot. create using matplotlib
+        var (float): variance of the NBD distribution. Must be used for NBD plots.
+
+    Returns:
+        ax (matplotlib.axes.Axes.axis): axes handle containing this plot
+
+    Raises:
+        ValueError: throws error if NBD tests are supplied without a variance
+    """
+
+    variance = var
+    percentile = 97.5
+    p_values = []
+
+    # Differentiate between N-tests and other consistency tests
+    if test_results[0].name == 'NBD N-Test' or test_results[0].name == 'Poisson N-Test':
+        legend_elements = [matplotlib.lines.Line2D([0], [0], marker='o', color='red', lw=0, label=r'p < 10e-5', markersize=10, markeredgecolor='k'),
+                           matplotlib.lines.Line2D([0], [0], marker='o', color='#FF7F50', lw=0, label=r'10e-5 $\leq$ p < 10e-4', markersize=10, markeredgecolor='k'),
+                           matplotlib.lines.Line2D([0], [0], marker='o', color='gold', lw=0, label=r'10e-4 $\leq$ p < 10e-3', markersize=10, markeredgecolor='k'),
+                           matplotlib.lines.Line2D([0], [0], marker='o', color='white', lw=0, label=r'10e-3 $\leq$ p < 0.0125', markersize=10, markeredgecolor='k'),
+                           matplotlib.lines.Line2D([0], [0], marker='o', color='skyblue', lw=0, label=r'0.0125 $\leq$ p < 0.025', markersize=10, markeredgecolor='k'),
+                           matplotlib.lines.Line2D([0], [0], marker='o', color='blue', lw=0, label=r'p $\geq$ 0.025', markersize=10, markeredgecolor='k')]
+        ax.legend(handles=legend_elements, loc=4, fontsize=13, edgecolor='k')
+        # Act on Negative binomial tests
+        if test_results[0].name == 'NBD N-Test':
+            if var is None:
+                raise ValueError("var must not be None if N-tests use the NBD distribution.")
+
+            for i in range(len(test_results)):
+                mean = test_results[i].test_distribution[1]
+                upsilon = 1.0 - ((variance - mean) / variance)
+                tau = (mean**2 /(variance - mean))
+                phigh97 = scipy.stats.nbinom.ppf((1 - percentile/100.)/2., tau, upsilon)
+                plow97 = scipy.stats.nbinom.ppf(1 - (1 - percentile/100.)/2., tau, upsilon)
+                low97 = test_results[i].observed_statistic - plow97
+                high97 = phigh97 - test_results[i].observed_statistic
+                ax.errorbar(test_results[i].observed_statistic, (len(test_results)-1) - i, xerr=numpy.array([[low97, high97]]).T, capsize=4, 
+                            color='slategray', alpha=1.0, zorder=0)
+                p_values.append(test_results[i].quantile[1] * 2.0) # Calculated p-values according to Meletti et al., (2021)
+
+                if p_values[i] < 10e-5:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='red', markersize=8, zorder=2)
+                if p_values[i] >= 10e-5 and p_values[i] < 10e-4:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='#FF7F50', markersize=8, zorder=2)
+                if p_values[i] >= 10e-4 and p_values[i] < 10e-3:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='gold', markersize=8, zorder=2)
+                if p_values[i] >= 10e-3 and p_values[i] < 0.0125:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='white', markersize=8, zorder=2)
+                if p_values[i] >= 0.0125 and p_values[i] < 0.025:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='skyblue', markersize=8, zorder=2)
+                if p_values[i] >= 0.025:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='blue', markersize=8, zorder=2)
+        # Act on Poisson N-test
+        if test_results[0].name == 'Poisson N-Test':
+            for i in range(len(test_results)):
+                plow97 = scipy.stats.poisson.ppf((1 - percentile/100.)/2., test_results[i].test_distribution[1])
+                phigh97 = scipy.stats.poisson.ppf(1 - (1 - percentile/100.)/2., test_results[i].test_distribution[1])
+                low97 = test_results[i].observed_statistic - plow97
+                high97 = phigh97 - test_results[i].observed_statistic
+                ax.errorbar(test_results[i].observed_statistic, (len(test_results)-1) - i, xerr=numpy.array([[low97, high97]]).T, capsize=4, 
+                            color='slategray', alpha=1.0, zorder=0)
+                p_values.append(test_results[i].quantile[1] * 2.0)
+                if p_values[i] < 10e-5:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='red', markersize=8, zorder=2)
+                elif p_values[i] >= 10e-5 and p_values[i] < 10e-4:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='#FF7F50', markersize=8, zorder=2)
+                elif p_values[i] >= 10e-4 and p_values[i] < 10e-3:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='gold', markersize=8, zorder=2)
+                elif p_values[i] >= 10e-3 and p_values[i] < 0.0125:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='white', markersize=8, zorder=2)
+                elif p_values[i] >= 0.0125 and p_values[i] < 0.025:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='skyblue', markersize=8, zorder=2)
+                elif p_values[i] >= 0.025:
+                    ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='blue', markersize=8, zorder=2)
+    # Operate on all other consistency tests
+    else:
+        for i in range(len(test_results)):
+            plow97 = numpy.percentile(test_results[i].test_distribution, 2.5)
+            phigh97 = numpy.percentile(test_results[i].test_distribution, 97.5)
+            low97 = test_results[i].observed_statistic - plow97
+            high97 = phigh97 - test_results[i].observed_statistic  
+            ax.errorbar(test_results[i].observed_statistic, (len(test_results)-1) -i, xerr=numpy.array([[low97, high97]]).T, capsize=4, 
+                        color='slategray', alpha=1.0, zorder=0)
+            p_values.append(test_results[i].quantile)
+
+            if p_values[i] < 10e-5:
+                ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='red', markersize=8, zorder=2)
+            elif p_values[i] >= 10e-5 and p_values[i] < 10e-4:
+                ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='#FF7F50', markersize=8, zorder=2)
+            elif p_values[i] >= 10e-4 and p_values[i] < 10e-3:
+                ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='gold', markersize=8, zorder=2)
+            elif p_values[i] >= 10e-3  and p_values[i] < 0.025:
+                ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='white', markersize=8, zorder=2)
+            elif p_values[i] >= 0.025 and p_values[i] < 0.05:
+                ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='skyblue', markersize=8, zorder=2)
+            elif p_values[i] >= 0.05:
+                ax.plot(test_results[i].observed_statistic, (len(test_results)-1) - i, marker='o', color='blue', markersize=8, zorder=2)
+
+        legend_elements = [
+            matplotlib.lines.Line2D([0], [0], marker='o', color='red', lw=0, label=r'p < 10e-5', markersize=10, markeredgecolor='k'),
+            matplotlib.lines.Line2D([0], [0], marker='o', color='#FF7F50', lw=0, label=r'10e-5 $\leq$ p < 10e-4', markersize=10, markeredgecolor='k'),
+            matplotlib.lines.Line2D([0], [0], marker='o', color='gold', lw=0, label=r'10e-4 $\leq$ p < 10e-3', markersize=10, markeredgecolor='k'),
+            matplotlib.lines.Line2D([0], [0], marker='o', color='white', lw=0, label=r'10e-3 $\leq$ p < 0.025', markersize=10, markeredgecolor='k'),
+            matplotlib.lines.Line2D([0], [0], marker='o', color='skyblue', lw=0, label=r'0.025 $\leq$ p < 0.05', markersize=10, markeredgecolor='k'),
+            matplotlib.lines.Line2D([0], [0], marker='o', color='blue', lw=0, label=r'p $\geq$ 0.05', markersize=10, markeredgecolor='k')]
+
+        ax.legend(handles=legend_elements, loc=4, fontsize=13, edgecolor='k') 
+
+    return ax        

csep/utils/readers.py

@@ -719,7 +719,7 @@ def _parse_datetime_to_zmap(date, time):
         out['second'] = dt.second
         return out
 
-def load_quadtree_forecast(ascii_fname):
+def quadtree_ascii_loader(ascii_fname):
     """ Load quadtree forecasted stored as ascii text file
 
         Note: This function is adapted form csep.forecasts.load_ascii
@@ -756,4 +756,33 @@ def load_quadtree_forecast(ascii_fname):
     # reshape rates into correct 2d format
     rates = data[:, -1].reshape(n_poly, n_mag_bins)
 
+    return rates, region, mws
+
+
+def quadtree_csv_loader(csv_fname):
+    """ Load quadtree forecasted stored as csv file
+
+        The format expects forecast as a comma separated file, in which first column corresponds to quadtree grid cell (quadkey).
+        The second and thrid columns indicate depth range.
+        The corresponding enteries in the respective row are forecast rates corresponding to the magnitude bins.
+        The first line of forecast is a header, and its format is listed here:
+            'Quadkey', depth_min, depth_max, Mag_0, Mag_1, Mag_2, Mag_3 , ....
+             Quadkey is a string. Rest of the values are floats.
+        For the purposes of defining region objects quadkey is used.
+
+        We assume that the starting value of magnitude bins are provided in the header.
+        Args:
+            csv_fname: file name of csep forecast in csv format
+        Returns:
+            rates, region, mws (numpy.ndarray, QuadtreeRegion2D, numpy.ndarray): rates, region, and magnitude bins needed
+                                                                                 to define QuadTree forecasts
+     """
+
+    data = numpy.genfromtxt(csv_fname, dtype='str', delimiter=',')
+    quadkeys = data[1:, 0]
+    mws = data[0, 3:]
+    rates = data[1:, 3:]
+    rates = rates.astype(float)
+    region = QuadtreeGrid2D.from_quadkeys(quadkeys, magnitudes=mws)
+
     return rates, region, mws

requirements.txt

@@ -1,4 +1,4 @@
-numpy<=1.21.5
+numpy
 scipy
 pandas
 matplotlib

requirements.yml

@@ -4,7 +4,7 @@ channels:
   - defaults
 dependencies:
   - python>=3.7
-  - numpy<=1.21.5
+  - numpy
   - pandas
   - scipy
   - matplotlib

setup.py

@@ -26,7 +26,7 @@ def get_version():
     long_description=long_description,
     long_description_content_type='text/markdown',
     install_requires = [
-        'numpy',
+        'numpy<=1.21.5',
         'scipy',
         'pandas',
         'matplotlib',