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

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

Fixing and working copy of the script
File size: 24.4 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[xn]] = newXvarvals
396        if len(Yvarvals.shape) == 1:
397            dimvarvalues[dimvariables['lat']] = newYvarvals
398
399for yn in axesvars['Y']:
400    Xvarvals = dimvarvalues[CFdims['lon']]
401    Yvarvals = dimvarvalues[dimvariables[yn]]
402    if len(Yvarvals.shape) == 1:
403        newXvarvals, newYvarvals = np.meshgrid(Xvarvals, Yvarvals) 
404        dimvarvalues[dimvariables[yn]] = newYvarvals
405
406# Retrieving surface data
407##
408
409# Diagnosted variables
410diagvars = {}
411
412# Coinident surface station
413sfccoinc = {}
414
415for ist in range(Nsfcstats):
416    jumpstation = False
417
418    lonv = sfclonvals[ist]
419    latv = sfclatvals[ist]
420    labelv = sfclabvals[ist]
421    heightv = sfchgtvals[ist]
422
423    stfilen = 'sfc_station_' + labelv + '.nc'
424
425    creation_sfcstation_file(stfilen, lonv, latv, heightv, labelv, utime,            \
426      sfcvariables.keys(), dt, simfilen)
427
428    onewnc = NetCDFFile(stfilen, 'a')
429
430    stationsji = np.ones((Nsfcstats,2), dtype=int)*gen.fillValueI
431    coincstats = np.ones((Nsfcstats), dtype=int)*gen.fillValueI
432
433    for sfcv in sfcvariables.keys():
434        fvn = sfcvariables[sfcv]
435        if not gen.searchInlist(NONcheckingvars,fvn) and                             \
436          not onewnc.variables.has_key(sfcv):
437            print gen.errormsg
438            print "  Newfile for the station '" + stfilen + "' does not content " +  \
439              " variable '" + sfcv + "' !!"
440            print '    variables:', onewnc.variables.keys()
441            quit(-1) 
442
443        if not gen.searchInlist(NONcheckingvars,fvn):
444            ogetvar = onc.variables[fvn]
445            getdims = ogetvar.dimensions
446        else:
447            getdims = list(onc.variables[sfcrefvar].dimensions)
448            getvdims = [dimvariables[getdims[0]], dimvariables[getdims[1]],          \
449              dimvariables[getdims[2]]]
450
451        # Looking for the X,Y axis for location of station within file
452        for dn in getdims:
453            axis = gen.dictionary_key_list(axesdims, dn)
454            if not dimvariables.has_key(dn):
455                print gen.errormsg
456                print "  dimension '" + dn + "' not ready in 'dimvariables' !!"
457                print '    add it to proceed'
458                quit(-1)
459            if axis == 'X':
460                lvals = dimvarvalues[dimvariables[dn]]
461                xdiff = (lvals-lonv)**2
462            elif axis == 'Y':
463                Lvals = dimvarvalues[dimvariables[dn]]
464                ydiff = (Lvals-latv)**2
465
466        # Looking for the closest point
467        diff = np.sqrt(xdiff + ydiff)
468        mindiff = np.min(diff)
469        minji = gen.index_mat(diff, mindiff)
470       
471        coincstats = (stationsji[:,0] - minji[0]) + (stationsji[:,1] - minji[1])
472        if np.any(coincstats == 0):
473            print '  Surface station coincidence by closest grid point from file !!'
474            iist = gen.index_vec(coincstats, 0)
475            print '  ' + labelv, '&', sfclabvals[iist]
476            onewnc.close()
477            sub.call('rm ' + stfilen, shell=True)
478            sub.call('cp ' + 'sfc_station_' + sfclabvals[iist] + '.nc ' + stfilen,   \
479              shell=True)
480            onewnc = NetCDFFile(stfilen, 'a')
481            ostlab = onewnc.variables['station']
482            Llab = len(labelv)
483            ostlab[0:Llab] = labelv
484            onewnc.sync()
485            onewnc.close()
486            if sfccoinc.has_key(sfclabvals[iist]):
487                lvals = sfccoinc[sfclabvals[iist]]
488                lvals.append(labelv)
489                sfccoinc[sfclabvals[iist]] = lvals
490            else:
491                sfccoinc[sfclabvals[iist]] = [labelv]
492            continue
493       
494        # Writting information to file
495        if sfcv == sfcrefvar:
496            ojvar = onewnc.variables['jpoint']
497            ojvar[:] = minji[0]
498            ojvar = onewnc.variables['ipoint']
499            ojvar[:] = minji[1]
500            ohvar = onewnc.variables['fheight']
501            ohvar[:] = dimvarvalues[dimvariables['H']][minji[0],minji[1]]
502            onewnc.sync()
503
504        # Computing a single time the diagnosted variables
505        if gen.searchInlist(NONcheckingvars,fvn):
506            if ist == 0:
507                if fvn[0:3] == 'WRF':
508                    fvn_diag = fvn[3:len(fvn)]
509                    if not Wcomputed[fvn_diag]:
510                        # Using a given series of prepared diagnostics
511                        Wcomputed[fvn_diag] = True
512                        vardiag = diag.W_diagnostic(fvn_diag, onc,                  \
513                          sfcvariables, sndvariables, getdims, getvdims, CFdims)
514                        diagvars[fvn] = vardiag.values
515                elif fvn[0:2] == 'C_':
516                    fvn_diag = fvn.split('_')[1]
517                    if not Ccomputed[fvn_diag]:
518                        # Using a given series of prepared diagnostics
519                        Ccomputed[fvn_diag] = True
520                        vardiag = C_diagnostic(fvn_diag, onc, sfcvariables,          \
521                          sndvariables, CFdims)
522                        diagvars[fvn] = vardiag.values
523                ogetvar = diagvars[fvn]
524            else:
525                ogetvar = diagvars[fvn]
526               
527        slicevar = []
528        for dn in getdims:
529            axis = gen.dictionary_key_list(axesdims, dn)
530            if axis == 'X': slicevar.append(minji[1])
531            elif axis == 'Y': slicevar.append(minji[0])
532            else: slicevar.append(slice(0,len(onc.dimensions[dn])))
533
534        # Writting data and slicing
535        onewvar = onewnc.variables[sfcv]
536        onewvar[:] = ogetvar[tuple(slicevar)]   
537        onewnc.sync()
538
539    if jumpstation: continue
540
541    onewnc.close()
542
543if len(sfccoinc.keys()) > 0:
544    print 'Coincident surface stations _______'
545    gen.printing_dictionary(sfccoinc)
546
547# Retrieving sounding data
548
549# Diagnosted variables
550diagvars = {}
551
552# Coinident sounding station
553sndcoinc = {}
554
555for ist in range(Nsndstats):
556    jumpstation = False
557
558    lonv = sndlonvals[ist]
559    latv = sndlatvals[ist]
560    labelv = sndlabvals[ist]
561    heightv = sndhgtvals[ist]
562
563    stfilen = 'snd_station_' + labelv + '.nc'
564
565    creation_sndstation_file(stfilen, lonv, latv, heightv, labelv, utime, UnitsPress,\
566      presstime, sndvariables.keys(), dt, dz, simfilen)
567
568    onewnc = NetCDFFile(stfilen, 'a')
569
570    stationsji = np.ones((Nsfcstats,2), dtype=int)*gen.fillValueI
571    coincstats = np.ones((Nsfcstats), dtype=int)*gen.fillValueI
572
573    for sndv in sndvariables.keys():
574        fvn = sndvariables[sndv]
575        if not gen.searchInlist(NONcheckingvars,fvn) and                             \
576          not onewnc.variables.has_key(sndv):
577            print gen.errormsg
578            print "  Newfile for the station '" + stfilen + "' does not content " +  \
579              " variable '" + sndv + "' !!"
580            print '    variables:', onewnc.variables.keys()
581            quit(-1) 
582
583        if not gen.searchInlist(NONcheckingvars,fvn):
584            ogetvar = onc.variables[fvn]
585            getdims = ogetvar.dimensions
586        else:
587            getdims = list(onc.variables[sndrefvar].dimensions)
588            getvdims = [dimvariables[getdims[0]], dimvariables[getdims[1]],          \
589              dimvariables[getdims[2]], dimvariables[getdims[3]]]
590
591        # Looking for the X,Y axis for location of station within file
592        for dn in getdims:
593            axis = gen.dictionary_key_list(axesdims, dn)
594            if not dimvariables.has_key(dn):
595                print gen.errormsg
596                print "  dimension '" + dn + "' not ready in 'dimvariables' !!"
597                print '    add it to proceed'
598                quit(-1)
599            if axis == 'X':
600                lvals = dimvarvalues[dimvariables[dn]]
601                xdiff = (lvals-lonv)**2
602            elif axis == 'Y':
603                Lvals = dimvarvalues[dimvariables[dn]]
604                ydiff = (Lvals-latv)**2
605
606        # Looking for the closest point
607        diff = np.sqrt(xdiff + ydiff)
608        mindiff = np.min(diff)
609        minji = gen.index_mat(diff, mindiff)
610       
611        coincstats = (stationsji[:,0] - minji[0]) + (stationsji[:,1] - minji[1])
612        if np.any(coincstats == 0):
613            print '  Sounding station coincidence by closest grid point from file !!'
614            iist = gen.index_vec(coincstats, 0)
615            print '  ' + labelv, '&', sndlabvals[iist]
616            onewnc.close()
617            sub.call('rm ' + stfilen, shell=True)
618            sub.call('cp ' + 'snd_station_' + sndlabvals[iist] + '.nc ' + stfilen,   \
619              shell=True)
620            onewnc = NetCDFFile(stfilen, 'a')
621            ostlab = onewnc.variables['station']
622            Llab = len(labelv)
623            ostlab[0:Llab] = labelv
624            onewnc.sync()
625            onewnc.close()
626            if sndcoinc.has_key(sndlabvals[iist]):
627                lvals = sfccoinc[sndlabvals[iist]]
628                lvals.append(labelv)
629                sndcoinc[sndlabvals[iist]] = lvals
630            else:
631                sndcoinc[sndlabvals[iist]] = [labelv]
632            continue
633       
634        # Writting information to file
635        if sndv == sndrefvar:
636            ojvar = onewnc.variables['jpoint']
637            ojvar[:] = minji[0]
638            ojvar = onewnc.variables['ipoint']
639            ojvar[:] = minji[1]
640            if dimvariables['H'] != 'None':
641                ohvar = onewnc.variables['fheight']
642                ohvar[:] = dimvarvalues[dimvariables['H']][minji[0],minji[1]]
643            onewnc.sync()
644
645        # Computing a single time the diagnosted variables
646        if gen.searchInlist(NONcheckingvars,fvn):
647            if ist == 0:
648                if fvn[0:3] == 'WRF':
649                    fvn_diag = fvn[3:len(fvn)]
650                    if not Wcomputed[fvn_diag]:
651                        # Using a given series of prepared diagnostics
652                        Wcomputed[fvn_diag] = True
653                        vardiag = diag.W_diagnostic(fvn_diag, onc,                   \
654                          sfcvariables, sndvariables, getdims, getvdims, CFdims)
655                        diagvars[fvn] = vardiag.values
656                elif fvn[0:2] == 'C_':
657                    fvn_diag = fvn.split('_')[1]
658                    if not Ccomputed[fvn_diag]:
659                        # Using a given series of prepared diagnostics
660                        Ccomputed[fvn_diag] = True
661                        vardiag = diag.C_diagnostic(fvn_diag, onc, sfcvariables,     \
662                          sndvariables, CFdims)
663                        diagvars[fvn] = vardiag.values
664
665                ogetvar = diagvars[fvn]
666            else:
667                ogetvar = diagvars[fvn]
668               
669        slicevar = []
670        for dn in getdims:
671            axis = gen.dictionary_key_list(axesdims, dn)
672            if axis == 'X': slicevar.append(minji[1])
673            elif axis == 'Y': slicevar.append(minji[0])
674            else: slicevar.append(slice(0,len(onc.dimensions[dn])))
675
676        # Writting data and slicing
677        onewvar = onewnc.variables[sndv]
678        print 'Lluis sndv:', sndv, ' shapes onewvar:', onewvar.shape, 'slice:', slicevar
679        onewvar[:] = ogetvar[tuple(slicevar)]   
680        onewnc.sync()
681
682    if jumpstation: continue
683
684    onewnc.close()
685
686if len(sndcoinc.keys()) > 0:
687    print 'Coincident sounding stations _______'
688    gen.printing_dictionary(sndcoinc)
689
Note: See TracBrowser for help on using the repository browser.