1 | # Python script to vertically nterpolate a 3D field from a netCDF file |
---|
2 | # L. Fita, LMD. Jussieu, July 2014 |
---|
3 | # |
---|
4 | ## export PATH=/u/lflmd/bin/gcc_Python-2.7.5/bin:${PATH} |
---|
5 | ## g.e. # vertical_interpolation.py -f /home/lluis/WRF/util/p_interp/wrfout_d01_1995-01-15_00:00:00 -o WRFp -i 100000.,92500.,85000.,70000.,60000.,50000.,40000.,30000.,20000.,10000. -k 'lin' -v WRFt,WRFght -d T:Time,Z:bottom_top,Y:south_north,X:west_east -D T:Times,Z:ZNU,Y:XLAT,X:XLONG |
---|
6 | ## g.e. # vertical_interpolation.py -f /media/data1/etudes/WRF_LMDZ/WaquaL/WRF_LMDZ/NPv31/wrfout/wrfout_d01_1980-03-01_00:00:00 -o WRFp -i 100000.,97500.,95000.,92500.,90000.,85000.,80000.,75000.,70000.,65000.,60000.,55000.,50000.,45000.,40000.,35000.,30000.,25000.,20000.,15000.,10000. -k 'lin' -v WRFt,U,V,WRFrh,WRFght -d T:Time,Z:bottom_top,Y:south_north,X:west_east -D T:Times,Z:ZNU,Y:XLAT,X:XLONG |
---|
7 | |
---|
8 | import numpy as np |
---|
9 | from netCDF4 import Dataset as NetCDFFile |
---|
10 | import os |
---|
11 | import re |
---|
12 | import nc_var_tools as ncvar |
---|
13 | from optparse import OptionParser |
---|
14 | |
---|
15 | def lonlat_creation(dimpn,nco,nwnc): |
---|
16 | """ Function to create longitude/latitude variables according to a dictionary |
---|
17 | dimpn: dictionary of varibale names for each dimension: 'T', 'Z', 'Y', 'x' |
---|
18 | nco: original netCDF object |
---|
19 | nwnc: new netCDF object |
---|
20 | """ |
---|
21 | fname = 'lonlat_creation' |
---|
22 | |
---|
23 | print 'Lluis:',dimpn['X'] |
---|
24 | |
---|
25 | if dimpn.has_key('X') and dimpn.has_key('Y'): |
---|
26 | lonobj = nco.variables[dimpn['X']] |
---|
27 | if len(lonobj.shape) == 3: |
---|
28 | newvar = nwnc.createVariable('lon', 'f8', ('y', 'x')) |
---|
29 | neattr = ncvar.basicvardef(newvar, 'lon', 'longitude', 'degrees East') |
---|
30 | newvar[:] = lonobj[0,:,:] |
---|
31 | |
---|
32 | newvar = nwnc.createVariable('lat', 'f8', ('y', 'x')) |
---|
33 | neattr = ncvar.basicvardef(newvar, 'lat', 'latitude', 'degrees Nord') |
---|
34 | newvar[:] = nco.variables[dimvns[1]][0,:,:] |
---|
35 | elif len(lonobj.shape) == 2: |
---|
36 | newvar = nwnc.createVariable('lon', 'f8', ('y', 'x')) |
---|
37 | neattr = ncvar.basicvardef(newvar, 'lon', 'longitude', 'degrees East') |
---|
38 | newvar[:] = lonobj[:] |
---|
39 | |
---|
40 | newvar = nwnc.createVariable('lat', 'f8', ('y', 'x')) |
---|
41 | neattr = ncvar.basicvardef(newvar, 'lat', 'latitude', 'degrees Nord') |
---|
42 | newvar[:] = nco.variables[dimvns[1]][:] |
---|
43 | else: |
---|
44 | print errormsg |
---|
45 | print ' ' + main + ': shape of longitude: ',lonobj.shape,' not ready !!' |
---|
46 | quit(-1) |
---|
47 | |
---|
48 | elif dimpn.has_key('X') and not dimpn.has_key('Y'): |
---|
49 | print warnmsg |
---|
50 | print ' ' + main + ": variable pressure with 'X', but not 'Y'!!" |
---|
51 | lonobj = nco.variables[dimpn['X']] |
---|
52 | if len(lonobj.shape) == 3: |
---|
53 | newvar = nwnc.createVariable('lon', 'f8', ('x')) |
---|
54 | neattr = ncvar.basicvardef(newvar, 'lon', 'longitude', 'degrees East') |
---|
55 | newvar[:] = lonobj[0,0,:] |
---|
56 | |
---|
57 | elif len(lonobj.shape) == 2: |
---|
58 | newvar = nwnc.createVariable('lon', 'f8', ('x')) |
---|
59 | neattr = ncvar.basicvardef(newvar, 'lon', 'longitude', 'degrees East') |
---|
60 | newvar[:] = lonobj[0,:] |
---|
61 | |
---|
62 | elif len(lonobj.shape) == 1: |
---|
63 | newvar = nwnc.createVariable('lon', 'f8', ('x')) |
---|
64 | neattr = ncvar.basicvardef(newvar, 'lon', 'longitude', 'degrees East') |
---|
65 | newvar[:] = lonobj[:] |
---|
66 | else: |
---|
67 | print errormsg |
---|
68 | print ' ' + main + ': shape of longitude: ',lonobj.shape,' not ready !!' |
---|
69 | quit(-1) |
---|
70 | |
---|
71 | elif not dimpn.has_key('X') and dimpn.has_key('Y'): |
---|
72 | print warnmsg |
---|
73 | print ' ' + main + ": variable pressure with 'Y', but not 'X'!!" |
---|
74 | latobj = nco.variables[dimpn['Y']] |
---|
75 | if len(latobj.shape) == 3: |
---|
76 | newvar = nwnc.createVariable('lat', 'f8', ('y')) |
---|
77 | neattr = ncvar.basicvardef(newvar, 'lat', 'latitude', 'degrees North') |
---|
78 | newvar[:] = latobj[0,0,:] |
---|
79 | |
---|
80 | elif len(latobj.shape) == 2: |
---|
81 | newvar = nwnc.createVariable('lat', 'f8', ('x')) |
---|
82 | neattr = ncvar.basicvardef(newvar, 'lat', 'latitude', 'degrees North') |
---|
83 | newvar[:] = latobj[0,:] |
---|
84 | |
---|
85 | elif len(latobj.shape) == 1: |
---|
86 | newvar = nwnc.createVariable('lat', 'f8', ('x')) |
---|
87 | neattr = ncvar.basicvardef(newvar, 'lat', 'latitude', 'degrees North') |
---|
88 | newvar[:] = latobj[:] |
---|
89 | else: |
---|
90 | print errormsg |
---|
91 | print ' ' + main + ': shape of latitude: ',latobj.shape,' not ready !!' |
---|
92 | quit(-1) |
---|
93 | else: |
---|
94 | print warnmsg |
---|
95 | print ' ' + main + ": variable pressure without 'X' and 'Y'!!" |
---|
96 | |
---|
97 | return |
---|
98 | |
---|
99 | ####### ###### ##### #### ### ## # |
---|
100 | |
---|
101 | main='vertical_interpolation.py' |
---|
102 | errormsg='ERROR -- error -- ERROR -- error' |
---|
103 | warnmsg='WARNING -- warning -- WARNING -- warning' |
---|
104 | |
---|
105 | parser = OptionParser() |
---|
106 | parser.add_option("-o", "--original_ZValues", dest="ozvals", |
---|
107 | help="variable name with the original z-values ('WRFp', WRF derived pressure values)", metavar="VALUE") |
---|
108 | parser.add_option("-i", "--interpolate_ZValues", dest="izvals", |
---|
109 | help="comma-separated list of values to interpolate ([zi1],[zi2],[...[ziN]])", metavar="VALUES") |
---|
110 | parser.add_option("-f", "--file", dest="file", |
---|
111 | help="netCDF file to use", metavar="FILE") |
---|
112 | parser.add_option("-k", "--interpolation_kind", dest="kfig", |
---|
113 | help="kind of vertical inerpolation ('lin', linear)", metavar="LABEL") |
---|
114 | parser.add_option("-v", "--variables", dest="vars", |
---|
115 | help="comma-separated list of name of the variables to inerpolate ('all', all variables)", |
---|
116 | metavar="LABEL") |
---|
117 | parser.add_option("-d", "--dimensions", dest="dims", |
---|
118 | help="comma-separated list of dimensions name (where applicable, T:[dimt],Z:[dimz],Y:[dimy],X:[dimx])", |
---|
119 | metavar="LABEL") |
---|
120 | parser.add_option("-D", "--dimension_vnames", dest="dimvns", |
---|
121 | help="comma-separated list of variables with the dimensions values (T:[dimt],Y:[dimy],X:[dimx])", metavar="LABEL") |
---|
122 | |
---|
123 | (opts, args) = parser.parse_args() |
---|
124 | |
---|
125 | ####### ###### ##### #### ### ## # |
---|
126 | |
---|
127 | # Variables wich do not exist in the file, but they will be computed |
---|
128 | NOfileVars = ['WRFght', 'WRFrh', 'WRFt'] |
---|
129 | |
---|
130 | ####### ####### |
---|
131 | ## MAIN |
---|
132 | ####### |
---|
133 | |
---|
134 | ofile = 'vertical_interpolation_' + opts.ozvals +'.nc' |
---|
135 | |
---|
136 | if not os.path.isfile(opts.file): |
---|
137 | print errormsg |
---|
138 | print ' ' + main + ' file: "' + opts.file + '" does not exist !!' |
---|
139 | quit(-1) |
---|
140 | |
---|
141 | ncobj = NetCDFFile(opts.file, 'r') |
---|
142 | |
---|
143 | if opts.ozvals == 'WRFp': |
---|
144 | ncdims = ncobj.variables['P'].dimensions |
---|
145 | else: |
---|
146 | ncdims = ncobj.dimensions |
---|
147 | |
---|
148 | |
---|
149 | dimns = opts.dims.split(',') |
---|
150 | |
---|
151 | # Variables with values of the dimensions |
---|
152 | dimvns = [] |
---|
153 | d4dimvns = {} |
---|
154 | |
---|
155 | for i in range(4): |
---|
156 | dimvid = opts.dimvns.split(',')[i].split(':')[0] |
---|
157 | dimvv = opts.dimvns.split(',')[i].split(':')[1] |
---|
158 | |
---|
159 | d4dimvns[dimvid] = dimvv |
---|
160 | if dimvid != 'Z': dimvns.append(dimvv) |
---|
161 | |
---|
162 | dimpresv = {} |
---|
163 | idim = 0 |
---|
164 | dimpresn = {} |
---|
165 | dimpresnv = {} |
---|
166 | d4dims = {} |
---|
167 | for dim in dimns: |
---|
168 | dimid = dim.split(':')[0] |
---|
169 | dimn = dim.split(':')[1] |
---|
170 | if not ncvar.searchInlist(ncdims, dimn): |
---|
171 | print warnmsg |
---|
172 | print ' ' + main + ' in file: "' + opts.file + '" dimension "' + dimn + \ |
---|
173 | '" does not exist !!' |
---|
174 | # quit(-1) |
---|
175 | else: |
---|
176 | dimpresv[dimid] = len(ncobj.dimensions[dimn]) |
---|
177 | dimpresn[dimid] = dimn |
---|
178 | dimpresnv[dimn] = len(ncobj.dimensions[dimn]) |
---|
179 | |
---|
180 | d4dims[dimid] = dimn |
---|
181 | idim = idim + 1 |
---|
182 | |
---|
183 | print 'Generic 4d dimensions:',d4dims |
---|
184 | |
---|
185 | dimvl = [] |
---|
186 | for dimid in ['T', 'Z', 'Y', 'X']: |
---|
187 | if dimpresv.has_key(dimid): |
---|
188 | dimvl.append(dimpresv[dimid]) |
---|
189 | |
---|
190 | dimv = np.array(dimvl) |
---|
191 | print 'dimensions pressure variable: ',dimpresnv |
---|
192 | |
---|
193 | # Interpolation values |
---|
194 | intervals = opts.izvals.split(',') |
---|
195 | intvals = ncvar.list_toVec(intervals, 'npfloat') |
---|
196 | Nintvals = len(intvals) |
---|
197 | |
---|
198 | # Obtaining original vertical coordinate |
---|
199 | if opts.ozvals == 'WRFp': |
---|
200 | print ' ' + main + ': Retrieving original pressure values from WRF ' + \ |
---|
201 | 'files as P + PB' |
---|
202 | zorigvals = np.zeros((dimv), dtype = np.float) |
---|
203 | zorigvals = ncobj.variables['P'][:] + ncobj.variables['PB'][:] |
---|
204 | zorigname = 'pressure' |
---|
205 | zorigsn = 'pressure' |
---|
206 | zorigln = 'pressure of the air' |
---|
207 | zorigu = 'hPa' |
---|
208 | else: |
---|
209 | zorigvals = ncobj.variables[opts.ozvals][:] |
---|
210 | zorigname = opts.ozvals |
---|
211 | varattrs = ncobj.variables[opts.ozvals].ncattrs() |
---|
212 | if ncvar.searchInlist(varattrs,'standard_name'): |
---|
213 | zorigsn = ncobj.variables[var].getncattr('standard_name') |
---|
214 | else: |
---|
215 | zorigsn = var |
---|
216 | |
---|
217 | if ncvar.searchInlist(varattrs,'long_name'): |
---|
218 | zorigln = ncobj.variables[var].getncattr('long_name') |
---|
219 | else: |
---|
220 | zorigln = ncobj.variables[var].getncattr('description') |
---|
221 | |
---|
222 | zorigu = ncobj.variables[var].getncattr('units') |
---|
223 | |
---|
224 | # Which is the sense of the original vertical coordinate? |
---|
225 | Ndimszorig = len(zorigvals.shape) |
---|
226 | |
---|
227 | if Ndimszorig == 4: |
---|
228 | Lzorig = zorigvals.shape[1] |
---|
229 | inczorig = zorigvals[0,Lzorig-1,0,0] - zorigvals[0,0,0,0] |
---|
230 | elif Ndimszorig == 3: |
---|
231 | Lzorig = zorigvals.shape[0] |
---|
232 | inczorig = zorigvals[Lzorig-1,0,0] - zorigvals[0,0,0] |
---|
233 | elif Ndimszorig == 2: |
---|
234 | # Assuming Time,z |
---|
235 | Lzorig = zorigvals.shape[1] |
---|
236 | inczorig = zorigvals[0,Lzorig-1] - zorigvals[0,0] |
---|
237 | else: |
---|
238 | print errormsg |
---|
239 | print ' ' + main + ': vertical coordinate shape:',zorigvals.shape,'not ready!!!' |
---|
240 | quit(-1) |
---|
241 | |
---|
242 | print ' ' + main + ': vertical coordinate to interpolate:',dimpresn['Z'] |
---|
243 | # Which are the desired variables |
---|
244 | if opts.vars == 'all': |
---|
245 | varns0 = ncobj.variables |
---|
246 | varns = [] |
---|
247 | # only getting that variables which contain dimns[1] |
---|
248 | for varn in varns0: |
---|
249 | vobj = ncobj.variables[varn] |
---|
250 | if ncvar.searchInlist(vobj.dimensions, dimpresn['Z']): |
---|
251 | varns.append(varn) |
---|
252 | print ' ' + main + ': Interpolation of all variables !!!' |
---|
253 | print ' ', varns |
---|
254 | else: |
---|
255 | varns = opts.vars.split(',') |
---|
256 | |
---|
257 | # Creation of output file |
---|
258 | ## |
---|
259 | newnc = NetCDFFile(ofile, 'w') |
---|
260 | |
---|
261 | # Creation of dimensions |
---|
262 | if dimpresn.has_key('X'): newdim = newnc.createDimension('x', dimpresv['X']) |
---|
263 | if dimpresn.has_key('Y'): newdim = newnc.createDimension('y', dimpresv['Y']) |
---|
264 | if dimpresn.has_key('Z'): newdim = newnc.createDimension('z', Nintvals) |
---|
265 | if dimpresn.has_key('T'): |
---|
266 | if not newnc.dimensions.has_key('time'): |
---|
267 | timeobj = ncobj.variables[d4dimvns['T']] |
---|
268 | if opts.ozvals == 'WRFp': |
---|
269 | newdim = newnc.createDimension('time', None) |
---|
270 | newvar = newnc.createVariable('time', 'f8', ('time')) |
---|
271 | |
---|
272 | timevals = np.zeros((dimv[0]), dtype=np.float) |
---|
273 | |
---|
274 | if d4dimvns['T'] == 'Times': |
---|
275 | timeu = 'hours since 1949-12-01 00:00:00' |
---|
276 | for it in range(dimv[0]): |
---|
277 | gdate = ncvar.datetimeStr_conversion(timeobj[it,:],'WRFdatetime',\ |
---|
278 | 'matYmdHMS') |
---|
279 | timevals[it] = ncvar.realdatetime1_CFcompilant(gdate, \ |
---|
280 | '19491201000000', 'hours') |
---|
281 | else: |
---|
282 | timevals = timeobj[:] |
---|
283 | timeu = timeobj.getncattr('units') |
---|
284 | |
---|
285 | newattr = ncvar.basicvardef(newvar, 'time', 'time', timeu) |
---|
286 | newvar[:] = timevals |
---|
287 | |
---|
288 | # Varibales dimension |
---|
289 | #ofunc = lonlat_creation(dimpresn,ncobj,newnc) |
---|
290 | ofunc = lonlat_creation(d4dimvns,ncobj,newnc) |
---|
291 | |
---|
292 | newvar = newnc.createVariable(zorigname, 'f8', ('z')) |
---|
293 | neattr = ncvar.basicvardef(newvar, zorigsn, zorigln, zorigu) |
---|
294 | newvar[:] = intvals |
---|
295 | |
---|
296 | # Looping along the variables |
---|
297 | for var in varns: |
---|
298 | print 'variable:',var |
---|
299 | |
---|
300 | if not ncvar.searchInlist(ncobj.variables,var) and not \ |
---|
301 | ncvar.searchInlist(NOfileVars,var): |
---|
302 | print errormsg |
---|
303 | print ' ' + main + ' in file: "' + opts.file + '" variable "' + var + \ |
---|
304 | '" does not exist !!' |
---|
305 | quit(-1) |
---|
306 | |
---|
307 | if ncvar.searchInlist(NOfileVars,var): |
---|
308 | if var == 'WRFght': |
---|
309 | print ' ' + main + ': computing geopotential height from WRF as ' + \ |
---|
310 | ' PH + PHB ...' |
---|
311 | varvals = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:] |
---|
312 | varsn = 'zg' |
---|
313 | varln = 'geopotential height' |
---|
314 | varu = 'gpm' |
---|
315 | varinterp = np.zeros((dimv[0], Nintvals, dimv[2], dimv[3]), dtype=np.float) |
---|
316 | newvar = newnc.createVariable(var, 'f4', ('time','z','y','x')) |
---|
317 | elif var == 'WRFrh': |
---|
318 | print ' ' + main + ': computing relative humidity from WRF as ' + \ |
---|
319 | " Tetens' equation (T,P) ..." |
---|
320 | p0=100000. |
---|
321 | p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] |
---|
322 | |
---|
323 | tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) |
---|
324 | qv = ncobj.variables['QVAPOR'][:] |
---|
325 | |
---|
326 | data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65)) |
---|
327 | data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1) |
---|
328 | varvals = qv/data2 |
---|
329 | |
---|
330 | varsn = 'rh' |
---|
331 | varln = 'relative humidity of the air' |
---|
332 | varu = '%' |
---|
333 | varinterp = np.zeros((dimv[0], Nintvals, dimv[2], dimv[3]), dtype=np.float) |
---|
334 | newvar = newnc.createVariable(var, 'f4', ('time','z','y','x')) |
---|
335 | elif var == 'WRFt': |
---|
336 | print ' ' + main + ': computing temperature from WRF as ' + \ |
---|
337 | ' inv_potT(T + 300) ...' |
---|
338 | p0=100000. |
---|
339 | p=ncobj.variables['P'][:] + ncobj.variables['PB'][:] |
---|
340 | |
---|
341 | varvals = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.) |
---|
342 | varsn = 'ta' |
---|
343 | varln = 'temperature of the air' |
---|
344 | varu = 'K' |
---|
345 | varinterp = np.zeros((dimv[0], Nintvals, dimv[2], dimv[3]), dtype=np.float) |
---|
346 | newvar = newnc.createVariable(var, 'f4', ('time','z','y','x')) |
---|
347 | else: |
---|
348 | print errormsg |
---|
349 | print ' ' + fname + ': variable "' + var + '" not ready !!!!' |
---|
350 | quit(-1) |
---|
351 | else: |
---|
352 | varobj = ncobj.variables[var] |
---|
353 | vardims = varobj.dimensions |
---|
354 | varattrs = varobj.ncattrs() |
---|
355 | if ncvar.searchInlist(varattrs,'standard_name'): |
---|
356 | varsn = ncobj.variables[var].getncattr('standard_name') |
---|
357 | else: |
---|
358 | varsn = var |
---|
359 | |
---|
360 | if ncvar.searchInlist(varattrs,'long_name'): |
---|
361 | varln = ncobj.variables[var].getncattr('long_name') |
---|
362 | else: |
---|
363 | varln = ncobj.variables[var].getncattr('description') |
---|
364 | |
---|
365 | varu = ncobj.variables[var].getncattr('units') |
---|
366 | |
---|
367 | # Getting variable values: |
---|
368 | if len(varobj.shape) == 4: |
---|
369 | varinterp = np.zeros((dimv[0], Nintvals, dimv[2], dimv[3]), dtype=np.float) |
---|
370 | newvar = newnc.createVariable(var, 'f4', ('time','z','y','x')) |
---|
371 | varvals = varobj[:] |
---|
372 | elif len(varobj.shape) <= 3 and len(varobj.shape) >= 1: |
---|
373 | varpdimvs = [] |
---|
374 | varpdimns = [] |
---|
375 | varpslice = [] |
---|
376 | for vdim in vardims: |
---|
377 | vardim = len(ncobj.dimensions[vdim]) |
---|
378 | if ncvar.searchInlist(ncdims,vdim): |
---|
379 | if vdim == dimpresn['Z']: |
---|
380 | varpdimvs.append(Nintvals) |
---|
381 | varpdimns.append('z') |
---|
382 | varpslice.append(slice(0,vardim)) |
---|
383 | else: |
---|
384 | varpdimvs.append(vardim) |
---|
385 | varpslice.append(slice(0,vardim)) |
---|
386 | if dimpresn.has_key('X') and vdim == dimpresn['X']: |
---|
387 | varpdimns.append('x') |
---|
388 | elif dimpresn.has_key('Y') and vdim == dimpresn['Y']: |
---|
389 | varpdimns.append('y') |
---|
390 | elif dimpresn.has_key('T') and vdim == dimpresn['T']: |
---|
391 | varpdimns.append('time') |
---|
392 | else: |
---|
393 | print errormsg |
---|
394 | print ' ' + main + ": dimension variable '" + vdim + \ |
---|
395 | "' is in pressure but it is not found?" |
---|
396 | print ' pressure dimensions:', ncdims |
---|
397 | quit(-1) |
---|
398 | else: |
---|
399 | # Dimension of the variable is not in the pressure variable |
---|
400 | varpdimvs.append(vardim) |
---|
401 | varpdimns.append(vdim) |
---|
402 | varpslice.append(slice(0,vardim)) |
---|
403 | if not newnc.dimensions.has_key(vdim) and not \ |
---|
404 | ncvar.searchInlist(dimpresn,vdim): |
---|
405 | print ' ' + main + ": dimension '" + vdim + "' not in the " + \ |
---|
406 | 'pressure variable adding it!' |
---|
407 | if ncobj.dimensions[vdim].isunlimited(): |
---|
408 | newnc.createDimension(vdim, None) |
---|
409 | else: |
---|
410 | newnc.createDimension(vdim, vardim) |
---|
411 | |
---|
412 | varinterp = np.zeros((varpdimvs), dtype=np.float) |
---|
413 | newvar = newnc.createVariable(var, 'f4', tuple(varpdimns)) |
---|
414 | |
---|
415 | varvals = varobj[tuple(varpslice)] |
---|
416 | else: |
---|
417 | print errormsg |
---|
418 | print ' ' + main + ': variable shape "', varobj.shape, '" not ready !!!!' |
---|
419 | quit(-1) |
---|
420 | |
---|
421 | # print 'variable:',var,'shape:',varvals.shape,'len:',len(varvals.shape) |
---|
422 | if len(varvals.shape) == 3 and len(zorigvals.shape) == 3: |
---|
423 | if (varvals.shape[1] + varvals.shape[2]) != (dimv[1] + dimv[2]): |
---|
424 | print warnmsg |
---|
425 | print ' ' + main + ': variable=', varvals.shape[1:3], \ |
---|
426 | 'with different shape!',dimv[1],dimv[2] |
---|
427 | if (varvals.shape[1] + varvals.shape[2]) - (dimv[1] + dimv[2]) == 1: |
---|
428 | print ' Assuming staggered variable' |
---|
429 | varvals0 = np.zeros((Nintvals, dimv[1], dimv[2]), dtype=np.float) |
---|
430 | |
---|
431 | if (varvals.shape[1] > dimv[1]): |
---|
432 | varvals0 = (varvals[0:dimv[1],:] + varvals[1:dimv[1]+1,:])/2. |
---|
433 | else: |
---|
434 | varvals0 = (varvals[:,0:dimv[2]] + varvals[:,1:dimv[2]+1])/2. |
---|
435 | varinterp = ncvar.interpolate_3D(zorigvals, varvals0, intvals, 'lin') |
---|
436 | else: |
---|
437 | varinterp = ncvar.interpolate_3D(zorigvals, varvals, intvals, 'lin') |
---|
438 | |
---|
439 | elif len(varvals.shape) == 4 and len(zorigvals.shape) == 4: |
---|
440 | if (varvals.shape[2] + varvals.shape[3]) != (dimv[2] + dimv[3]): |
---|
441 | print warnmsg |
---|
442 | print ' ' + main + ': variable=', varvals.shape[2:4], \ |
---|
443 | 'with different shape!',dimv[2],dimv[3] |
---|
444 | if (varvals.shape[2] + varvals.shape[3]) - (dimv[2] + dimv[3]) == 1: |
---|
445 | print ' Assuming staggered variable' |
---|
446 | varvals0 = np.zeros((dimv[0],dimv[1],dimv[2],dimv[3]), \ |
---|
447 | dtype=np.float) |
---|
448 | if (varvals.shape[2] > dimv[2]): |
---|
449 | varvals0=(varvals[:,:,0:dimv[2],:]+varvals[:,:,1:dimv[2]+1,:])/2. |
---|
450 | else: |
---|
451 | varvals0=(varvals[:,:,:,0:dimv[3]]+varvals[:,:,:,1:dimv[3]+1])/2. |
---|
452 | |
---|
453 | for it in range(dimv[0]): |
---|
454 | ncvar.percendone(it, dimv[0], 5, 'interpolated') |
---|
455 | |
---|
456 | varinterp[it,:,:,:] = ncvar.interpolate_3D(zorigvals[it,:,:,:], \ |
---|
457 | varvals0[it,:,:,:], intvals, 'lin') |
---|
458 | else: |
---|
459 | for it in range(dimv[0]): |
---|
460 | ncvar.percendone(it, dimv[0], 5, 'interpolated') |
---|
461 | |
---|
462 | varinterp[it,:,:,:] = ncvar.interpolate_3D(zorigvals[it,:,:,:], \ |
---|
463 | varvals[it,:,:,:], intvals, 'lin') |
---|
464 | # print varinterp[it,:,:,:] |
---|
465 | # quit() |
---|
466 | elif len(varvals.shape) == 2 and len(zorigvals.shape) == 2: |
---|
467 | zdimid = varpdimns.index('z') |
---|
468 | varinterp = ncvar.interpolate_2D(zorigvals, varvals, zdimid, intvals, 'lin') |
---|
469 | |
---|
470 | elif len(varvals.shape) == 3 and len(zorigvals.shape) == 2: |
---|
471 | print 'Lluis here variable with an extra dimension in comparison to pressure' |
---|
472 | zdimid = varpdimns.index('z') |
---|
473 | presdn = [] |
---|
474 | for idi in range(2): |
---|
475 | presdn.append(dimpresn.keys()[idi].lower()) |
---|
476 | for dimid in range(3): |
---|
477 | if not ncvar.searchInlist(presdn,varpdimns): loopd = dimid |
---|
478 | # Looping along the extra dimension |
---|
479 | if loopd == 0: |
---|
480 | for il in range(varvals.shape[loopd]): |
---|
481 | varinterp[il,:,:] = ncvar.interpolate_2D(zorigvals, varvals[il,:,:], \ |
---|
482 | zdimid, intvals, 'lin') |
---|
483 | elif loopd == 1: |
---|
484 | for il in range(varvals.shape[loopd]): |
---|
485 | varinterp[:,il,:] = ncvar.interpolate_2D(zorigvals, varvals[:,il,:], \ |
---|
486 | zdimid, intvals, 'lin') |
---|
487 | elif loopd == 2: |
---|
488 | for il in range(varvals.shape[loopd]): |
---|
489 | varinterp[:,:,il] = ncvar.interpolate_2D(zorigvals, varvals[:,:,il], \ |
---|
490 | zdimid, intvals, 'lin') |
---|
491 | else: |
---|
492 | # print errormsg |
---|
493 | print warnmsg |
---|
494 | print ' ' + main + ': dimension of values:',varvals.shape, \ |
---|
495 | ' not ready to interpolate using:',zorigvals.shape,'!!!' |
---|
496 | print ' skipping variable' |
---|
497 | # quit(-1) |
---|
498 | |
---|
499 | print 'v_interp Lluis dims newvar:',newvar.dimensions |
---|
500 | print 'v_interp Lluis shapes: newvar',newvar.shape,'varinterp:',varinterp.shape |
---|
501 | newvar[:] = varinterp |
---|
502 | ncvar.basicvardef(newvar, varsn, varln, varu) |
---|
503 | |
---|
504 | newnc.sync() |
---|
505 | |
---|
506 | # Global attributes |
---|
507 | ## |
---|
508 | |
---|
509 | atvar = ncvar.set_attribute(newnc, 'program', 'vertical_inerpolation.py') |
---|
510 | atvar = ncvar.set_attribute(newnc, 'version', '1.0') |
---|
511 | atvar = ncvar.set_attribute(newnc, 'author', 'Fita Borrell, Lluis') |
---|
512 | atvar = ncvar.set_attribute(newnc, 'institution', 'Laboratoire Meteorologie ' + \ |
---|
513 | 'Dynamique') |
---|
514 | atvar = ncvar.set_attribute(newnc, 'university', 'Universite Pierre et Marie ' + \ |
---|
515 | 'Curie -- Jussieu') |
---|
516 | atvar = ncvar.set_attribute(newnc, 'centre', 'Centre national de la recherche ' + \ |
---|
517 | 'scientifique') |
---|
518 | atvar = ncvar.set_attribute(newnc, 'city', 'Paris') |
---|
519 | atvar = ncvar.set_attribute(newnc, 'original_file', opts.file) |
---|
520 | |
---|
521 | gorigattrs = ncobj.ncattrs() |
---|
522 | |
---|
523 | for attr in gorigattrs: |
---|
524 | attrv = ncobj.getncattr(attr) |
---|
525 | atvar = ncvar.set_attribute(newnc, attr, attrv) |
---|
526 | |
---|
527 | ncobj.close() |
---|
528 | newnc.close() |
---|
529 | |
---|
530 | print main + ': successfull writting of verticaly interpolated file "'+ofile+'" !!!' |
---|