diff --git a/.github/parm/use_case_groups.json b/.github/parm/use_case_groups.json index 119576e30..30822690d 100644 --- a/.github/parm/use_case_groups.json +++ b/.github/parm/use_case_groups.json @@ -127,7 +127,7 @@ { "category": "s2s", "index_list": "1-3", - "run": false + "run": true }, { "category": "s2s", @@ -157,7 +157,7 @@ { "category": "s2s", "index_list": "11", - "run": false + "run": true }, { "category": "s2s", diff --git a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_Blocking/Blocking.py b/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_Blocking/Blocking.py deleted file mode 120000 index 29f05d48e..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_Blocking/Blocking.py +++ /dev/null @@ -1 +0,0 @@ -../UserScript_obsERA_obsOnly_Blocking/Blocking.py \ No newline at end of file diff --git a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_Blocking/Blocking_WeatherRegime_util.py b/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_Blocking/Blocking_WeatherRegime_util.py deleted file mode 120000 index 7839a8807..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_Blocking/Blocking_WeatherRegime_util.py +++ /dev/null @@ -1 +0,0 @@ -../UserScript_obsERA_obsOnly_Blocking/Blocking_WeatherRegime_util.py \ No newline at end of file diff --git a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime.conf b/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime.conf index ab36a19d6..d06418294 100644 --- a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime.conf +++ b/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime.conf @@ -2,9 +2,8 @@ [config] # All steps, including pre-processing: # PROCESS_LIST = RegridDataPlane(regrid_obs), PcpCombine(daily_mean_obs), UserScript(script_wr), StatAnalysis(sanal_wrclass), StatAnalysis(sanal_wrfreq) -# Weather Regime Analysis only: -#PROCESS_LIST = UserScript(script_wr), StatAnalysis(sanal_wrclass), StatAnalysis(sanal_wrfreq) -PROCESS_LIST = UserScript(script_wr) +# Weather Regime Analysis and stat_analysis: +PROCESS_LIST = UserScript(script_wr), StatAnalysis(sanal_wrclass), StatAnalysis(sanal_wrfreq) # time looping - options are INIT, VALID, RETRO, and REALTIME # If set to INIT or RETRO: @@ -178,6 +177,13 @@ FCST_NUM_PCS = {OBS_NUM_PCS} OBS_WR_FREQ = 7 FCST_WR_FREQ = {OBS_WR_FREQ} +# These variables control reordering the forecast weather regime to match the +# observations if their orders are different +# It is recommended to set this to False if this is the first time running the +# case +REORDER_FCST = True +FCST_ORDER = 1,3,4,2,5,6 + # Type, name and directory of Output File for weather regime classification # Type options are text or netcdf OBS_WR_OUTPUT_FILE_TYPE = text diff --git a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime/Blocking_WeatherRegime_util.py b/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime/Blocking_WeatherRegime_util.py deleted file mode 120000 index 7839a8807..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime/Blocking_WeatherRegime_util.py +++ /dev/null @@ -1 +0,0 @@ -../UserScript_obsERA_obsOnly_Blocking/Blocking_WeatherRegime_util.py \ No newline at end of file diff --git a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime/WeatherRegime.py b/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime/WeatherRegime.py deleted file mode 120000 index a9f379087..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_fcstGFS_obsERA_WeatherRegime/WeatherRegime.py +++ /dev/null @@ -1 +0,0 @@ -../UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime.py \ No newline at end of file diff --git a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking.py b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking.py deleted file mode 100644 index 9dd2efc62..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking.py +++ /dev/null @@ -1,414 +0,0 @@ -import os -import numpy as np -import datetime -import bisect -from scipy import stats -from scipy.signal import argrelextrema -#from metplus.util import config_metplus, get_start_end_interval_times, get_lead_sequence -#from metplus.util import get_skip_times, skip_time, is_loop_by_init, ti_calculate -from Blocking_WeatherRegime_util import read_nc_met - -class BlockingCalculation(): - """Contains the programs to calculate Blocking via the Pelly-Hoskins Method - """ - def __init__(self,label): - - self.blocking_anomaly_var = os.environ.get(label+'_BLOCKING_ANOMALY_VAR','Z500_ANA') - self.blocking_var = os.environ.get(label+'_BLOCKING_VAR','Z500') - self.smoothing_pts = int(os.environ.get(label+'_SMOOTHING_PTS',9)) - lat_delta_in = os.environ.get(label+'_LAT_DELTA','-5,0,5') - self.lat_delta = list(map(int,lat_delta_in.split(","))) - self.n_s_limits = int(os.environ.get(label+'_NORTH_SOUTH_LIMITS',30)) - self.ibl_dist = int(os.environ.get(label+'_IBL_DIST',7)) - self.ibl_in_gibl = int(os.environ.get(label+'_IBL_IN_GIBL',15)) - self.gibl_overlap = int(os.environ.get(label+'_GIBL_OVERLAP',10)) - self.block_time = int(os.environ.get(label+'_BLOCK_TIME',5)) ###Should fix so it supports other times" - self.block_travel = int(os.environ.get(label+'_BLOCK_TRAVEL',45)) - self.block_method = os.environ.get(label+'_BLOCK_METHOD','PH') - - # Check data requirements - if self.smoothing_pts % 2 == 0: - print('ERROR: Smoothing Radius must be an odd number given in grid points, Exiting...') - exit() - - - def run_CBL(self,cblinfiles,nseasons,dseasons): - - z500_anom_4d,lats,lons,timedict = read_nc_met(cblinfiles,self.blocking_anomaly_var,nseasons,dseasons) - - #Create Latitude Weight based for NH - cos = lats - cos = cos * np.pi / 180.0 - way = np.cos(cos) - way = np.sqrt(way) - weightf = np.repeat(way[:,np.newaxis],360,axis=1) - - ####Find latitude of maximum high-pass STD (CBL) - yrlen = len(z500_anom_4d[:,0,0,0]) - mstd = np.nanstd(z500_anom_4d,axis=1) - mhweight = mstd * weightf - cbli = np.argmax(mhweight,axis=1) - CBL = np.zeros((yrlen,len(lons))) - for j in np.arange(0,yrlen,1): - CBL[j,:] = lats[cbli[j,:]] - - ###Apply Moving Average to Smooth CBL Profiles - lt = len(lons) - CBLf = np.zeros((yrlen,len(lons))) - m=int((self.smoothing_pts-1)/2.0) - for i in np.arange(0,len(CBL[0,:]),1): - ma_indices = np.arange(i-m,i+m+1) - ma_indices = np.where(ma_indices >= lt,ma_indices-lt,ma_indices) - CBLf[:,i] = np.nanmean(CBL[:,ma_indices],axis=1).astype(int) - - return CBLf,lats,lons,mhweight,timedict - - - def run_mod_blocking1d(self,a,cbl,lat,lon,meth): - lat_d = self.lat_delta - dc = (90 - cbl).astype(int) - db = self.n_s_limits - BI = np.zeros((len(a[:,0,0]),len(lon))) - blon = np.zeros((len(a[:,0,0]),len(lon))) - if meth=='PH': - # loop through days - for k in np.arange(0,len(a[:,0,0]),1): - blontemp=0 - q=0 - BI1=np.zeros((len(lat_d),len(lon))) - for l in lat_d: - blon1 = np.zeros(len(lon)) - d0 = dc-l - dn = ((dc - 1*db/2) - l).astype(np.int64) - ds = ((dc + 1*db/2) - l).astype(np.int64) - GHGS = np.zeros(len(cbl)) - GHGN = np.zeros(len(cbl)) - for jj in np.arange(0,len(cbl),1): - GHGN[jj] = np.mean(a[k,dn[jj]:d0[jj]+1,jj]) - GHGS[jj] = np.mean(a[k,d0[jj]:ds[jj]+1,jj]) - BI1[q,:] = GHGN-GHGS - q = q +1 - BI1 = np.max(BI1,axis=0) - block = np.where((BI1>0))[0] - blon1[block]=1 - blontemp = blontemp + blon1 - BI[k,:] = BI1 - blon[k,:] = blontemp - - return blon,BI - - - def run_Calc_IBL(self,cbl,iblinfiles,nseasons, dseasons): - - z500_daily,lats,lons,timedict = read_nc_met(iblinfiles,self.blocking_var,nseasons,dseasons) - - #Initilize arrays for IBLs and the blocking index - # yr, day, lon - yrlen = len(z500_daily[:,0,0,0]) - blonlong = np.zeros((yrlen,len(z500_daily[0,:,0,0]),len(lons))) - BI = np.zeros((yrlen,len(z500_daily[0,:,0,0]),len(lons))) - - #Using long-term mean CBL and acsessing module of mod.py - cbl = np.nanmean(cbl,axis=0) - for i in np.arange(0,yrlen,1): - blon,BI[i,:,:] = self.run_mod_blocking1d(z500_daily[i,:,:,:],cbl,lats,lons,self.block_method) - blonlong[i,:,:] = blon - - return blonlong,timedict - - - def run_Calc_GIBL(self,ibl,lons): - - #Initilize GIBL Array - GIBL = np.zeros(np.shape(ibl)) - - #####Loop finds IBLs within 7 degree of each other creating one group. Finally - ##### A GIBL exist if it has more than 15 IBLs - crit = self.ibl_in_gibl - - for i in np.arange(0,len(GIBL[:,0,0]),1): - for k in np.arange(0,len(GIBL[0,:,0]),1): - gibli = np.zeros(len(GIBL[0,0,:])) - thresh = crit/2.0 - a = ibl[i,k,:] - db = self.ibl_dist - ibli = np.where(a==1)[0] - if len(ibli)==0: - continue - idiff = ibli[1:] - ibli[:-1] - bt=0 - btlon = ibli[0] - ct = 1 - btfin = [] - block = ibli - block = np.append(block,block+360) - for ll in np.arange(1,len(block),1): - diff = np.abs(block[ll] - block[ll-1]) - if diff == 1: - bt = [block[ll]] - btlon = np.append(btlon,bt) - ct = ct + diff - if diff <= thresh and diff != 1: - bt = np.arange(block[ll-1]+1,block[ll]+1,1) - btlon = np.append(btlon,bt) - ct = ct + diff - if diff > thresh or ll==(len(block)-1): - if ct >= crit: - btfin = np.append(btfin,btlon) - ct=1 - ct = 1 - btlon = block[ll] - if len(btfin)/2 < crit : - btfin = [] - if len(btfin)==0: - continue - gibl1 = btfin - temp = np.where(gibl1>=360)[0] - gibl1[temp] = gibl1[temp] - 360 - gibli[gibl1.astype(int)] = 1 - GIBL[i,k,:] = gibli - - return GIBL - - - def Remove(self,duplicate): - final_list = [] - for num in duplicate: - if num not in final_list: - final_list.append(num) - return final_list - - - def run_Calc_Blocks(self,ibl,GIBL,lon,tsin): - - crit = self.ibl_in_gibl - - ##Count up the blocked longitudes for each GIBL - c = np.zeros((GIBL.shape)) - lonlen = len(lon) - sz = [] - mx = [] - min = [] - - for y in np.arange(0,len(GIBL[:,0,0]),1): - for k in np.arange(0,len(GIBL[0,:,0]),1): - a = GIBL[y,k] # Array of lons for each year,day - ct=1 - ai = np.where(a==1)[0] - ai = np.append(ai,ai+360) - temp = np.where(ai>=360)[0] - bi=list(ai) - bi = np.array(bi) - bi[temp] = bi[temp]-360 - # Loop through the lons that are part of the GIBL - for i in np.arange(0,len(ai)-1,1): - diff = ai[i+1] - ai[i] - c[y,k,bi[i]] = ct - if diff==1: - ct=ct+1 - else: - sz = np.append(sz,ct) - ct=1 - - ########## - finding where the left and right limits of the block are - ################ - for i in np.arange(0,len(c[:,0,0]),1): - for k in np.arange(0,len(c[0,:,0]),1): - maxi = argrelextrema(c[i,k],np.greater,mode='wrap')[0] - mini = np.where(c[i,k]==1)[0] - if c[i,k,lonlen-1]!=0 and c[i,k,0]!=0: - mm1 = mini[-1] - mm2 = mini[:-1] - mini = np.append(mm1,mm2) - mx = np.append(mx,maxi) - min = np.append(min,mini) - - locy, locd, locl = np.where(c==crit) - - A = np.zeros(lonlen) - A = np.expand_dims(A,axis=0) - - ################# - Splitting up each GIBL into its own array - ################### - - for i in np.arange(0,len(locy),1): - m = locy[i] #year - n = locd[i] #day - o = locl[i] #long where 15 - mm = int(mx[i]) - mn = min[i] - temp1 = GIBL[m,n] - temp2 = np.zeros(lonlen) - if mn>mm: - diff = int(mm - c[m,n,mm] + 1) - lons = lon[diff] - place1 = np.arange(lons,lonlen,1) - place2 = np.arange(0,mm+1,1) - bl = np.append(place2,place1).astype(int) - if temp1[lonlen-1] ==1 and mm>200: - lons = lon[mm] - beg = mm - c[m,n,mm] + 1 - bl = np.arange(beg,mm+1,1).astype(int) - if mm>mn: #temp1[359] ==0: - lons = lon[mm] - beg = mm - c[m,n,mm] + 1 - bl = np.arange(beg,mm+1,1).astype(int) - temp2[bl] = 1 - temp2 = np.expand_dims(temp2,axis=0) - A = np.concatenate((A,temp2),axis=0) - - A = A[1:] - - ######### - Getting rid of non-consectutve Time steps which would prevent blocking - ################# - dd=[] - dy = [] - dA = A[0] - dA = np.expand_dims(dA,axis=0) - ct=0 - for i in np.arange(1,len(locy),1): - dd1 = locd[i-1] - dd2 = locd[i] - if dd2-dd1 > 2: - ct = 0 - continue - if ct == 0: - dd = np.append(dd,locd[i-1]) - dy = np.append(dy,locy[i-1]) - temp2 = np.expand_dims(A[i-1],axis=0) - dA = np.concatenate((dA,temp2),axis=0) - ct = ct + 1 - if dd2-dd1<=2: - dd=np.append(dd,locd[i]) - dy = np.append(dy,locy[i]) - temp2 = np.expand_dims(A[i],axis=0) - dA = np.concatenate((dA,temp2),axis=0) - ct = ct + 1 - - dA=dA[1:] - dAfin = dA - - ############ - Finding center longitude of block - ############## - middle=[] - for l in np.arange(0,len(dAfin),1): - temp = np.where(dAfin[l]==1)[0] - if len(temp) % 2 == 0: - temp = np.append(temp,0) - midtemp = np.median(temp) - middle = np.append(middle,midtemp) - - - #####Track blocks. Makes sure that blocks overlap with at least 10 longitude points on consecutive - overlap = self.gibl_overlap - btime = self.block_time - fin = [[]] - finloc = [[]] - ddcopy=dd*1.0 - noloc=[] - failloc = [[]] - for i in np.arange(0,len(c[:,0,0]),1): - yri = np.where(dy==i)[0] - B = [[]] - ddil =1.0 * ddcopy[yri] - dyy = np.where(dy==i)[0] - rem = [] - for dk in ddil: - if len(ddil) < btime: - continue - ddil = np.append(ddil[0]-1,ddil) - diff = np.diff(ddil) - diffB=[] - dB =1 - cnt = 1 - while dB<=2: - diffB = np.append(diffB,ddil[cnt]) - dB = diff[cnt-1] - if ddil[cnt]==ddil[-1]: - dB=5 - cnt=cnt+1 - diffB = np.array(self.Remove(diffB)) - locb = [] - for ll in diffB: - locb = np.append(locb,np.where((dy==i) & (dd==ll))[0]) - ddil=ddil[1:] - locbtemp = 1.0*locb - ree=np.empty(0,dtype=int) - for hh in np.arange(0,len(noloc),1): - ree = np.append(ree,np.where(locb == noloc[hh])[0]) - ree.astype(int) - locbtemp = np.delete(locbtemp,ree) - locb=locbtemp * 1.0 - datemp = dAfin[locb.astype(int)] - blocktemp = [[datemp[0]]] - locbi = np.array([locb[0]]) - ll1=0 - pass1 = 0 - ai=[0] - add=0 - for ll in np.arange(0,len(locb)-1,1): - if ((dd[locb[ll+1].astype(int)] - dd[locb[ll1].astype(int)]) >=1) & ((dd[locb[ll+1].astype(int)] - dd[locb[ll1].astype(int)]) <=2): - add = datemp[ll1] + datemp[ll+1] - ai = np.where(add==2)[0] - if len(ai)>overlap: - locbi=np.append(locbi,locb[ll+1]) - ll1=ll+1 - add=0 - if (len(ai)4: - noloc = np.append(noloc,locbi) - finloc = finloc + [locbi] - for jj in locbi: - rem = np.append(rem,np.where(dyy==jj)[0]) - ddil = np.delete(ddcopy[yri],rem.astype(int)) - if len(locbi)<=4: - noloc = np.append(noloc,locbi) - if len(locbi)<=2: - failloc=failloc+[locbi] - for jj in locbi: - rem = np.append(rem,np.where(dyy==jj)[0]) - ddil = np.delete(ddcopy[yri],rem.astype(int)) - - blocks = finloc[1:] - noblock = failloc[1:] - - ############Get rid of blocks that travel downstream########## - ######If center of blocks travel downstream further than 45 degrees longitude, - ######cancel block moment it travels out of this limit - newblock = [[]] - newnoblock=[[]] - distthresh = self.block_travel - for bb in blocks: - diffb = [] - start = middle[bb[0].astype(int)] - for bbs in bb: - diffb = np.append(diffb, start - middle[bbs.astype(int)]) - loc = np.where(np.abs(diffb) > distthresh)[0] - if len(loc)==0: - newblock = newblock +[bb] - if len(loc)>0: - if len(bb[:loc[0]]) >4: - newblock = newblock + [bb[:loc[0]]] - if len(bb[:loc[0]]) <=2: - noblock = noblock + [bb] - - blocks = newblock[1:] - - #Create a final array with blocking longitudes. Similar to IBL/GIBL, but those that pass the duration threshold - blockfreq = np.zeros((len(ibl[:,0,0]),len(ibl[0,:,0]),360)) - savecbl=[] - savemiddle = [] - saveyr=[] - numblock=0 - for i in np.arange(0,len(blocks),1): - temp = blocks[i] - numblock=numblock+1 - for j in temp: - j=int(j) - daycomp = int(dd[j]) - yearcomp = int(dy[j]) - saveyr = np.append(saveyr,dy[j]) - middlecomp = middle[j].astype(int) - mc1 = int(round(middlecomp / 2.5)) - blockfreq[yearcomp,daycomp] = blockfreq[yearcomp,daycomp] + dAfin[j] - ct = ct + 1 - - return blockfreq diff --git a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking_WeatherRegime_util.py b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking_WeatherRegime_util.py deleted file mode 100644 index d4b6be3f9..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking_WeatherRegime_util.py +++ /dev/null @@ -1,119 +0,0 @@ -import os -import netCDF4 -import numpy as np -import datetime - - -def parse_steps(): - - steps_param_fcst = os.environ.get('FCST_STEPS','') - steps_list_fcst = steps_param_fcst.split("+") - - steps_param_obs = os.environ.get('OBS_STEPS','') - steps_list_obs = steps_param_obs.split("+") - - return steps_list_fcst, steps_list_obs - - -def write_mpr_file(data_obs,data_fcst,lats_in,lons_in,time_obs,time_fcst,mname,desc,fvar,funit,flev,ovar,ounit,olev,maskname,obslev,outfile): - - dlength = len(lons_in) - bdims = data_obs.shape - - index_num = np.arange(0,dlength,1)+1 - - # Get the length of the model, FCST_VAR, FCST_LEV, OBS_VAR, OBS_LEV, VX_MASK - mname_len = str(max([5,len(mname)])+3) - desc_len = str(max([4,len(mname)])+1) - mask_len = str(max([7,len(maskname)])+3) - fvar_len = str(max([8,len(fvar)])+3) - funit_len = str(max([8,len(funit)])+3) - flev_len = str(max([8,len(flev)])+3) - ovar_len = str(max([7,len(ovar)])+3) - ounit_len = str(max([8,len(ounit)])+3) - olev_len = str(max([7,len(olev)])+3) - - format_string = '%-7s %-'+mname_len+'s %-'+desc_len+'s %-12s %-18s %-18s %-12s %-17s %-17s %-'+fvar_len+'s ' \ - '%-'+funit_len+'s %-'+flev_len+'s %-'+ovar_len+'s %-'+ounit_len+'s %-'+olev_len+'s %-10s %-'+mask_len+'s ' \ - '%-13s %-13s %-13s %-13s %-13s %-13s %-9s\n' - format_string2 = '%-7s %-'+mname_len+'s %-'+desc_len+'s %-12s %-18s %-18s %-12s %-17s %-17s %-'+fvar_len+'s ' \ - '%-'+funit_len+'s %-'+flev_len+'s %-'+ovar_len+'s %-'+ounit_len+'s %-'+olev_len+'s %-10s %-'+mask_len+'s ' \ - '%-13s %-13s %-13s %-13s %-13s %-13s %-9s %-10s %-10s %-10s %-12.4f %-12.4f %-10s %-10s %-12.4f %-12.4f ' \ - '%-10s %-10s %-10s %-10s\n' - - # Write the file - for y in range(bdims[0]): - for dd in range(bdims[1]): - if time_fcst['valid'][y][dd]: - ft_stamp = time_fcst['lead'][y][dd]+'L_'+time_fcst['valid'][y][dd][0:8]+'_' \ - +time_fcst['valid'][y][dd][9:15]+'V' - mpr_outfile_name = outfile+'_'+ft_stamp+'.stat' - with open(mpr_outfile_name, 'w') as mf: - mf.write(format_string % ('VERSION', 'MODEL', 'DESC', 'FCST_LEAD', 'FCST_VALID_BEG', 'FCST_VALID_END', - 'OBS_LEAD', 'OBS_VALID_BEG', 'OBS_VALID_END', 'FCST_VAR', 'FCST_UNITS', 'FCST_LEV', 'OBS_VAR', - 'OBS_UNITS', 'OBS_LEV', 'OBTYPE', 'VX_MASK', 'INTERP_MTHD', 'INTERP_PNTS', 'FCST_THRESH', - 'OBS_THRESH', 'COV_THRESH', 'ALPHA', 'LINE_TYPE')) - for dpt in range(dlength): - mf.write(format_string2 % ('V9.1',mname,desc,time_fcst['lead'][y][dd],time_fcst['valid'][y][dd], - time_fcst['valid'][y][dd],time_obs['lead'][y][dd],time_obs['valid'][y][dd], - time_obs['valid'][y][dd],fvar,funit,flev,ovar,ounit,olev,'ADPUPA',maskname, - 'NEAREST','1','NA','NA','NA','NA','MPR',str(dlength),str(index_num[dpt]),'NA', - lats_in[dpt],lons_in[dpt],obslev,'NA',data_fcst[y,dd,dpt],data_obs[y,dd,dpt],'NA','NA', - 'NA','NA')) - - -def read_nc_met(infiles,invar,nseasons,dseasons): - - print("Reading in Data") - - # Check to make sure that everything is not set to missing: - if all('missing' == fn for fn in infiles): - raise Exception('No input files found as given, check paths to input files') - - #Find the first non empty file name so I can get the variable sizes - locin = next(sub for sub in infiles if sub != 'missing') - indata = netCDF4.Dataset(locin) - lats = indata.variables['lat'][:] - lons = indata.variables['lon'][:] - invar_arr = indata.variables[invar][:] - indata.close() - - var_3d = np.empty([len(infiles),len(invar_arr[:,0]),len(invar_arr[0,:])]) - init_list = [] - valid_list = [] - lead_list = [] - - for i in range(0,len(infiles)): - - #Read in the data - if (infiles[i] != 'missing'): - indata = netCDF4.Dataset(infiles[i]) - new_invar = indata.variables[invar][:] - - init_time_str = indata.variables[invar].getncattr('init_time') - valid_time_str = indata.variables[invar].getncattr('valid_time') - lead_dt = datetime.datetime.strptime(valid_time_str,'%Y%m%d_%H%M%S') - datetime.datetime.strptime(init_time_str,'%Y%m%d_%H%M%S') - leadmin,leadsec = divmod(lead_dt.total_seconds(), 60) - leadhr,leadmin = divmod(leadmin,60) - lead_str = str(int(leadhr)).zfill(2)+str(int(leadmin)).zfill(2)+str(int(leadsec)).zfill(2) - indata.close() - else: - new_invar = np.empty((1,len(var_3d[0,:,0]),len(var_3d[0,0,:])),dtype=np.float) - init_time_str = '' - valid_time_str = '' - lead_str = '' - new_invar[:] = np.nan - init_list.append(init_time_str) - valid_list.append(valid_time_str) - lead_list.append(lead_str) - var_3d[i,:,:] = new_invar - - var_4d = np.reshape(var_3d,[nseasons,dseasons,len(var_3d[0,:,0]),len(var_3d[0,0,:])]) - - # Reshape time arrays and store them in a dictionary - init_list_2d = np.reshape(init_list,[nseasons,dseasons]) - valid_list_2d = np.reshape(valid_list,[nseasons,dseasons]) - lead_list_2d = np.reshape(lead_list,[nseasons,dseasons]) - time_dict = {'init':init_list_2d,'valid':valid_list_2d,'lead':lead_list_2d} - - return var_4d,lats,lons,time_dict diff --git a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking_driver.py b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking_driver.py index 91eefde4c..0f854d956 100755 --- a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking_driver.py +++ b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_Blocking/Blocking_driver.py @@ -6,10 +6,10 @@ import netCDF4 import warnings -from Blocking import BlockingCalculation +from metcalcpy.contributed.blocking_weather_regime.Blocking import BlockingCalculation +from metcalcpy.contributed.blocking_weather_regime.Blocking_WeatherRegime_util import parse_steps, write_mpr_file from metplotpy.contributed.blocking_s2s import plot_blocking as pb from metplotpy.contributed.blocking_s2s.CBL_plot import create_cbl_plot -from Blocking_WeatherRegime_util import parse_steps, write_mpr_file def main(): diff --git a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/Blocking_WeatherRegime_util.py b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/Blocking_WeatherRegime_util.py deleted file mode 120000 index 7839a8807..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/Blocking_WeatherRegime_util.py +++ /dev/null @@ -1 +0,0 @@ -../UserScript_obsERA_obsOnly_Blocking/Blocking_WeatherRegime_util.py \ No newline at end of file diff --git a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime.py b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime.py deleted file mode 100644 index 98ae69598..000000000 --- a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime.py +++ /dev/null @@ -1,257 +0,0 @@ -import os -import numpy as np -from pylab import * -from sklearn.cluster import KMeans -import scipy -import netCDF4 as nc4 -from scipy import stats,signal -from numpy import ones,vstack -from numpy.linalg import lstsq -from eofs.standard import Eof - - -class WeatherRegimeCalculation(): - """Contains the programs to perform a Weather Regime Analysis - """ - def __init__(self,label): - - self.wrnum = int(os.environ.get(label+'_WR_NUMBER',6)) - self.numi = int(os.environ.get(label+'_NUM_CLUSTERS',20)) - self.NUMPCS = int(os.environ.get(label+'_NUM_PCS',10)) - self.wr_tstep = int(os.environ.get(label+'_WR_FREQ',7)) - self.wr_outfile_type = os.environ.get(label+'_WR_OUTPUT_FILE_TYPE','text') - self.wr_outfile_dir = os.environ.get('WR_OUTPUT_FILE_DIR',os.environ['SCRIPT_OUTPUT_BASE']) - self.wr_outfile = os.environ.get(label+'_WR_OUTPUT_FILE',label+'_WeatherRegime') - - - def get_cluster_fraction(self, m, label): - return (m.labels_==label).sum()/(m.labels_.size*1.0) - - - def weights_detrend(self,lats,lons,indata): - - arr_shape = indata.shape - - ##Set up weight array - cos = lats * np.pi / 180.0 - way = np.cos(cos) - if len(lats.shape) == 1: - weightf = np.repeat(way[:,np.newaxis],len(lons),axis=1) - else: - weightf = way - - #Remove trend and seasonal cycle - atemp = np.zeros(arr_shape) - for i in np.arange(0,len(indata[0,:,0,0]),1): - atemp[:,i] = signal.detrend(atemp[:,i],axis=0) - atemp[:,i] = (indata[:,i] - np.nanmean(indata[:,i],axis=0)) * weightf - - a = atemp - atemp=0 - - #Reshape into time X space - a1 = np.reshape(a,[len(a[:,0,0,0])*len(a[0,:,0,0]),len(lons)*len(lats)]) - - return a,a1 - - - def run_elbow(self,a1): - - k=KMeans(n_clusters=self.wrnum, random_state=0) #Initilize cluster centers - - #Calculate sum of squared distances for clusters 1-15 - kind = np.arange(1,self.numi,1) - Sum_of_squared_distances = [] - K = range(1,self.numi) - for k in K: - km = KMeans(n_clusters=k) - km = km.fit(a1) - Sum_of_squared_distances.append(km.inertia_) - - # Calculate the bend of elbow - points = [(K[0],Sum_of_squared_distances[0]),(K[-1],Sum_of_squared_distances[-1])] - x_coords, y_coords = zip(*points) - A = vstack([x_coords,ones(len(x_coords))]).T - m, c = lstsq(A, y_coords,rcond=None)[0] - line = m*kind + c - curve = Sum_of_squared_distances - curve=np.array(curve)*10**-10 - line = line*10**-10 - - d=[] - for i in np.arange(0,self.numi-1,1): - p1=np.array([K[0],curve[0]]) - p2=np.array([K[-1],curve[-1]]) - p3=np.array([K[i],curve[i]]) - d=np.append(d,np.cross(p2-p1,p3-p1)/np.linalg.norm(p2-p1)) - - mi = np.where(d==d.min())[0] - print('Optimal Cluster # = '+str(mi+1)+'') - - return K,d,mi,line,curve - - - def Calc_EOF(self,eofin): - - #Remove trend and seasonal cycle - for d in np.arange(0,len(eofin[0,:,0,0]),1): - eofin[:,d] = signal.detrend(eofin[:,d],axis=0) - eofin[:,d] = eofin[:,d] - np.nanmean(eofin[:,d],axis=0) - - #Reshape into time X space - arr_shape = eofin.shape - arrdims = len(arr_shape) - eofin = np.reshape(eofin,(np.prod(arr_shape[0:arrdims-2]),arr_shape[arrdims-2]*arr_shape[arrdims-1])) - - # Use EOF solver and get PCs and EOFs - solver = Eof(eofin,center=False) - pc = solver.pcs(npcs=self.NUMPCS,pcscaling=1) - eof = solver.eofsAsCovariance(neofs=self.NUMPCS,pcscaling=1) - eof = np.reshape(eof,(self.NUMPCS,arr_shape[arrdims-2],arr_shape[arrdims-1])) - - #Get variance fractions - variance_fractions = solver.varianceFraction(neigs=self.NUMPCS) * 100 - print(variance_fractions) - - return eof, pc, self.wrnum, variance_fractions - - - def reconstruct_heights(self,eof,pc,reshape_arr): - - rssize = len(reshape_arr) - eofshape = eof.shape - eosize = len(eofshape) - - #reconstruction. If NUMPCS=nt, then a1=a0 - eofs=np.reshape(eof,(self.NUMPCS,eofshape[eosize-2]*eofshape[eosize-1])) - a1=np.matmul(pc,eofs) - - return a1 - - - def run_K_means(self,a1,timedict,arr_shape): - - arrdims = len(arr_shape) - - k=KMeans(n_clusters=self.wrnum, random_state=0) - - #fit the K-means algorithm to the data - f=k.fit(a1) - - #Obtain the cluster anomalies - y=f.cluster_centers_ - - #Obtain cluster labels for each day [Reshape to Year,day] - wr = f.labels_ - wr = np.reshape(wr,arr_shape[0:arrdims-2]) - - yf = np.reshape(y,[self.wrnum,arr_shape[arrdims-2],arr_shape[arrdims-1]]) # reshape cluster anomalies to latlon - - #Get frequency of occurrence for each cluster - perc=np.zeros(self.wrnum) - for ii in np.arange(0,self.wrnum,1): - perc[ii] = self.get_cluster_fraction(f,ii) - - #Sort % from low to high - ii = np.argsort(perc) - print(perc[ii]) - - #Reorder - perc = perc[ii] - input=yf[ii] - ii = ii[::-1] - - #Reorder from max to min and relabel - wrc = wr*1.0/1.0 - for i in np.arange(0,self.wrnum,1): - wrc[wr==ii[i]] = i+1 - - perc = perc[::-1] - input = input[::-1] - - #Save Label data [YR,DAY] - # Make some conversions first - wrc_shape = wrc.shape - len1d = wrc.size - valid_time_1d = np.reshape(timedict['valid'],len1d) - yr_1d = [] - mth_1d = [] - day_1d = [] - for vt1 in valid_time_1d: - yr_1d.append(vt1[0:4]) - mth_1d.append(vt1[4:6]) - day_1d.append(vt1[6:8]) - wrc_1d = np.reshape(wrc,len1d) - - # netcdf file - if self.wr_outfile_type=='netcdf': - wr_full_outfile = os.path.join(self.wr_outfile_dir,self.wr_outfile+'.nc') - - if os.path.isfile(wr_full_outfile): - os.remove(wr_full_outfile) - - # Create CF compliant time unit - rdate = datetime.datetime(int(yr_1d[0]),int(mth_1d[0]),int(day_1d[0]),0,0,0) - cf_diffdays = np.zeros(len(yr_1d)) - ymd_arr = np.empty(len(yr_1d)) - for dd in range(len(yr_1d)): - loopdate = datetime.datetime(int(yr_1d[dd]),int(mth_1d[dd]),int(day_1d[dd]),0,0,0) - cf_diffdays[dd] = (loopdate - rdate).days - ymd_arr[dd] = yr_1d[dd]+mth_1d[dd]+day_1d[dd] - - nc = nc4.Dataset(wr_full_outfile, 'w') - nc.createDimension('time', len(mth_1d)) - nc.Conventions = "CF-1.7" - nc.title = "Weather Regime Classification" - nc.institution = "NCAR DTCenter" - nc.source = "Weather Regime METplus use-case" - - ncti = nc.createVariable('time','d',('time')) - nc.variables['time'].long_name = "time" - nc.variables['time'].standard_name = "time" - nc.variables['time'].units = "days since "+rdate.strftime('%Y-%m-%d %H:%M:%S') - nc.variables['time'].calendar = "gregorian" - - ncdate = nc.createVariable('date','i',('time')) - nc.variables['date'].long_name = "date YYYYMMDD" - - ncnum = nc.createVariable('wrnum','i',('time'),fill_value=-9999.0) - nc.variables['wrnum'].long_name = "weather_regime_number" - - ncti[:] = cf_diffdays - ncdate[:] = ymd_arr - ncnum[:] = wrc_1d - nc.close() - - # text file - if self.wr_outfile_type=='text': - wr_full_outfile = os.path.join(self.wr_outfile_dir,self.wr_outfile+'.txt') - - if os.path.isfile(wr_full_outfile): - os.remove(wr_full_outfile) - - otdata = np.array([yr_1d, mth_1d, day_1d, wrc_1d]) - otdata = otdata.T - - with open(wr_full_outfile, 'w+') as datafile_id: - np.savetxt(datafile_id, otdata, fmt=['%6s','%3s','%4s','%6s'], header='Year Month Day WeatherRegime') - - return input, self.wrnum, perc, wrc - - - def compute_wr_freq(self, WR): - - ######## Simple Count ########## - WRfreq = np.zeros((self.wrnum,len(WR[:,0]),len(WR[0,:])-self.wr_tstep+1)) - - for yy in np.arange(0,len(WRfreq[0,:,0]),1): - d1=0;d2=self.wr_tstep - for dd in np.arange(len(WRfreq[0,0,:])): - temp = WR[yy,d1:d2] - for ww in np.arange(self.wrnum): - WRfreq[ww,yy,dd] = len(np.where(temp==ww+1)[0]) - d1=d1+1;d2=d2+1 - - dlen_plot = len(WRfreq[0,0,:]) - - return WRfreq, dlen_plot, self.wr_tstep diff --git a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime_driver.py b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime_driver.py index 790d0c1fc..e91a3082a 100755 --- a/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime_driver.py +++ b/parm/use_cases/model_applications/s2s/UserScript_obsERA_obsOnly_WeatherRegime/WeatherRegime_driver.py @@ -5,9 +5,9 @@ import netCDF4 import warnings -from WeatherRegime import WeatherRegimeCalculation +from metcalcpy.contributed.blocking_weather_regime.WeatherRegime import WeatherRegimeCalculation +from metcalcpy.contributed.blocking_weather_regime.Blocking_WeatherRegime_util import parse_steps, read_nc_met, write_mpr_file, reorder_fcst_regimes from metplotpy.contributed.weather_regime import plot_weather_regime as pwr -from Blocking_WeatherRegime_util import parse_steps, read_nc_met, write_mpr_file def main(): @@ -136,6 +136,15 @@ def main(): z500_fcst.shape) if ("KMEANS" in steps_list_obs) and ("KMEANS" in steps_list_fcst): + # Check to see if reordering the data so that the weather regime patterns match between + # the forecast and observations, is needed + #TODO: make this automated based on spatial correlations + reorder_fcst = os.environ.get('REORDER_FCST','False').lower() + fcst_order_str = os.environ['FCST_ORDER'].split(',') + fcst_order = [int(fo) for fo in fcst_order_str] + if reorder_fcst == 'true': + kmeans_fcst,perc_fcst,wrc_fcst = reorder_fcst_regimes(kmeans_fcst,perc_fcst,wrc_fcst,wrnum_fcst,fcst_order) + # Write matched pair output for weather regime classification modname = os.environ.get('MODEL_NAME','GFS') maskname = os.environ.get('MASK_NAME','FULL')