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

Last change on this file since 1682 was 1679, checked in by lfita, 8 years ago

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