Merge pull request #1 from lee1043/add_taylor_diagram

Add taylor diagram

pcmdi_metrics/graphics/__init__.py

@@ -12,3 +12,4 @@
 from .portrait_plot.portrait_plot_lib import portrait_plot  # noqa
 from .portrait_plot.portrait_plot_lib import prepare_data  # noqa
 from .bar_chart.lib import BarChart  # noqa
+from .taylor_diagram.taylor_diagram import TaylorDiagram  # noqa

pcmdi_metrics/graphics/share/Metrics.py

@@ -5,7 +5,7 @@
 import numpy as np
 import pandas as pd
 
-from pcmdi_metrics.graphics import read_mean_clim_json_files
+from pcmdi_metrics.graphics import combine_ref_dicts, read_mean_clim_json_files
 
 
 class Metrics:
@@ -54,6 +54,7 @@ def __init__(self, files):
             self.df_dict,
             self.var_list,
             self.var_unit_list,
+            self.var_ref_dict,
             self.regions,
             self.stats,
         ) = read_mean_clim_json_files(files)
@@ -145,6 +146,7 @@ def merge(self, metrics_obj):
 
         result.var_list = list(set(self.var_list + metrics_obj.var_list))
         result.var_unit_list = list(set(self.var_unit_list + metrics_obj.var_unit_list))
+        result.var_ref_dict = combine_ref_dicts(self.var_ref_dict, metrics_obj.var_ref_dict)
         result.regions = list(set(self.regions + metrics_obj.regions))
         result.stats = list(set(self.stats + metrics_obj.stats))
 

pcmdi_metrics/graphics/share/plot_utils.py

@@ -1,5 +1,6 @@
 import os
 import urllib.request
+from collections import defaultdict
 
 import matplotlib.pyplot as plt
 import requests
@@ -74,3 +75,33 @@ def download_archived_results(path, local_dir):
     except Exception:
         print(path, 'not exist in ', url_head_github_repo)
         pass
+
+
+def combine_ref_dicts(d1, d2):
+    """
+    Combine two dictionaries for reference datasets of variables, raising warning message if different reference datasets were used.
+    Below code is revised from https://stackoverflow.com/a/5946322
+
+    Parameters
+    ----------
+    d1: dict
+    d2: dict
+
+    Return
+    ------
+    dd: merged dict
+    """
+    # Merge dicts
+    dd = defaultdict(list)
+    for d in (d1, d2):  # you can list as many input dicts as you want here
+        for key, value in d.items():
+            dd[key].append(value)
+    # Check consistency in content
+    for key in dd:
+        if len(list(set(dd[key]))) == 1:
+            dd[key] = dd[key][0]
+        else:
+            print('Warning: differnt reference datasets detected for ' + key + ': ', dd[key])
+    # Convert outcome to normal dict and return
+    dd = dict(dd)
+    return dd

pcmdi_metrics/graphics/share/read_json_mean_clim.py

