source: lmdz_wrf/trunk/tools/get_stations.py

Last change on this file was 1725, checked in by lfita, 8 years ago

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