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 <huiky...@argo.jpl.nasa.gov> Authored: Mon Apr 17 21:39:08 2017 -0700 Committer: huikyole <huiky...@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')
