You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@climate.apache.org by hu...@apache.org on 2017/04/18 21:54:06 UTC
[1/2] climate git commit: CLIMATE-909 - An example script to evaluate
joint PDF of precipitation intensity and duration
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 <hu...@argo.jpl.nasa.gov>
Authored: Mon Apr 17 21:39:08 2017 -0700
Committer: huikyole <hu...@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')
[2/2] climate git commit: Merge branch 'CLIMATE-909'
Posted by hu...@apache.org.
Merge branch 'CLIMATE-909'
Project: http://git-wip-us.apache.org/repos/asf/climate/repo
Commit: http://git-wip-us.apache.org/repos/asf/climate/commit/03f3541a
Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/03f3541a
Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/03f3541a
Branch: refs/heads/master
Commit: 03f3541a4c75813b87000a453f3a1b8f6692ed9a
Parents: 169f182 2c4b4be
Author: huikyole <hu...@argo.jpl.nasa.gov>
Authored: Tue Apr 18 14:53:53 2017 -0700
Committer: huikyole <hu...@argo.jpl.nasa.gov>
Committed: Tue Apr 18 14:53:53 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(-)
----------------------------------------------------------------------