Feature 47 diff index

metplotpy/difficulty_index/README_plot_difficulty_index.rst

@@ -0,0 +1,58 @@
+Difficulty Index
+================
+Written by Bill Campbell and Liz Satterfield (NRL)
+Modified by Lindsay Blank (NCAR)
+Date Modified: October 21, 2020
+
+Background
+==========
+
+The overall aim of this work is to graphically represent the expected difficulty of a decision based on a set of forecasts (ensemble) of, e.g., significant wave height as a function of space and time. There are two basic factors that can make a decision difficult. The first factor is the proximity of the ensemble mean forecast to a decision threshold, e.g. 12 ft seas. If the ensemble mean is either much lower or much higher than the threshold, the decision is easier; if it is closer to the threshold, the decision is harder. The second factor is the forecast precision, or ensemble spread. The greater the spread around the ensemble mean, the more likely it is that there will be ensemble members both above and below the decision threshold, making the decision harder. (A third factor that we will not address here is undiagnosed systematic error, which adds uncertainty in a similar way to ensemble spread.) The challenge is combing these factors into a continuous function that allows the user to assess relative risk.
+
+
+Pre-requisites:
+===============
+Python packages:
+---------------
+- numpy
+- matplotlib version 2.2.2
+- scipy
+
+Python 3.6
+
+Input files:
+------------
+**Need to find out from NRL**
+
+Input Required:
+===============
+The function to plot the difficulty index is as follows: 
+
+def plot_field(field, lats, lons, vmin, vmax,
+        xlab, ylab, cmap, clab, title):
+
+The following is required to plot the difficulty index:
+#. **field:** Difficulty index values (see METcalcpy resources for how to calculate the difficulty index)
+#. **lats:** Latitude values
+#. **lons:** Longitude values
+#. **vmin:** Minimum value on the colorbar
+#. **vmax:** Maximum value on the colorbar
+#. **xlab:** x-axis label
+#. **ylab:** y-axis label
+#. **cmap:** Color map for plot
+#. **clab:** Label for colorbar
+#. **title:** Plot title
+
+
+Output generated:
+=================
+Calling the plot_field function in plot_difficulty_index.py will return a difficulty index plot.
+
+How to run:
+==========
+Run a 'pip install -e .' in $METplotpy to add metplotpy to the path. To call plot_difficulty_index, add the following import statement to the script:
+
+'from metplotpy.difficulty_index.plot_difficulty_index import plot_field'
+
+To see an example, please look at the following:
+$METplotpy/test/difficulty_index/example_difficulty_index.py

metplotpy/difficulty_index/mycolormaps.py

