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 <huiky...@argo.jpl.nasa.gov> Authored: Tue Aug 16 08:50:23 2016 -0700 Committer: huikyole <huiky...@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)
