source: lmdz_wrf/trunk/tools/get_stations.py @ 1701

Last change on this file since 1701 was 1701, checked in by lfita, 7 years ago

Adding more general generation of 2D X,Y variables
File size: 24.3 KB
Line 
1## Getting statistics from a series of surface stations and soundings from any given netCDF file
2# L. Fita, CIMA, November 2017
3#
4import numpy as np
5#import matplotlib as mpl
6#mpl.use('Agg')
7#from matplotlib.pylab import *
8#import matplotlib.pyplot as plt
9#from mpl_toolkits.basemap import Basemap
10import os
11from netCDF4 import Dataset as NetCDFFile
12import nc_var_tools as ncvar
13import generic_tools as gen
14#import drawing_tools as drw
15import diag_tools as diag
16
17# Getting configuration
18from config_get_stations import *
19
20####### ###### ##### #### ### ## #
21
22# Gneric variables prepared to be computed
23Cvars = ['C_td', 'C_wd', 'C_ws']
24
25# WRF specific diagnostics
26WRFvars = ['WRFp', 'WRFta', 'WRFtd', 'WRFtime', 'WRFua', 'WRFuas', 'WRFva', 'WRFvas',\
27  'WRFzg']
28
29# Variables not to check their existence inside file
30NONcheckingvars = Cvars + WRFvars + ['None']
31
32def creation_sfcstation_file(filen, lv, Lv, hv, lab, tunits, sfcvars, dimt, ifilen):
33    """ Function to create the structure of the surface station file
34      filen: name of the file
35      lv: value of the longitude of the station
36      Lv: value of the latitude of the station
37      hv: value of the height of the station
38      lab: label of the station
39      tunits: uints of time
40      sfcvars: variables to be included
41      dimt: quantity of time-steps
42      ifilen: name of the file from which data is retrieved
43    """
44    fname = 'creation_sfcstation_file'
45
46    Lstring = 256
47
48    onewnc = NetCDFFile(filen, 'w')
49    # Dimensions
50    newdim = onewnc.createDimension('time', None)
51    newdim = onewnc.createDimension('Lstring', Lstring)
52
53    # Variable-dimensions
54    newvar = onewnc.createVariable('time', 'f8', ('time'))
55    ncvar.basicvardef(newvar, 'time', 'Time', tunits)
56    ncvar.set_attribute(newvar, 'calendar', 'standard')
57    newvar.setncattr('axis', 'T')
58    newvar.setncattr('_CoordinateAxisType', 'Time')
59    onewnc.sync()
60
61    # station Variables
62    Llab = len(lab)
63    newvar = onewnc.createVariable('station', 'c', ('Lstring'))
64    newvar[0:Llab] = lab[:]
65    ncvar.basicvardef(newvar, 'station', 'Name of the station', '-')
66
67    newvar = onewnc.createVariable('lon', 'f8')
68    newvar[:] = lv
69    ncvar.basicvardef(newvar, 'lon', 'Longitude', 'degrees_west')
70    newvar.setncattr('axis', 'X')
71    newvar.setncattr('_CoordinateAxisType', 'Lon')
72   
73    newvar = onewnc.createVariable('lat', 'f8')
74    newvar[:] = Lv
75    ncvar.basicvardef(newvar, 'lat', 'Latitude', 'degrees_north')
76    newvar.setncattr('axis', 'Y')
77    newvar.setncattr('_CoordinateAxisType', 'Lat')
78   
79    newvar = onewnc.createVariable('height', 'f8')
80    newvar[:] = hv
81    ncvar.basicvardef(newvar, 'height', 'Height', 'm')
82    newvar.setncattr('axis', 'Z')
83    newvar.setncattr('_CoordinateAxisType', 'Height')
84    onewnc.sync()
85
86    newvar = onewnc.createVariable('flon', 'f8')
87    newvar[:] = lv
88    ncvar.basicvardef(newvar, 'file_lon', 'Longitude closest grid-point from file',  \
89      'degrees_west')
90    newvar.setncattr('axis', 'X')
91    newvar.setncattr('_CoordinateAxisType', 'Lon')
92   
93    newvar = onewnc.createVariable('flat', 'f8')
94    newvar[:] = Lv
95    ncvar.basicvardef(newvar, 'file_lat', 'Latitude closest grid-point from file',   \
96      'degrees_north')
97    newvar.setncattr('axis', 'Y')
98    newvar.setncattr('_CoordinateAxisType', 'Lat')
99   
100    newvar = onewnc.createVariable('fheight', 'f8')
101    newvar[:] = hv
102    ncvar.basicvardef(newvar, 'file_height', 'Height closest grid-point from file',  \
103      'm')
104    newvar.setncattr('axis', 'Z')
105    newvar.setncattr('_CoordinateAxisType', 'Height')
106
107    newvar = onewnc.createVariable('ipoint', 'i')
108    newvar[:] = 0
109    ncvar.basicvardef(newvar, 'file_i', 'x-axis closest grid-point from file', '-')
110
111    newvar = onewnc.createVariable('jpoint', 'i')
112    newvar[:] = 0
113    ncvar.basicvardef(newvar, 'file_j', 'y-axis closest grid-point from file', '-')
114
115    onewnc.sync()
116
117    # Variables
118    NOvarvals = np.ones((dimt), dtype=np.float)*gen.fillValueF
119    for vn in sfcvars:
120        CFvalues = gen.variables_values(vn)
121        newvar = onewnc.createVariable(vn, 'f4', ('time'), fill_value=gen.fillValueF)
122        newvar[:] = NOvarvals[:]
123        ncvar.basicvardef(newvar, CFvalues[0], CFvalues[4].replace('|',' '),         \
124          CFvalues[5])
125
126    # Global attributes
127    ncvar.add_global_PyNCplot(onewnc, 'get_stations.py', fname, '0.1')
128    ncvar.set_attribute(newvar, 'data_origin', 'SMN')
129
130    onewnc.sync()
131    onewnc.close()
132
133    return
134
135def creation_sndstation_file(filen, lv, Lv, hv, lab, tunits, punits, ptime, sndvars, \
136  dimt, dimz, ifilen):
137    """ Function to create the structure of the surface station file
138      filen: name of the file
139      lv: value of the longitude of the station
140      Lv: value of the latitude of the station
141      hv: value of the height of the station
142      lab: label of the station
143      tunits: uints of time
144      punits: uints of pressure
145      ptime: does pressure evolves in time? (True/False)
146      sndvars: variables to be included
147      dimt: quantity of time-steps
148      dimz: quantity of vertical levels
149      ifilen: name of the file from which data is retrieved
150    """
151    fname = 'creation_sndstation_file'
152
153    Lstring = 256
154
155    onewnc = NetCDFFile(filen, 'w')
156    # Dimensions
157    newdim = onewnc.createDimension('time', None)
158    newdim = onewnc.createDimension('plev', dimz)
159    newdim = onewnc.createDimension('Lstring', Lstring)
160
161    # Variable-dimensions
162    newvar = onewnc.createVariable('time', 'f8', ('time'))
163    ncvar.basicvardef(newvar, 'time', 'Time', tunits)
164    ncvar.set_attribute(newvar, 'calendar', 'standard')
165    newvar.setncattr('axis', 'T')
166    newvar.setncattr('_CoordinateAxisType', 'Time')
167
168    if ptime:
169        newvar = onewnc.createVariable('plev', 'f8', ('time','plev'))
170        ncvar.basicvardef(newvar, 'plev', 'air pressure', punits)
171        ncvar.set_attribute(newvar, 'down', 'up')
172        newvar.setncattr('axis', 'Z')
173        newvar.setncattr('_CoordinateAxisType', 'pressure')
174    else:
175        newvar = onewnc.createVariable('plev', 'f8', ('plev'))
176        ncvar.basicvardef(newvar, 'plev', 'air pressure', punits)
177        ncvar.set_attribute(newvar, 'down', 'up')
178        newvar.setncattr('axis', 'Z')
179        newvar.setncattr('_CoordinateAxisType', 'pressure')
180    onewnc.sync()
181
182    # station Variables
183    Llab = len(lab)
184    newvar = onewnc.createVariable('station', 'c', ('Lstring'))
185    newvar[0:Llab] = lab[:]
186    ncvar.basicvardef(newvar, 'station', 'Name of the station', '-')
187
188    newvar = onewnc.createVariable('lon', 'f8')
189    newvar[:] = lv
190    ncvar.basicvardef(newvar, 'lon', 'Longitude', 'degrees_west')
191    newvar.setncattr('axis', 'X')
192    newvar.setncattr('_CoordinateAxisType', 'Lon')
193   
194    newvar = onewnc.createVariable('lat', 'f8')
195    newvar[:] = Lv
196    ncvar.basicvardef(newvar, 'lat', 'Latitude', 'degrees_north')
197    newvar.setncattr('axis', 'Y')
198    newvar.setncattr('_CoordinateAxisType', 'Lat')
199   
200    newvar = onewnc.createVariable('height', 'f8')
201    newvar[:] = hv
202    ncvar.basicvardef(newvar, 'height', 'Height', 'm')
203    newvar.setncattr('axis', 'Z')
204    newvar.setncattr('_CoordinateAxisType', 'Height')
205    onewnc.sync()
206
207    newvar = onewnc.createVariable('flon', 'f8')
208    newvar[:] = lv
209    ncvar.basicvardef(newvar, 'file_lon', 'Longitude closest grid-point from file',  \
210      'degrees_west')
211    newvar.setncattr('axis', 'X')
212    newvar.setncattr('_CoordinateAxisType', 'Lon')
213   
214    newvar = onewnc.createVariable('flat', 'f8')
215    newvar[:] = Lv
216    ncvar.basicvardef(newvar, 'file_lat', 'Latitude closest grid-point from file',   \
217      'degrees_north')
218    newvar.setncattr('axis', 'Y')
219    newvar.setncattr('_CoordinateAxisType', 'Lat')
220   
221    newvar = onewnc.createVariable('fheight', 'f8')
222    newvar[:] = hv
223    ncvar.basicvardef(newvar, 'file_height', 'Height closest grid-point from file',  \
224      'm')
225    newvar.setncattr('axis', 'Z')
226    newvar.setncattr('_CoordinateAxisType', 'Height')
227
228    newvar = onewnc.createVariable('ipoint', 'i')
229    newvar[:] = 0
230    ncvar.basicvardef(newvar, 'file_i', 'x-axis closest grid-point from file', '-')
231
232    newvar = onewnc.createVariable('jpoint', 'i')
233    newvar[:] = 0
234    ncvar.basicvardef(newvar, 'file_j', 'y-axis closest grid-point from file', '-')
235
236    onewnc.sync()
237
238    # Variables
239    NOvarvals = np.ones((dimt,dimz), dtype=np.float)*gen.fillValueF
240    for vn in sndvars:
241        if not onewnc.variables.has_key(vn):
242            CFvalues = gen.variables_values(vn)
243            newvar = onewnc.createVariable(vn, 'f4', ('time', 'plev'),               \
244              fill_value=gen.fillValueF)
245            newvar[:] = NOvarvals[:]
246            ncvar.basicvardef(newvar, CFvalues[0], CFvalues[4].replace('|',' '),     \
247              CFvalues[5])
248
249    # Global attributes
250    ncvar.add_global_PyNCplot(onewnc, 'get_stations.py', fname, '0.1')
251    ncvar.set_attribute(newvar, 'data_origin', 'SMN')
252
253    onewnc.sync()
254    onewnc.close()
255
256    return
257
258
259#######    #######
260## MAIN
261    #######
262
263Wcomputed = {}
264for vn in diag.Wavailablediags:
265    Wcomputed[vn] = False
266
267Ccomputed = {}
268for vn in diag.Cavailablediags:
269    Ccomputed[vn] = False
270
271yrref = ReferenceDate[0:4]
272monref = ReferenceDate[4:6]
273dayref = ReferenceDate[6:8]
274horref = ReferenceDate[8:10]
275minref = ReferenceDate[10:12]
276secref = ReferenceDate[12:14]
277
278utime = UnitsTime + ' since ' + yrref + '-' + monref + '-' + dayref + ' ' + horref + \
279  ':' + minref + ':' + secref
280
281# surface stations
282if sfcstatfile != 'None':
283    osfcstatf = open(sfcstatfile, 'r')
284
285    sfclonvals = []
286    sfclatvals = []
287    sfclabvals = []
288    sfchgtvals = []
289
290    for line in osfcstatf:
291        if line[0:1] != comchar and len(line) > 1:
292            vals = line.replace('\n','').split(sepchar)
293            if not gen.searchInlist(missvalS, vals[sfcloncol]) and                       \
294              not gen.searchInlist(missvalS, vals[sfclatcol]):
295                stlon = np.float(vals[sfcloncol])
296                stlat = np.float(vals[sfclatcol])
297                if stlon >= nlon and stlon <= xlon and stlat >= nlat and stlat <= xlat:
298                    sfclonvals.append(stlon)
299                    sfclatvals.append(stlat)
300                    sfclabvals.append(vals[sfclabcol])
301                    if gen.searchInlist(missvalS, vals[sfchgtcol]):
302                        if vals[sfchgtcol] == '0': 
303                            sfchgtvals.append(np.float(vals[sfchgtcol]))
304                        else: 
305                            sfchgtvals.append(gen.fillValueF)
306                    else:
307                        sfchgtvals.append(np.float(vals[sfchgtcol]))
308   
309    osfcstatf.close()
310    Nsfcstats = len(sfclonvals)
311    print '  Number of surface stations: ', Nsfcstats
312else:
313    Nsfcstats = 0
314    print '  No surface stations !!'
315
316# sounding stations
317if sndstatfile != 'None':
318    osndstatf = open(sndstatfile, 'r')
319
320    sndlonvals = []
321    sndlatvals = []
322    sndlabvals = []
323    sndhgtvals = []
324
325    for line in osndstatf:
326        if line[0:1] != comchar and len(line) > 1:
327            vals = line.replace('\n','').split(sepchar)
328            if vals[sndloncol] != missvalS and vals[sndlatcol] != missvalS:
329                stlon = np.float(vals[sndloncol])
330                stlat = np.float(vals[sndlatcol])
331                if stlon >= nlon and stlon <= xlon and stlat >= nlat and stlat <= xlat:
332                    sndlonvals.append(stlon)
333                    sndlatvals.append(stlat)
334                    sndlabvals.append(vals[sndlabcol])
335                    sndhgtvals.append(np.float(vals[sndhgtcol]))
336
337    osndstatf.close()
338    Nsndstats = len(sndlonvals)
339    print '  Number of sounding stations: ', Nsndstats
340else:
341    Nsndstats = 0
342    print '  No sounding stations !!'
343
344
345if not os.path.isfile(simfilen):
346    print gen.errormsg
347    print "  file '" + simfilen + "' does not exist !!"
348    quit(-1)
349
350onc = NetCDFFile(simfilen, 'r')
351olistv = onc.variables.keys()
352olistv.sort()
353
354# Getting basic values for each dimension
355dimvarvalues = {}
356for dimn in dimvariables.keys():
357    varn = dimvariables[dimn]
358    if not gen.searchInlist(NONcheckingvars, varn) and not gen.searchInlist(olistv, varn):
359        print gen.errormsg
360        print "  file '" + simfilen + "' does not have variable '" + varn + "' !!"
361        print '    available ones: ', olistv
362        quit(-1)
363
364    if not gen.searchInlist(NONcheckingvars, varn):
365        # Except for temporal variables, we don't want dimensions with time-axis
366        #  (assuming fixed)
367        if varn == 'WRFtime':
368            varvalues, CFtu = ncvar.compute_WRFtime(onc.variables['Times'],          \
369              ReferenceDate, UnitsTime)
370            dt = varvalues.shape[0]
371            dimvarvalues[varn] = varvalues
372        else:
373            ovar = onc.variables[varn]
374            slicevar = {}
375            if not gen.searchInlist(axesvars['T'], varn):
376                for dn in ovar.dimensions:
377                    if not gen.searchInlist(axesdims['T'], dn):
378                        slicevar[dn] = -1
379                    else:
380                        slicevar[dn] = 0
381            else:
382                slicevar[dn] = -1
383                dt = len(onc.dimensions[dn])
384            slicevar, vardims = ncvar.SliceVarDict(ovar, slicevar)
385            dimvarvalues[varn] = ovar[tuple(slicevar)]
386
387    if dimn == Pressdimref: dz = len(onc.dimensions[dimn])
388
389# Looking for 2D 'lon', 'lat' related variables
390for xn in axesvars['X']:
391    Xvarvals = dimvarvalues[dimvariables[xn]]
392    Yvarvals = dimvarvalues[CFdims['lat']]
393    if len(Xvarvals.shape) == 1:
394        newXvarvals, newYvarvals = np.meshgrid(Xvarvals, Yvarvals) 
395        dimvarvalues[dimvariables[dimvariables[xn]] = newXvarvals
396for yn in axesvars['Y']:
397    Xvarvals = dimvarvalues[CFdims['lon']]
398    Yvarvals = dimvarvalues[dimvariables[yn]]
399    if len(Yvarvals.shape) == 1:
400        newXvarvals, newYvarvals = np.meshgrid(Xvarvals, Yvarvals) 
401        dimvarvalues[dimvariables[dimvariables[yn]] = newYvarvals
402
403# Retrieving surface data
404##
405
406# Diagnosted variables
407diagvars = {}
408
409# Coinident surface station
410sfccoinc = {}
411
412for ist in range(Nsfcstats):
413    jumpstation = False
414
415    lonv = sfclonvals[ist]
416    latv = sfclatvals[ist]
417    labelv = sfclabvals[ist]
418    heightv = sfchgtvals[ist]
419
420    stfilen = 'sfc_station_' + labelv + '.nc'
421
422    creation_sfcstation_file(stfilen, lonv, latv, heightv, labelv, utime,            \
423      sfcvariables.keys(), dt, simfilen)
424
425    onewnc = NetCDFFile(stfilen, 'a')
426
427    stationsji = np.ones((Nsfcstats,2), dtype=int)*gen.fillValueI
428    coincstats = np.ones((Nsfcstats), dtype=int)*gen.fillValueI
429
430    for sfcv in sfcvariables.keys():
431        fvn = sfcvariables[sfcv]
432        if not gen.searchInlist(NONcheckingvars,fvn) and                             \
433          not onewnc.variables.has_key(sfcv):
434            print gen.errormsg
435            print "  Newfile for the station '" + stfilen + "' does not content " +  \
436              " variable '" + sfcv + "' !!"
437            print '    variables:', onewnc.variables.keys()
438            quit(-1) 
439
440        if not gen.searchInlist(NONcheckingvars,fvn):
441            ogetvar = onc.variables[fvn]
442            getdims = ogetvar.dimensions
443        else:
444            getdims = list(onc.variables[sfcrefvar].dimensions)
445            getvdims = [dimvariables[getdims[0]], dimvariables[getdims[1]],          \
446              dimvariables[getdims[2]]]
447
448        # Looking for the X,Y axis for location of station within file
449        for dn in getdims:
450            axis = gen.dictionary_key_list(axesdims, dn)
451            if not dimvariables.has_key(dn):
452                print gen.errormsg
453                print "  dimension '" + dn + "' not ready in 'dimvariables' !!"
454                print '    add it to proceed'
455                quit(-1)
456            if axis == 'X':
457                lvals = dimvarvalues[dimvariables[dn]]
458                xdiff = (lvals-lonv)**2
459            elif axis == 'Y':
460                Lvals = dimvarvalues[dimvariables[dn]]
461                ydiff = (Lvals-latv)**2
462
463        # Looking for the closest point
464        diff = np.sqrt(xdiff + ydiff)
465        mindiff = np.min(diff)
466        minji = gen.index_mat(diff, mindiff)
467       
468        coincstats = (stationsji[:,0] - minji[0]) + (stationsji[:,1] - minji[1])
469        if np.any(coincstats == 0):
470            print '  Surface station coincidence by closest grid point from file !!'
471            iist = gen.index_vec(coincstats, 0)
472            print '  ' + labelv, '&', sfclabvals[iist]
473            onewnc.close()
474            sub.call('rm ' + stfilen, shell=True)
475            sub.call('cp ' + 'sfc_station_' + sfclabvals[iist] + '.nc ' + stfilen,   \
476              shell=True)
477            onewnc = NetCDFFile(stfilen, 'a')
478            ostlab = onewnc.variables['station']
479            Llab = len(labelv)
480            ostlab[0:Llab] = labelv
481            onewnc.sync()
482            onewnc.close()
483            if sfccoinc.has_key(sfclabvals[iist]):
484                lvals = sfccoinc[sfclabvals[iist]]
485                lvals.append(labelv)
486                sfccoinc[sfclabvals[iist]] = lvals
487            else:
488                sfccoinc[sfclabvals[iist]] = [labelv]
489            continue
490       
491        # Writting information to file
492        if sfcv == sfcrefvar:
493            ojvar = onewnc.variables['jpoint']
494            ojvar[:] = minji[0]
495            ojvar = onewnc.variables['ipoint']
496            ojvar[:] = minji[1]
497            ohvar = onewnc.variables['fheight']
498            ohvar[:] = dimvarvalues[dimvariables['H']][minji[0],minji[1]]
499            onewnc.sync()
500
501        # Computing a single time the diagnosted variables
502        if gen.searchInlist(NONcheckingvars,fvn):
503            if ist == 0:
504                if fvn[0:3] == 'WRF':
505                    fvn_diag = fvn[3:len(fvn)]
506                    if not Wcomputed[fvn_diag]:
507                        # Using a given series of prepared diagnostics
508                        Wcomputed[fvn_diag] = True
509                        vardiag = diag.W_diagnostic(fvn_diag, onc,                  \
510                          sfcvariables, sndvariables, getdims, getvdims, CFdims)
511                        diagvars[fvn] = vardiag.values
512                elif fvn[0:2] == 'C_':
513                    fvn_diag = fvn.split('_')[1]
514                    if not Ccomputed[fvn_diag]:
515                        # Using a given series of prepared diagnostics
516                        Ccomputed[fvn_diag] = True
517                        vardiag = C_diagnostic(fvn_diag, onc, sfcvariables,          \
518                          sndvariables, CFdims)
519                        diagvars[fvn] = vardiag.values
520                ogetvar = diagvars[fvn]
521            else:
522                ogetvar = diagvars[fvn]
523               
524        slicevar = []
525        for dn in getdims:
526            axis = gen.dictionary_key_list(axesdims, dn)
527            if axis == 'X': slicevar.append(minji[1])
528            elif axis == 'Y': slicevar.append(minji[0])
529            else: slicevar.append(slice(0,len(onc.dimensions[dn])))
530
531        # Writting data and slicing
532        onewvar = onewnc.variables[sfcv]
533        onewvar[:] = ogetvar[tuple(slicevar)]   
534        onewnc.sync()
535
536    if jumpstation: continue
537
538    onewnc.close()
539
540if len(sfccoinc.keys()) > 0:
541    print 'Coincident surface stations _______'
542    gen.printing_dictionary(sfccoinc)
543
544# Retrieving sounding data
545
546# Diagnosted variables
547diagvars = {}
548
549# Coinident sounding station
550sndcoinc = {}
551
552for ist in range(Nsndstats):
553    jumpstation = False
554
555    lonv = sndlonvals[ist]
556    latv = sndlatvals[ist]
557    labelv = sndlabvals[ist]
558    heightv = sndhgtvals[ist]
559
560    stfilen = 'snd_station_' + labelv + '.nc'
561
562    creation_sndstation_file(stfilen, lonv, latv, heightv, labelv, utime, UnitsPress,\
563      presstime, sndvariables.keys(), dt, dz, simfilen)
564
565    onewnc = NetCDFFile(stfilen, 'a')
566
567    stationsji = np.ones((Nsfcstats,2), dtype=int)*gen.fillValueI
568    coincstats = np.ones((Nsfcstats), dtype=int)*gen.fillValueI
569
570    for sndv in sndvariables.keys():
571        fvn = sndvariables[sndv]
572        if not gen.searchInlist(NONcheckingvars,fvn) and                             \
573          not onewnc.variables.has_key(sndv):
574            print gen.errormsg
575            print "  Newfile for the station '" + stfilen + "' does not content " +  \
576              " variable '" + sndv + "' !!"
577            print '    variables:', onewnc.variables.keys()
578            quit(-1) 
579
580        if not gen.searchInlist(NONcheckingvars,fvn):
581            ogetvar = onc.variables[fvn]
582            getdims = ogetvar.dimensions
583        else:
584            getdims = list(onc.variables[sndrefvar].dimensions)
585            getvdims = [dimvariables[getdims[0]], dimvariables[getdims[1]],          \
586              dimvariables[getdims[2]], dimvariables[getdims[3]]]
587
588        # Looking for the X,Y axis for location of station within file
589        for dn in getdims:
590            axis = gen.dictionary_key_list(axesdims, dn)
591            if not dimvariables.has_key(dn):
592                print gen.errormsg
593                print "  dimension '" + dn + "' not ready in 'dimvariables' !!"
594                print '    add it to proceed'
595                quit(-1)
596            if axis == 'X':
597                lvals = dimvarvalues[dimvariables[dn]]
598                xdiff = (lvals-lonv)**2
599            elif axis == 'Y':
600                Lvals = dimvarvalues[dimvariables[dn]]
601                ydiff = (Lvals-latv)**2
602
603        # Looking for the closest point
604        diff = np.sqrt(xdiff + ydiff)
605        mindiff = np.min(diff)
606        minji = gen.index_mat(diff, mindiff)
607       
608        coincstats = (stationsji[:,0] - minji[0]) + (stationsji[:,1] - minji[1])
609        if np.any(coincstats == 0):
610            print '  Sounding station coincidence by closest grid point from file !!'
611            iist = gen.index_vec(coincstats, 0)
612            print '  ' + labelv, '&', sndlabvals[iist]
613            onewnc.close()
614            sub.call('rm ' + stfilen, shell=True)
615            sub.call('cp ' + 'snd_station_' + sndlabvals[iist] + '.nc ' + stfilen,   \
616              shell=True)
617            onewnc = NetCDFFile(stfilen, 'a')
618            ostlab = onewnc.variables['station']
619            Llab = len(labelv)
620            ostlab[0:Llab] = labelv
621            onewnc.sync()
622            onewnc.close()
623            if sndcoinc.has_key(sndlabvals[iist]):
624                lvals = sfccoinc[sndlabvals[iist]]
625                lvals.append(labelv)
626                sndcoinc[sndlabvals[iist]] = lvals
627            else:
628                sndcoinc[sndlabvals[iist]] = [labelv]
629            continue
630       
631        # Writting information to file
632        if sndv == sndrefvar:
633            ojvar = onewnc.variables['jpoint']
634            ojvar[:] = minji[0]
635            ojvar = onewnc.variables['ipoint']
636            ojvar[:] = minji[1]
637            if dimvariables['H'] != 'None':
638                ohvar = onewnc.variables['fheight']
639                ohvar[:] = dimvarvalues[dimvariables['H']][minji[0],minji[1]]
640            onewnc.sync()
641
642        # Computing a single time the diagnosted variables
643        if gen.searchInlist(NONcheckingvars,fvn):
644            if ist == 0:
645                if fvn[0:3] == 'WRF':
646                    fvn_diag = fvn[3:len(fvn)]
647                    if not Wcomputed[fvn_diag]:
648                        # Using a given series of prepared diagnostics
649                        Wcomputed[fvn_diag] = True
650                        vardiag = diag.W_diagnostic(fvn_diag, onc,                   \
651                          sfcvariables, sndvariables, getdims, getvdims, CFdims)
652                        diagvars[fvn] = vardiag.values
653                elif fvn[0:2] == 'C_':
654                    fvn_diag = fvn.split('_')[1]
655                    if not Ccomputed[fvn_diag]:
656                        # Using a given series of prepared diagnostics
657                        Ccomputed[fvn_diag] = True
658                        vardiag = diag.C_diagnostic(fvn_diag, onc, sfcvariables,     \
659                          sndvariables, CFdims)
660                        diagvars[fvn] = vardiag.values
661
662                ogetvar = diagvars[fvn]
663            else:
664                ogetvar = diagvars[fvn]
665               
666        slicevar = []
667        for dn in getdims:
668            axis = gen.dictionary_key_list(axesdims, dn)
669            if axis == 'X': slicevar.append(minji[1])
670            elif axis == 'Y': slicevar.append(minji[0])
671            else: slicevar.append(slice(0,len(onc.dimensions[dn])))
672
673        # Writting data and slicing
674        onewvar = onewnc.variables[sndv]
675        print 'Lluis sndv:', sndv, ' shapes onewvar:', onewvar.shape, 'slice:', slicevar
676        onewvar[:] = ogetvar[tuple(slicevar)]   
677        onewnc.sync()
678
679    if jumpstation: continue
680
681    onewnc.close()
682
683if len(sndcoinc.keys()) > 0:
684    print 'Coincident sounding stations _______'
685    gen.printing_dictionary(sndcoinc)
686
Note: See TracBrowser for help on using the repository browser.