climate-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From huiky...@apache.org
Subject [1/2] climate git commit: CLIMATE-909 - An example script to evaluate joint PDF of precipitation intensity and duration
Date Tue, 18 Apr 2017 21:54:06 GMT
Repository: climate
Updated Branches:
  refs/heads/master 169f182b2 -> 03f3541a4


CLIMATE-909 - An example script to evaluate joint PDF of precipitation intensity and duration

- A new example, GPM_WRF24_JPDF_comparison, has been added.
- Joint PDF metric has been debugged.
- A new plotting module, draw_precipitation_JPDF, has been added.


Project: http://git-wip-us.apache.org/repos/asf/climate/repo
Commit: http://git-wip-us.apache.org/repos/asf/climate/commit/2c4b4bec
Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/2c4b4bec
Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/2c4b4bec

Branch: refs/heads/master
Commit: 2c4b4becffccad1ee8dbc5eb6681124ef046bb92
Parents: 169f182
Author: huikyole <huikyole@argo.jpl.nasa.gov>
Authored: Mon Apr 17 21:39:08 2017 -0700
Committer: huikyole <huikyole@argo.jpl.nasa.gov>
Committed: Mon Apr 17 21:39:08 2017 -0700

----------------------------------------------------------------------
 examples/GPM_WRF24_JPDF_comparison.py | 103 +++++++++++++++++++++++++++++
 ocw/data_source/local.py              |   2 +-
 ocw/metrics.py                        |   5 +-
 ocw/plotter.py                        |  68 +++++++++++++++++++
 4 files changed, 176 insertions(+), 2 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/climate/blob/2c4b4bec/examples/GPM_WRF24_JPDF_comparison.py
