From 3b698337d685d3f8c5f800622f15513d9209b09f Mon Sep 17 00:00:00 2001 From: Gabriel Wolf Date: Mon, 17 Sep 2018 09:43:16 +0100 Subject: [PATCH] Update of get_data and calc_ref_state: running version to read data --- calc_ref_state.py | 6 +-- get_data.py | 94 ++++++++++++++++++++++++++++++++++------------ get_data.pyc | Bin 0 -> 3791 bytes 3 files changed, 72 insertions(+), 28 deletions(-) create mode 100644 get_data.pyc diff --git a/calc_ref_state.py b/calc_ref_state.py index b0a60d5..7af441f 100644 --- a/calc_ref_state.py +++ b/calc_ref_state.py @@ -26,6 +26,7 @@ variables = {'s':{},'theta':{},'grd':[],'valid':[]} # Input 1 #################################################################################### # Physical constants +# input : None # output : r_earth, rho0, grav, gbuo # r_earth : earth radis [m] # rho0 : reference density [kg m**-3] @@ -33,10 +34,7 @@ variables = {'s':{},'theta':{},'grd':[],'valid':[]} # gbuo : # PHYSICAL CONSTATS MISSING # ############################################################################################ -r_earth = 6.4e06 -rho0 = 1027.0 -grav = 9.81 -gbuo = -grav/rho0 +# constants stored in myconstants.py # step a ##################################################################################### # Read necessary data to calculate reference state diff --git a/get_data.py b/get_data.py index a9dbb04..a61b8fc 100644 --- a/get_data.py +++ b/get_data.py @@ -1,5 +1,5 @@ -def get_data_woa2009(variables): - # get_data_woa2009(variables) ############################################################ +def get_data_woa2009(variables,grid_order=('lon','lat','z','time',)): + # get_data_woa2009(variables,grid_order) ################################################# # written by : Gabriel Wolf, g.a.wolf@reading.ac.uk # adapted from get_woa2009_data.m of Remi Tailleux # last modified : 13.09.2018 @@ -16,14 +16,19 @@ def get_data_woa2009(variables): # Johnson, 2010: World Ocean Atlas 2009, Volume 2: Salinity. S. Levitus, # Ed. NOAA Atlas NESDIS 69, 184 pp. # -> PDF : https://www.nodc.noaa.gov/OC5/indpub.html#woa09 - # ####################################################################################### + # ######################################################################################## print 'Use of get_data_woa2009 in get_data.py' dir_data = '/glusterfs/inspect/users/xg911182/data/WOA2009/' grid_names = ('lon','lat','depth','time',) # import modules import numpy as np - from netCDF4 import Dataset # to read netcdf files - from mydata_classes import Grid + from netCDF4 import Dataset # to read netcdf files + from mydata_classes import Grid # class containing grid information + import myconstants as my_const # my own defined constants + from mycalc_ocean import calc_theta_from_temp + # define grid_order (used for data permutation) + find_index = lambda searchlist, elem: [[i for i, x in enumerate(searchlist) if x == e] for e in elem] + grid_order = np.squeeze(np.array(find_index(('lon','lat','z','time',),grid_order))) # Allocation data = {} grd = [] @@ -42,35 +47,76 @@ def get_data_woa2009(variables): Uz = nc_fid.variables[grid_names[2]].units Ut = nc_fid.variables[grid_names[3]].units grd = Grid(LON,LAT,Z,TIME,len(LON),len(LAT),len(Z),len(TIME), Ulon, Ulat, Uz, Ut) + # define data permutation to match input dimensions + nc_fid.variables['t_an'].dimensions + # read salinity data + if 's' in list_allkeys: + print 'Load salinity data' + fn = dir_data + 'salinity_annual_1deg.nc' + nc_fid = Dataset(fn,'r') + # read salinity data + DATA_d = nc_fid.variables['s_an'] + MASK = np.array(DATA_d)!=DATA_d._FillValue + data_perm = np.squeeze(np.asarray(find_index(DATA_d.dimensions,[grid_names[i] for i in grid_order]))) + varname = 's' + data[varname] = {} + data[varname]['val'] = np.array(DATA_d).transpose(data_perm) + data[varname]['units'] = DATA_d.units + data[varname]['fill_value'] = DATA_d._FillValue + data[varname]['standard_name'] = DATA_d.standard_name + data[varname]['valid'] = MASK.transpose(data_perm) # Read temperature data if 'temp' or 'theta' in list_allkeys: print 'Load temperature data' fn = dir_data + 'temperature_annual_1deg.nc' nc_fid = Dataset(fn,'r') # read temperature data - DATA_d = nc_fid.variables['t_an'] + DATA_d = nc_fid.variables['t_an'] + data_perm = np.squeeze(np.asarray(find_index(DATA_d.dimensions,[grid_names[i] for i in grid_order]))) + TEMP = np.array(DATA_d) + MASK = TEMP!=DATA_d._FillValue if 'temp' in list_allkeys: varname = 'temp' - data[varname] = {} - data[varname]['val'] = np.array(DATA_d) - data[varname]['units'] = DATA_d.units - data[varname]['fill_value'] = DATA_d._FillValue + data[varname] = {} + data[varname]['val'] = TEMP.transpose(data_perm) + data[varname]['units'] = DATA_d.units + data[varname]['fill_value'] = DATA_d._FillValue + data[varname]['standard_name'] = DATA_d.standard_name + data[varname]['valid'] = MASK.transpose(data_perm) if 'theta' in list_allkeys: + # compute pot. temp. for valid points + print DATA_d.dimensions + print data_perm + THETA = TEMP.transpose(data_perm) + SA = data['s']['val'] + SA = SA[data['s']['valid']] + print DATA_d.dimensions[data_perm[0]] + print DATA_d.dimensions[data_perm[1]] + print DATA_d.dimensions[data_perm[2]] + print DATA_d.dimensions[data_perm[3]] + wvar = nc_fid.variables[DATA_d.dimensions[data_perm[0]]] + xvar = nc_fid.variables[DATA_d.dimensions[data_perm[1]]] + yvar = nc_fid.variables[DATA_d.dimensions[data_perm[2]]] + zvar = nc_fid.variables[DATA_d.dimensions[data_perm[3]]] + w3d, x3d, y3d, z3d = np.meshgrid(xvar,wvar,yvar,zvar) + p = y3d*(my_const.rho0*my_const.grav/1e4) # in dbar + print 'Doesnt work automatically - p_xyz' + del w3d, x3d, y3d, z3d, wvar, xvar, yvar, zvar + pr = p*0 + print 'SA.shape: ',SA.shape + print 'THETA.shape',THETA[data[varname]['valid']].shape + print 'p.shape, ',p[data[varname]['valid']].shape + dummy = input('Press enter to continue') + MASK_mesh = data[varname]['valid'] + THETA[MASK_mesh] = calc_theta_from_temp(SA,THETA[MASK_mesh],p[MASK_mesh],pr[MASK_mesh]) + varname = 'theta' + data[varname] = {} + data[varname]['val'] = THETA + data[varname]['units'] = DATA_d.units + data[varname]['fill_value'] = DATA_d._FillValue + data[varname]['standard_name'] = 'sea_water_potential_temperature' + data[varname]['valid'] = MASK.transpose(data_perm) print 'MISSING: Calculate pot. temp. from temperature data' - # read salinity data - if 's' in list_allkeys: - print 'Load salinity data' - fn = dir_data + 'salinity_annual_1deg.nc' - nc_fid = Dataset(fn,'r') - # read salinity data - DATA_d = nc_fid.variables['s_an'] - varname = 's' - data[varname] = {} - data[varname]['val'] = np.array(DATA_d) - data[varname]['units'] = DATA_d.units - data[varname]['fill_value'] = DATA_d._FillValue - - # return data if 'grd' in list_allkeys: return grd, data diff --git a/get_data.pyc b/get_data.pyc new file mode 100644 index 0000000000000000000000000000000000000000..f1b2ee3d6179dddc5e49fe8fa73d956e4c939559 GIT binary patch literal 3791 zcmb_f&2k&Z5$;`*APN3cq)5;bwUkH+mK{Q}l1gN$TxLjFQbnRw0F_iZcDB}-0lDDr z0@_)K1Z$TU`jB_Xo1}8f6Qn9{;6olDha8iydjS#_ow5%iu(LfqJ>N`!-8};i{d=nT z&;R`8k9F$(4deYe9y5zZ6-W)!cNTpGd4y`^Q7DrkG(juM-(Hef=$|A!G4VrPw;_Cl zX9~|Y9`mSx*%m#0rfkY=QgXy z6mM;s(n&fjQ2aM5#U5v9DVw5ngwiQevy79$0!sjVnwmg>3v*N|G4{-*k@?=p?4^+_ zy^*;AoCWBz`2ht3Hzz5(GN5Bf)-4R^7>sp`13Cs{-Bn5plwxbx0~VX71pdr07c9V~ zT$t+%^DLxT5eoqCDxCpjf+Ff+t3@g`r$IP?bEf`bt0hp!bS?iVT|$&9Q7K#M^ROhV zEKzfY4u>dP7JVKz>=Wc-XqtnB;O%@fOIeATPJ!_!o8J7&z6j_FVJF!K9M-f*=?C<< z7am!n371!7Eg37R6^sLWvS(aUJVEx6&n~HJxdL(MO}P%mGF?JirNY>9bgl?>ILUtLy&(+&$Y|GymwXx_{c+v z{8(XL1{Tk*GsTD=N2QdlQi^X+v5M?vjido)N{wOqd7sN*u+Jqh*yr>|Kd;gASzJ)1 z4CK6{XS7m=O~0dO)Ka!~Nzdq|?51QbN(B>v3A3AYHcVgRw7NyvI%OM_-KKPd(ot%n zoX{@Pg1SkI3)vk~x4_5gEjnGN(+xVkO{X{L)S=TgIt4Js=V&1-lUk>zAf|a0)@lQ? zu!cM@DJCh`--iCqdBUFSlOmMVlf0#KJ zOaH8dQ2_kO8d4v^u{YNXTubMN?}I{$qqRs4n`(H~a7EP0ALA0)HMh~eHri=5oPC|R z%1b==sO5cf@7_blkIr?Q?XIc7w6Y&`Orqn4srZp;>v~e@7#*9+@&3d6_wWDalL{kL z{_y1%BUGxGLRG%Wk`XjYx08b^Uv&O|@mH1szWjDD$>=0k@Rj#8_Ep|nK(FovbvHTC z*jpoRg>It5cAYOo2E{c*5pcPg6p#bg47-l~XyzpD)eYOhO;GrseZ#ab82fBXf_zlMk0tdW1 z9012RPRIDsz5~`AcB~U?*NdW#7r6J8-rtPs;y{cW0C8eA zd|u4wu;o8TU@j)c%tQC0rGVm0^_~)0!fSx<{II);Rk^W zbw}eqg=OMJ%8Qj7!3QzL7#v{~c0g2;8D=uVQVC zqa)OXleC;Ve#H3^9Jm2XC`1^IGls?$b}5ZRB|;fVBhT zO^59^9nN?Wd%+T)%1ppTbycOyjJ`X-p_7U$FGCW#kkCjgF1tu4+dE(UvMQN94Mo_M z1E3yw#%NRJ(n{F>unVU+7KqYh-sh5R)wLH@8A;&Gn~v5ed6>L7_PP?LL)|qAllIn& zeec@EiK9q<)b1o@{t;bW=gGVjPjp+c(y?GEKbAA_VvG~dBurz}7NT-i zxGKZ0E41VcFj~Wpl#53n%O04_>J0Qk8Xp)gGEE3!z($lX!%^LB^lJ0Pm%GBti>*DW zu&)GW@A=CoBE07R57^7sND|r8>y!MI2-42h-jV#oe zw!@oqIc)FkkmP@>`iEx!{^AXozoMZlBUZsG*pt>QK4|T6Yt4qnnzTpIFIqF!ihSm* zh2H0xJq`O29