I'm make a code, but big problem,, It takes too long to calculate.
In short, the definition of this EDI(Byun and wilhite 1999) EDI is accumulate and weighting precipitation index
pr = precipitation, accumulate days : 365 days(1 years) after accumulate days is redifined.
Big problem is redifined code..
my code is belows Main code: return_redifine_EP_MEP_DS
Sub code : return_EP_v2, return_MEP_v2 , return_DEP_v2
Additional explanation
im using the python iris and numpy etc.,
EP_cube.data or MEP_cube.data . ".data" is arrays
Data sets. CMIP5 models precipitation(time, latitude, longitude)
origin_pr is precipitation and shape (time, latitude, longitude): 1971~2000
EP_cube is accumulate precipitation and (time, latitude, longitude): 1972~2000
MEP_cube is clim mean(each calendar day)of EP_cube(time, latitude, longitude): 365 days
DEP_cube is EP_cube minus clim mean EP_cube(MEP) , (time, latitude, longitud): 1972~2000
sy,ey is climatology years
day is each models 1years day(example: HadGEM2-AO is 360days, ACCESS1-0: 365 days)
def return_redifine_EP_MEP_DS(origin_pr,EP_cube,DEP_cube,sy,ey,day):
origin_DS = day
origin_day= day
DS_cube = EP_cube - EP_cube.data + day
DS = origin_DS
o_DEP_cube =DEP_cube.copy()
return_time = 0
for j in range(0,DEP_cube.data.shape[1]):
for k in range(0,DEP_cube.data.shape[2]):
for i in range(1,DEP_cube.data.shape[0]):
if DEP_cube.data[i,j,k] <0 and DEP_cube.data[i-1,j,k] < 0:
DS = DS + 1
else:
DS = origin_DS
if DS != origin_DS:
EP_cube.data[i,j,k] = return_EP_v2(origin_pr[i-DS+origin_day:i+origin_day,j,k], DS)
day_of_year = DEP_cube[i,j,k].coord('day_of_year').points
MEP_cube.data[day_of_year-1,j,k] = return_MEP_v2(EP_cube[:,j,k], sy, ey,day_of_year)
DEP_cube.data[i,j,k] = return_DEP_v2(EP_cube[i,j,k], MEP_cube[day_of_year-1,j,k])
return EP_cube, MEP_cube, DEP_cube, DS_cube
and sub code below
def return_EP_v2(origin_pr,DS):
weights = np.arange(DS,0,-1) ## 1/1, 1/2, 1/3 ....
EP_data = np.sum(origin_pr.data / weights)
return EP_data
def return_MEP_v2(EP_cube,sy,ey,day_of_year):
### Below mean is extract years and want days(julian day)
ex_EP_cube = EP_cube.extract(iris.Constraint(year = lambda t: sy<=t<=ey, day_of_year = day_of_year ))
### Below mean is clim mean of EP_cube
MEP_data = ex_EP_cube.collapsed('day_of_year',iris.analysis.MEAN).data
return MEP_data
def return_DEP_v2(EP_cube, MEP_cube):
DEP_data = EP_cube.data - MEP_cube.data
return DEP_data