----------------------------------------------------------------------
diff --git a/examples/GPM_WRF24_JPDF_comparison.py b/examples/GPM_WRF24_JPDF_comparison.py
new file mode 100644
index 0000000..20b070e
--- /dev/null
+++ b/examples/GPM_WRF24_JPDF_comparison.py
@@ -0,0 +1,103 @@
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#    http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+
+from ocw.dataset import Bounds
+import ocw.data_source.local as local
+import ocw.dataset_processor as dsp
+import ocw.metrics as metrics
+import ocw.plotter as plotter
+
+import numpy as np
+import numpy.ma as ma
+from pickle import load, dump
+
+"""
+This is an example of calculating the joint probability distribution function of rainfall
intensity and duration for the Northern Great Plains using GPM IMERG data for June/01/2015
+""" 
+
+""" Step 1: Load the GPM and WRF24 datasets with spatial filter """
+
+GPM_dataset_filtered = local.load_GPM_IMERG_files_with_spatial_filter(
+                           file_path='./data/GPM_2015_summer/', 
+                           filename_pattern=['*2015*.HDF5'],
+                           user_mask_file='Bukovsky_regions.nc',
+                           mask_variable_name='Bukovsky',
+                           user_mask_values=[10],
+                           longitude_name='lon',
+                           latitude_name='lat')
+
+WRF_dataset = local.load_WRF_2d_files_RAIN(
+                  file_path='./data/WRF24_2010_summer/', 
+                  filename_pattern=['wrf2dout*'])
+
+""" Step 2: Load the spatial filter (Bukovsky region mask) """
+
+Bukovsky_mask = Bounds(boundary_type='user',
+                       user_mask_file='Bukovsky_regions.nc',
+                       mask_variable_name='Bukovsky',
+                       longitude_name='lon',
+                       latitude_name='lat')  
+
+""" Step 3: Spatial subset the WRF data (for Northern Great Plains, user_mask_values=[10])
"""
+
+WRF_dataset_filtered = dsp.subset(WRF_dataset, 
+                           Bukovsky_mask, user_mask_values=[10])
+
+""" Step 4: Analyze the wet spells """
+duration1, peak1, total1 = metrics.wet_spell_analysis(GPM_dataset_filtered, 
+                               threshold=0.1, nyear=1, dt=0.5)
+duration2, peak2, total2 = metrics.wet_spell_analysis(WRF_dataset_filtered.values, 
+                               threshold=0.1, nyear=1, dt=0.5)
+
+""" Step 5: Calculate the joint PDF(JPDF) of spell_duration and peak_rainfall """
+
+histo2d_GPM = metrics.calc_joint_histogram(data_array1=duration1, data_array2=peak1, 
+                   bins_for_data1=np.append(np.arange(25)+0.5,[48.5,120.5]),
+                   bins_for_data2=[0.1,0.2,0.5,1.0,2.0,5.0,10,20,30]) 
+histo2d_GPM = histo2d_GPM/np.sum(histo2d_GPM)*100.
+
+histo2d_WRF = metrics.calc_joint_histogram(data_array1=duration2, data_array2=peak2, 
+                   bins_for_data1=np.append(np.arange(25)+0.5,[48.5,120.5]),
+                   bins_for_data2=[0.1,0.2,0.5,1.0,2.0,5.0,10,20,30]) 
+histo2d_WRF = histo2d_WRF/np.sum(histo2d_WRF)*100.
+
+""" Step 6: Visualize the JPDF """
+
+
+plot_level = np.array([0., 0.01,0.05, 0.1, 0.2, 0.5,1,2,3,5,10,25])
+plotter.draw_precipitation_JPDF(plot_data=np.transpose(histo2d_GPM), 
+          plot_level=plot_level, title='', 
+          x_ticks=[0.5, 4.5, 9.5, 14.5, 19.5, 23.5], x_names=['1','5','10','15','20','24'],
+          y_ticks=np.arange(9), y_names=['0.1','0.2','0.5','1.0','2.0','5.0','10','20','30'],

+          output_file='GPM_JPDF_example')
+plotter.draw_precipitation_JPDF(plot_data=np.transpose(histo2d_WRF), 
+          plot_level=plot_level, title='', 
+          x_ticks=[0.5, 4.5, 9.5, 14.5, 19.5, 23.5], x_names=['1','5','10','15','20','24'],
+          y_ticks=np.arange(9), y_names=['0.1','0.2','0.5','1.0','2.0','5.0','10','20','30'],

+          output_file='WRF24_JPDF_example')
+
+overlap = metrics.calc_histogram_overlap(histo2d_GPM, histo2d_WRF)
+plot_level = np.array([-21, -3, -1, -0.5, -0.2, -0.1, -0.05, 0, 0.05, 0.1, 0.2, 0.5, 1, 3,
21])
+cbar_ticks = [-2, -0.5, -0.1, 0.1, 0.5, 2]
+cbar_label = [str(i) for i in cbar_ticks]
+plotter.draw_precipitation_JPDF(plot_data=np.transpose(histo2d_WRF - histo2d_GPM), 
+          plot_level=plot_level, title='overlap %d ' %overlap +r'%', diff_plot=True,
+          x_ticks=[0.5, 4.5, 9.5, 14.5, 19.5, 23.5], x_names=['1','5','10','15','20','24'],
+          y_ticks=np.arange(9), y_names=['0.1','0.2','0.5','1.0','2.0','5.0','10','20','30'],
+          output_file='GPM_WRF24_JPDF_comparison',  
+          cbar_ticks=cbar_ticks, cbar_label=cbar_label)
+

