[2/2] climate git commit: CLIMATE-850 - Precipitation intensity and duration analyzer using hourly data

[hu...@apache.org]

CLIMATE-850 - Precipitation intensity and duration analyzer using hourly data


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

Branch: refs/heads/master
Commit: 817c854baf52a19d25412da15eecfa6cd9f7abda
Parents: 43b04e4 111d3a4
Author: huikyole <hu...@argo.jpl.nasa.gov>
Authored: Tue Aug 16 08:50:23 2016 -0700
Committer: huikyole <hu...@argo.jpl.nasa.gov>
Committed: Tue Aug 16 08:50:23 2016 -0700

----------------------------------------------------------------------
 ocw/metrics.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/climate/blob/817c854b/ocw/metrics.py
----------------------------------------------------------------------
diff --cc ocw/metrics.py
index 0d79928,e29dc74..8e61799
--- a/ocw/metrics.py
+++ b/ocw/metrics.py
@@@ -347,44 -347,46 +347,87 @@@ def calc_rmse(target_array, reference_a
  
      return (ma.mean((calc_bias(target_array, reference_array))**2))**0.5 
  
 +def calc_histogram_overlap(hist1, hist2):
 +    ''' from Lee et al. (2014)
 +    :param hist1: a histogram array
 +    :type hist1: :class:'numpy.ndarray'
 +    :param hist2: a histogram array with the same size as hist1
 +    :type hist2: :class:'numpy.ndarray'
 +    '''
 +   
 +    hist1_flat = hist1.flatten()
 +    hist2_flat = hist2.flatten()
 +
 +    if len(hist1_flat) != len(hist2_flat):
 +        err = "The two histograms have different sizes"
 +        raise ValueError(err) 
 +    overlap = 0.
 +    for ii in len(hist1_flat):
 +        overlap = overlap + numpy.min(hist1_flat[ii], hist2_flat[ii])
 +    return overlap
 +
 +def calc_joint_histogram(data_array1, data_array2, bins_for_data1, bins_for_data2):
 +    ''' Calculate a joint histogram of two variables in data_array1 and data_array2
 +    :param data_array1: the first variable
 +    :type data_array1: :class:'numpy.ma.core.MaskedArray'
 +    :param data_array2: the second variable
 +    :type data_array2: :class:'numpy.ma.core.MaskedArray'
 +    :param bins_for_data1: histogram bin edges for data_array1
 +    :type bins_for_data1: :class:'numpy.ndarray'
 +    :param bins_for_data2: histogram bin edges for data_array2
 +    :type bins_for_data2: :class:'numpy.ndarray'
 +    '''
 +    if ma.count_masked(data_array1)!=0 or ma.count_masked(data_array2)!=0:
 +        index = numpy.where((data_array1.mask == False) & (data_array2.mask==False)) 
 +        new_array1 = data_array1[index]
 +        new_array2 = data_array2[index] 
 +    else:
 +        new_array1 = data_array1.flatten()
 +        new_array2 = data_array2.flatten()
 +
 +    histo2d, xedge, yedge = numpy.histogram2d(new_array1, new_array2, bins=[bins_for_data1, bins_for_data2])
 +    return histo2d
 + 
+ def wet_spell_analysis(reference_array, threshold=0.1, nyear=1, dt=3.):
+     ''' Characterize wet spells using sub-daily (hourly) data
+ 
+     :param reference_array: an array to be analyzed
+     :type reference_array: :class:'numpy.ma.core.MaskedArray'
+ 
+     :param threshold: the minimum amount of rainfall [mm/hour] 
+     :type threshold: 'float'
+ 
+     :param nyear: the number of discontinous periods 
+     :type nyear: 'int'
+ 
+     :param dt: the temporal resolution of reference_array
+     :type dt: 'float'
+     '''
+     nt = reference_array.shape[0]
+     if reference_array.ndim == 3:
+         reshaped_array = reference_array.reshape[nt, reference_array.size/nt]
+     else:
+         reshaped_array = reference_array
+     xy_indices = np.where(reshaped_array.mask[0,:] == False)[0]
+ 
+     nt_each_year = nt/nyear 
+     spell_duration = []
+     peak_rainfall = []
+     total_rainfall = []
+    
+     for index in xy_indices:
+         for iyear in np.arange(nyear):
+             data0_temp = reshaped_array[nt_each_year*iyear:nt_each_year*(iyear+1),
+                                         index] 
+             # time indices when precipitation rate is smaller than the threshold [mm/hr]
+             t_index = np.where((data0_temp <= threshold) & (data0_temp.mask ==False))[0]
+             t_index = np.insert(t_index, 0, 0)
+             t_index = t_index + nt_each_year*iyear
+             for it in np.arange(t_index.size-1):
+                 if t_index[it+1] - t_index[it] >1:
+                     data1_temp = data0_temp[t_index[it]+1:t_index[it+1]]
+                     if not ma.is_masked(data1_temp):
+                         spell_duration.append((t_index[it+1]-t_index[it]-1)*dt)
+                         peak_rainfall.append(data1_temp.max())
+                         total_rainfall.append(data1_temp.sum())
+     return np.array(spell_duration), np.array(peak_rainfall), np.array(total_rainfall)