@@ -0,0 +1,77 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr  8 16:57:07 2020
+
+@author: campbell
+"""
+import numpy as np
+import matplotlib.colors as colors
+from scipy.interpolate import pchip
+
+
+def spectral(n=100):
+    """
+    Visually linear colormap based on linspecer Matlab map.
+    A vast improvement on the jet colormap
+    """
+    if n == 1:
+        cmap = np.array([0.2005, 0.5593, 0.7380])
+    elif n == 2:
+        cmap = np.array([[0.2005, 0.5593, 0.7380],
+                         [0.9684, 0.4799, 0.2723]])
+    else:
+        cmapp = np.array([
+                [158, 1, 66],     # dark red
+                [213, 62, 79],    # red
+                [244, 109, 67],   # dark orange
+                [253, 174, 97],   # orange
+                [254, 224, 139],  # light orange
+                [255, 255, 191],  # light yellow
+                [230, 245, 152],  # yellow green
+                [171, 221, 164],  # light green
+                [102, 194, 165],  # green
+                [50, 136, 189],   # blue
+                [94, 79, 162],    # violet
+                ])
+        x = np.linspace(0, n-1, np.shape(cmapp)[0])
+        xi = np.arange(n)
+        interp = pchip(x, cmapp, axis=0)
+        cmap = np.flipud(interp(xi)) / 255
+        cmap = colors.ListedColormap(cmap, name='spectral')
+
+    return cmap
+
+
+def stoplight(n=100):
+    """
+    An attempta at a visually nice stoplight colormap
+    """
+    if n == 1:
+        cmap = np.array([0.2005, 0.5593, 0.7380])
+    elif n == 2:
+        cmap = np.array([[0.2005, 0.5593, 0.7380],
+                         [0.9684, 0.4799, 0.2723]])
+    else:
+        cmapp = np.array([
+                # [90, 1, 33],     # very dark red
+                [158, 1, 66],     # dark red
+                [213, 62, 79],    # red
+                [244, 109, 67],   # dark orange
+                [253, 174, 97],   # orange
+                [254, 224, 139],  # light orange
+                [255, 255, 191],  # light yellow
+                [230, 245, 152],  # yellow green
+                [171, 221, 164],  # light green
+                # [102, 194, 165],  # green
+                [80, 165, 110],   # dark green
+                [40, 115, 50],   # very dark green
+                # [50, 136, 189],   # blue
+                # [94, 79, 162],    # violet
+                ])
+        x = np.linspace(0, n-1, np.shape(cmapp)[0])
+        xi = np.arange(n)
+        interp = pchip(x, cmapp, axis=0)
+        cmap = np.flipud(interp(xi)) / 255
+        cmap = colors.ListedColormap(cmap, name='spectral')
+
+    return cmap

metplotpy/difficulty_index/plot_difficulty_index.py

@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+"""
+Program name: plot_difficulty_index.py
+
+Plot forecast decision difficulty indices.
+Plot a set of decision difficulty indices for forecasts of postive
+definite quantities such as wind speed and wave height.
+
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from metcalcpy.calc_difficulty_index import forecast_difficulty as di
+from metcalcpy.piecewise_linear import PiecewiseLinear as plin
+import metplotpy.difficulty_index.mycolormaps as mcmap
+
+__author__ = 'Bill Campbell (NRL) and Lindsay Blank (NCAR)'
+__version__ = '0.1.0'
+__email__ = '[email protected]'
+
+# Enforce positive definiteness of quantities such as standard deviations
+EPS = np.finfo(np.float32).eps
+# Only allow 2D fields for now
+FIELD_DIM = 2
+
+def plot_field(field, lats, lons, vmin, vmax,
+        xlab, ylab, cmap, clab, title):
+    """
+    Parameters
+    ----------
+    field : 2D numpy array
+        Field you want to plot (difficulty index)
+    lats : 1D numpy array
+        Latitude values
+    lons : 1D numpy array
+        Longitude values
+    vmin : float
+        Minimum value on the colorbar
+    vmax : float
+        Maximum value on the colorbar
+    xlab : String
+        x-axis label
+    ylab : String
+        y-axis label
+    cmap: Matplotlib Colormap Class Object
+        Color map for plot
+    clab: String
+        Label for colorbar
+    title: String
+        Plot title
+
+    Returns
+    -------
+    fig : plot
+        Difficulty index plot
+    """
+
+    X, Y = np.meshgrid(lons, lats, indexing='ij')
+    fig, ax = plt.subplots(figsize=(8, 5))
+    plt.title(title)
+    if cmap is None:
+        cmap = mcmap.stoplight()
+
+    plt.pcolormesh(X, Y, field.T, shading='interp', cmap=cmap)
+    cbar = plt.colorbar(orientation='horizontal', aspect=30)
+    cbar.set_label(clab)
+    plt.clim(vmin=vmin, vmax=vmax)
+    plt.xlabel(xlab)
+    plt.ylabel(ylab)
+
+    return fig
+
+if __name__ == "__main__":
+    pass

test/difficulty_index/example_difficulty_index.py

