How to adjust the cumulative period and weigh in python code(effective drought index(EDI:BYUN, WILHITE 1999)

138 views Asked by At

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

  1. im using the python iris and numpy etc.,

  2. EP_cube.data or MEP_cube.data . ".data" is arrays

  3. 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
0

There are 0 answers