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

Last change on this file since 1695 was 1694, checked in by lfita, 8 years ago

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