@@ -0,0 +1,151 @@
+#!/usr/bin/env conda run -n diff_index_plot python
+# -*- coding: utf-8 -*-
+"""
+Load fieldijn from npz file created with save_ensemble_data.py
+helper function, compute ensemble mean and spread, compute
+difficulty index for a set of thresholds, plot and save the results.
+Created on Tue Mar 17 08:06:07 2020
+Last modified Mon Apr 6 11:30:30 2020
+@author: campbell
+
+Taken from original test_difficulty_index.py but replacing with METcalcpy and METplotpy.
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from metcalcpy.calc_difficulty_index import forecast_difficulty as di
+from metcalcpy.calc_difficulty_index import EPS
+from metcalcpy.piecewise_linear import PiecewiseLinear as plin
+import metplotpy.difficulty_index.mycolormaps as mcmap
+from metplotpy.difficulty_index.plot_difficulty_index import plot_field
+
+
+def load_data(filename):
+    """Load ensemble data from file"""
+    loaded = np.load(filename)
+    lats, lons = (loaded['lats'], loaded['lons'])
+    fieldijn = np.ma.masked_invalid(
+        np.ma.masked_array(
+            data=loaded['data'],
+            mask=loaded['mask']))
+
+    return lats, lons, fieldijn
+
+
+def compute_stats(field):
+    """Compute mean and std dev"""
+    mu = np.mean(field, axis=-1)
+    sigma = np.std(field, axis=-1, ddof=1)
+
+    return mu, sigma
+
+
+def compute_difficulty_index(field, mu, sigma, thresholds):
+    """
+    Compute difficulty index for an ensemble forecast given
+    a set of thresholds, returning a dictionary of fields.
+    """
+    dij = {}
+    for threshold in thresholds:
+        dij[threshold] =\
+            di(sigma, mu, threshold, field, sigma_over_mu_ref=EPS)
+
+    return dij
+
+
+def plot_difficulty_index(dij, lats, lons, thresholds):
+    """
+    Plot the difficulty index for a set of thresholds,
+    returning a dictionary of figures
+    """
+    plt.close('all')
+    myparams = {'figure.figsize': (8, 5),
+                'figure.max_open_warning': 40}
+    plt.rcParams.update(myparams)
+    figs = {}
+    units = 'feet'
+    cmap = mcmap.stoplight()
+    for threshold in thresholds:
+        if np.max(dij[threshold]) <= 1.0:
+            vmax = 1.0
+        else:
+            vmax = 1.5
+        figs[threshold] =\
+            plot_field(dij[threshold],
+                          lats, lons, vmin=0.0, vmax=vmax, cmap=cmap,
+                          xlab='Longitude \u00b0E', ylab='Latitude',
+                          clab='thresh={} {}'.format(threshold, units),
+                          title='Forecast Decision Difficulty Index')
+
+    return figs
+
+
+def save_difficulty_figures(figs, save_thresh, units='feet'):
+    """
+    Save subset of difficulty index figures.
+    """
+    fig_fmt = 'png'
+    fig_basename = './swh_North_Pacific_difficulty_index_'
+    for thresh in save_thresh:
+        thresh_str = '{:.2f}'.format(thresh).replace('.', '_')
+        fig_name = (fig_basename + thresh_str +
+                    '_' + units + '.' + fig_fmt)
+        print('Saving {}...\n'.format(fig_name))
+        figs[thresh].savefig(fig_name, format=fig_fmt)
+
+
+def plot_statistics(mu, sigma, lats, lons, units='feet'):
+    """Plot ensemble mean and spread, returning figure handles"""
+    cmap = mcmap.spectral()
+    mu_fig =\
+        plot_field(mu, lats, lons, cmap=cmap, clab=units,
+                      vmin=0.0, vmax=np.nanmax(mu),
+                      xlab='Longitude \u00b0E',
+                      ylab='Latitude',
+                      title='Forecast Ensemble Mean')
+    sigma_fig =\
+        plot_field(sigma, lats, lons, cmap=cmap, clab=units,
+                      vmin=0.0, vmax=np.nanmax(sigma),
+                      xlab='Longitude \u00b0E',
+                      ylab='Latitude',
+                      title='Forecast Ensemble Std')
+
+    return mu_fig, sigma_fig
+
+
+def save_stats_figures(mu_fig, sigma_fig):
+    """
+    Save ensemble mean and spread figures.
+    """
+    fig_fmt = 'png'
+    fig_basename = './swh_North_Pacific_5dy_'
+    mu_name = fig_basename + 'mean.' + fig_fmt
+    print('Saving {}...\n'.format(mu_name))
+    mu_fig.savefig(mu_name, format=fig_fmt)
+    sigma_name = fig_basename + 'std.' + fig_fmt
+    print('Saving {}...\n'.format(sigma_name))
+    sigma_fig.savefig(sigma_name, format=fig_fmt)
+
+
+def main():
+    """
+    Load fieldijn from npz file created with save_ensemble_data.py
+    helper function, compute ensemble mean and spread, compute
+    difficulty index for a set of thresholds, plot and save the results.
+    """
+    filename = './swh_North_Pacific_5dy_ensemble.npz'
+    lats, lons, fieldijn = load_data(filename)
+    muij, sigmaij = compute_stats(fieldijn)
+    thresholds = np.arange(4.0, 16.0, 1.0)
+    dij = compute_difficulty_index(fieldijn, muij, sigmaij, thresholds)
+    figs = plot_difficulty_index(dij, lats, lons, thresholds)
+    save_thresh = np.arange(9.0, 13.0, 1.0)
+    save_difficulty_figures(figs, save_thresh)
+    units = 'feet'
+    mu_fig, sigma_fig =\
+        plot_statistics(muij, sigmaij, lats, lons, units=units)
+    save_stats_figures(mu_fig, sigma_fig)
+
+
+if __name__ == '__main__':
+    main()

test/difficulty_index/test_difficulty_index_plotting.py

@@ -0,0 +1,19 @@
+import example_difficulty_index as edi  
+from metcalcpy.compare_images import CompareImages
+import pytest
+import warnings
+
+def test_difficulty_index_plot():
+    """
+    Compare the expected image (diff_index_expected.png) to the latest
+    one generated.  Invoke the plot_difficulty_index.py
+    script to generate the difficulty_index.png
+    """
+
+    warnings.filterwarnings("ignore", category=DeprecationWarning)
+
+    file1 = 'swh_North_Pacific_5dy_ensemble.npz'
+    edi.main()
+    comparison = CompareImages('diff_index_expected.png', 'swh_North_Pacific_difficulty_index_9_00_feet.png')
+    assert comparison.mssim == 1
+