@@ -24,33 +24,37 @@ def read_mean_clim_json_files(
                  (Rows: models, Columns: variables (i.e., 2d array)
     - `var_list`: list of string, all variables from JSON files
     - `var_unit_list`: list of string, all variables and its units from JSON files
+    - `var_ref_list`: list of string, list for default reference dataset for each variable
     - `regions`: list of string, regions
     - `stats`: list of string, statistics
     """
     # Read JSON and get unit for each variable
     results_dict = {}  # merged dict by reading all JSON files
     var_list = []
     var_unit_list = []
+    var_ref_dict = {}
     regions_all = []
 
     for json_file in json_list:
         if debug:
             print("json_file:", json_file)
         with open(json_file) as fj:
-            dict_temp = json.load(fj)
-        var = dict_temp["Variable"]["id"]
+            dict_temp = json.load(fj)  # e.g., load contents of precipitation json file
+        var = dict_temp["Variable"]["id"]  # e.g., 'pr'
         if "level" in list(dict_temp["Variable"].keys()):
             var += "-" + str(int(dict_temp["Variable"]["level"] / 100.0))  # Pa to hPa
         results_dict[var] = dict_temp
         unit = extract_unit(var, results_dict[var])
         var_unit = var + " [" + unit + "]"
         var_list.append(var)
         var_unit_list.append(var_unit)
+        var_ref_dict[var] = extract_ref(var, results_dict[var])
         regions_all.extend(extract_region(var, results_dict[var]))
         if stats is None:
             stats = extract_stat(var, results_dict[var])
     if debug:
         print("var_unit_list:", var_unit_list)
+        print("var_ref_dict:", var_ref_dict)
 
     if regions is None:
         regions = list(set(regions_all))  # Remove duplicates
@@ -81,7 +85,7 @@ def read_mean_clim_json_files(
                     results_dict, var_list, region, stat, season, mip, debug
                 )
 
-    return df_dict, var_list, var_unit_list, regions, stats
+    return df_dict, var_list, var_unit_list, var_ref_dict, regions, stats
 
 
 def extract_unit(var, results_dict_var):
@@ -90,6 +94,12 @@ def extract_unit(var, results_dict_var):
     return units
 
 
+def extract_ref(var, results_dict_var):
+    model_list = sorted(list(results_dict_var["RESULTS"].keys()))
+    ref = results_dict_var["RESULTS"][model_list[0]]["default"]["source"]
+    return ref
+
+
 def extract_region(var, results_dict_var):
     region_list, stat_list = extract_region_stat(var, results_dict_var)
     return region_list

pcmdi_metrics/graphics/taylor_diagram/__init__.py

@@ -0,0 +1 @@
+# init file

pcmdi_metrics/graphics/taylor_diagram/taylor_diagram.py

@@ -0,0 +1,136 @@
+def TaylorDiagram(
+        stddev, corrcoef, refstd, fig, colors,
+        normalize=True,
+        labels=None, markers=None, markersizes=None, zorders=None,
+        ref_label=None, smax=None):
+
+    """Plot a Taylor diagram.
+    This code was adpated from the ILAMB code by Nathan Collier found here:
+    https://github.com/rubisco-sfa/ILAMB/blob/master/src/ILAMB/Post.py#L80,
+    which was revised by Jiwoo Lee to enable more flexible customization of the plot.
+
+    The original code was written by Yannick Copin that can be found here:
+    https://gist.github.com/ycopin/3342888
+
+    Parameters
+    ----------
+    stddev : numpy.ndarray
+        an array of standard deviations
+    corrcoeff : numpy.ndarray
+        an array of correlation coefficients
+    refstd : float
+        the reference standard deviation
+    fig : matplotlib figure
+        the matplotlib figure
+    colors : array
+        an array of colors for each element of the input arrays
+    normalize : bool, optional
+        disable to skip normalization of the standard deviation
+    labels : list, optional
+        list of text for labels
+    markers : list, optional
+        list of marker type
+    markersizes : list, optional
+        list of integer for marker size
+    zorders : list, optional
+        list of integer for zorder
+    ref_label : str, optional
+        label for reference data
+    smax : int or float, optional
+        maximum of axis range for (normalized) standard deviation
+
+    Return
+    ------
+    fig : matplotlib figure
+        the matplotlib figure
+    ax : matplotlib axis
+        the matplotlib axis
+    """
+    import mpl_toolkits.axisartist.floating_axes as FA
+    import mpl_toolkits.axisartist.grid_finder as GF
+    import numpy as np
+    from matplotlib.projections import PolarAxes
+
+    # define transform
+    tr = PolarAxes.PolarTransform()
+
+    # correlation labels
+    rlocs = np.concatenate((np.arange(10) / 10., [0.95, 0.99]))
+    tlocs = np.arccos(rlocs)
+    gl1 = GF.FixedLocator(tlocs)
+    tf1 = GF.DictFormatter(dict(zip(tlocs, map(str, rlocs))))
+
+    # standard deviation axis extent
+    if normalize:
+        stddev = stddev / refstd
+        refstd = 1.
+    smin = 0
+    if smax is None:
+        smax = max(2.0, 1.1 * stddev.max())
+
+    # add the curvilinear grid
+    ghelper = FA.GridHelperCurveLinear(tr,
+                                       extremes=(0, np.pi / 2, smin, smax),
+                                       grid_locator1=gl1,
+                                       tick_formatter1=tf1)
+    ax = FA.FloatingSubplot(fig, 111, grid_helper=ghelper)
+    fig.add_subplot(ax)
+
+    # adjust axes
+    ax.axis["top"].set_axis_direction("bottom")
+    ax.axis["top"].toggle(ticklabels=True, label=True)
+    ax.axis["top"].major_ticklabels.set_axis_direction("top")
+    ax.axis["top"].label.set_axis_direction("top")
+    ax.axis["top"].label.set_text("Correlation")
+    ax.axis["left"].set_axis_direction("bottom")
+    if normalize:
+        ax.axis["left"].label.set_text("Normalized standard deviation")
+    else:
+        ax.axis["left"].label.set_text("Standard deviation")
+    ax.axis["right"].set_axis_direction("top")
+    ax.axis["right"].toggle(ticklabels=True)
+    ax.axis["right"].major_ticklabels.set_axis_direction("left")
+    ax.axis["bottom"].set_visible(False)
+    ax.grid(True)
+
+    ax = ax.get_aux_axes(tr)
+    # Plot data
+    corrcoef = corrcoef.clip(-1, 1)
+    for i in range(len(corrcoef)):
+        # customize label
+        if labels is not None:
+            label = labels[i]
+        else:
+            label = None
+        # customize marker
+        if markers is not None:
+            marker = markers[i]
+        else:
+            marker = 'o'
+        # customize marker size
+        if markersizes is not None:
+            ms = markersizes[i]
+        else:
+            ms = 8
+        # customize marker order
+        if zorders is not None:
+            zorder = zorders[i]
+        else:
+            zorder = None
+        # customize end
+        ax.plot(np.arccos(corrcoef[i]), stddev[i], marker, color=colors[i], mew=0, ms=ms, label=label, zorder=zorder)
+
+    # Add reference point and stddev contour
+    l, = ax.plot([0], refstd, 'k*', ms=12, mew=0, label=ref_label)
+    t = np.linspace(0, np.pi / 2)
+    r = np.zeros_like(t) + refstd
+    ax.plot(t, r, 'k--')
+
+    # centralized rms contours
+    rs, ts = np.meshgrid(np.linspace(smin, smax),
+                         np.linspace(0, np.pi / 2))
+    rms = np.sqrt(refstd**2 + rs**2 - 2 * refstd * rs * np.cos(ts))
+    contours = ax.contour(ts, rs, rms, 5, colors='k', alpha=0.4)
+    ax.clabel(contours, fmt='%1.1f')
+
+    return fig, ax