http://git-wip-us.apache.org/repos/asf/climate/blob/2c4b4bec/ocw/data_source/local.py
----------------------------------------------------------------------
diff --git a/ocw/data_source/local.py b/ocw/data_source/local.py
index 78cf3cc..b494d84 100644
--- a/ocw/data_source/local.py
+++ b/ocw/data_source/local.py
@@ -739,7 +739,7 @@ def load_GPM_IMERG_files_with_spatial_filter(file_path=None,
         file_object = h5py.File(file)
         values0 = ma.transpose(ma.masked_less(
             file_object['Grid'][variable_name][:], 0.))
-        values_masked = values0[y_index, x_index]
+        values_masked = values0[y_index, x_index] 
         values_masked = ma.expand_dims(values_masked, axis=0)
         if ifile == 0:
             values = values_masked

http://git-wip-us.apache.org/repos/asf/climate/blob/2c4b4bec/ocw/metrics.py
----------------------------------------------------------------------
diff --git a/ocw/metrics.py b/ocw/metrics.py
index cd027f0..3a6477e 100644
--- a/ocw/metrics.py
+++ b/ocw/metrics.py
@@ -421,7 +421,10 @@ def wet_spell_analysis(reference_array, threshold=0.1, nyear=1, dt=3.):
         reshaped_array = reference_array.reshape([nt, reference_array.size / nt])
     else:
         reshaped_array = reference_array
-    xy_indices = numpy.where(reshaped_array.mask[0, :] == False)[0]
+    if ma.count_masked(reshaped_array[0,:]) != 0:
+        xy_indices = numpy.where(reshaped_array.mask[0, :] == False)[0]
+    else:
+        xy_indices = numpy.arange(reshaped_array.shape[1])
 
     nt_each_year = nt / nyear
     spell_duration = []

http://git-wip-us.apache.org/repos/asf/climate/blob/2c4b4bec/ocw/plotter.py
----------------------------------------------------------------------
diff --git a/ocw/plotter.py b/ocw/plotter.py
index 55302a9..4896a64 100755
--- a/ocw/plotter.py
+++ b/ocw/plotter.py
@@ -18,6 +18,8 @@
 from tempfile import TemporaryFile
 import matplotlib as mpl
 import matplotlib.pyplot as plt
+from matplotlib.colors import BoundaryNorm
+from matplotlib import rcParams, cm
 from matplotlib.patches import Polygon
 from matplotlib.collections import PatchCollection
 from mpl_toolkits.basemap import Basemap
@@ -1206,3 +1208,69 @@ def draw_plot_to_compare_trends(obs_data, ens_data, model_data,
     if data_labels:
         ax.set_xticklabels(data_labels)
     fig.savefig('%s.%s' % (fname, fmt), bbox_inches='tight') 
+
+def draw_precipitation_JPDF (plot_data, plot_level, x_ticks, x_names,y_ticks,y_names, 
+               output_file, title=None, diff_plot=False, cmap = cm.BrBG, 
+               cbar_ticks=[0.01, 0.10, 0.5, 2, 5, 25], 
+               cbar_label=['0.01', '0.10', '0.5', '2', '5', '25']):
+    '''
+    :param plot_data: a numpy array of data to plot (dimY, dimX)
+    :type plot_data: :class:'numpy.ndarray'
+    :param plot_level: levels to plot
+    :type plot_level: :class:'numpy.ndarray'
+    :param x_ticks: x values where tick makrs are located
+    :type x_ticks: :class:'numpy.ndarray'                 
+    :param x_names: labels for the ticks on x-axis (dimX)
+    :type x_names: :class:'list'
+    :param y_ticks: y values where tick makrs are located
+    :type y_ticks: :class:'numpy.ndarray'                 
+    :param y_names: labels for the ticks on y-axis (dimY)
+    :type y_names: :class:'list'
+    :param output_file: name of output png file
+    :type output_file: :mod:'string'
+    :param title: (Optional) title of the plot
+    :type title: :mod:'string'
+    :param diff_plot: (Optional) if true, a difference plot will be generated
+    :type diff_plot: :mod:'bool'
+    :param cbar_ticks: (Optional) tick marks for the color bar
+    :type cbar_ticks: :class:'list'
+    :param cbar_label: (Optional) lables for the tick marks of the color bar
+    :type cbar_label: :class:'list'
+    '''
+
+    if diff_plot:
+        cmap = cm.RdBu_r
+
+    fig = plt.figure()
+    sb = fig.add_subplot(111)
+    dimY, dimX = plot_data.shape
+    plot_data2 = np.zeros([dimY,dimX]) # sectioned array for plotting
+
+    # fontsize 
+    rcParams['axes.labelsize'] = 8
+    rcParams['xtick.labelsize'] = 8
+    rcParams['ytick.labelsize'] = 8
+    # assign the values in plot_level to plot_data
+    for iy in range(dimY):
+         for ix in range(dimX):
+             if plot_data[iy,ix] <= np.min(plot_level):
+                 plot_data2[iy,ix] = -1.
+             else:
+                 plot_data2[iy,ix] = plot_level[np.where(plot_level <= plot_data[iy,ix])].max()
+    sb.set_xticks(x_ticks)
+    sb.set_xticklabels(x_names)
+    sb.set_yticks(y_ticks)
+    sb.set_yticklabels(y_names)
+
+    norm = BoundaryNorm(plot_level[1:], cmap.N)
+    a=sb.pcolor(plot_data2, edgecolors = 'k', linewidths=0.5, cmap = cmap, norm = norm)
+    a.cmap.set_under('w')
+    sb.set_xlabel('Spell duration [hrs]')
+    sb.set_ylabel('Peak rainfall [mm/hr]')
+    cax = fig.add_axes([0.2, -0.06, 0.6, 0.02])
+    cbar = plt.colorbar(a, cax=cax, orientation = 'horizontal', extend='both')
+    cbar.set_ticks(cbar_ticks)
+    cbar.set_ticklabels(cbar_label)
+    if title:
+        fig.suptitle(title)
+    fig.savefig(output_file, dpi=600,bbox_inches='tight')


Mime
View raw message