source: lmdz_wrf/trunk/tools/ORforcing_reconstruct.py @ 2453

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

New version for matrix reconstruction taking as reference vector values
File size: 11.6 KB
Line 
1# Reconstructing ORCHIDEE's forcing files as matrices for a gien area
2# L. Fita, CIMA. December 2017
3#
4import numpy as np
5from netCDF4 import Dataset as NetCDFFile
6import os
7import re
8import numpy.ma as ma
9# Importing generic tools file 'generic_tools.py'
10import generic_tools as gen
11import nc_var_tools as ncvar
12import subprocess as sub
13import module_ForSci as Sci
14from optparse import OptionParser
15
16parser = OptionParser()
17parser.add_option("-F", "--filename", dest="fn", help="name of files (overwrites -f option)", \
18  metavar="VALUE")
19parser.add_option("-f", "--fileHEader", dest="fh", help="header of files",           \
20  metavar="VALUE")
21parser.add_option("-i", "--indices", dest="indn", help="name of the variable indices", \
22  metavar="VALUE")
23parser.add_option("-L", "--latitude", dest="latn", help="name of the variable latitiude", \
24  metavar="VALUE")
25parser.add_option("-l", "--longitude", dest="lonn", help="name of the variable longitiude", \
26  metavar="VALUE")
27parser.add_option("-t", "--TransposeVariables", dest="tvars", help="whether variables should be trasposed", \
28  metavar="VALUE")
29parser.add_option("-v", "--Variables", dest="varns", help="',' separated list of variables", \
30  metavar="VALUE")
31parser.add_option("-y", "--year", dest="year", help="year to process",               \
32  metavar="VALUE")
33
34(opts, args) = parser.parse_args()
35
36####### ###### ##### #### ### ## #
37
38if opts.fn is None:
39    filen = opts.fh + opts.year + '.nc'
40else:
41    filen = opts.fn
42
43# Variable whcih provides the indices of a 1D vector from the dimy, dimx space
44indvar = opts.indn
45
46# 2D longitude, latitude matrices
47lonvar = opts.lonn
48latvar = opts.latn
49
50# Range to retrieve
51Xmin=-90.25
52Xmin='all'
53Xmax=-33.25
54Ymin=-67.25
55Ymax=15.25
56
57# Variables to reconstruct
58#variable = 'all'
59variable = opts.varns
60
61# Resolution
62resX = 0.5
63resY = 0.5
64
65# Minimum difference from matrix to localized point
66maxdiff = 0.05
67
68# Projection of the matrix
69matProj = 'latlon'
70
71#######    #######
72## MAIN
73    #######
74main = 'ORforcing_reconstruct.py'
75fname = 'ORforcing_reconstruct.py'
76
77onc = NetCDFFile(filen, 'r')
78
79availProj = ['latlon']
80
81ofilen = 'reconstruct_matrix_' + opts.year + '.nc'
82
83ncvars = onc.variables.keys()
84
85varstocheck = [indvar, lonvar, latvar]
86for vn in varstocheck:
87    if not gen.searchInlist(ncvars, vn):
88        print gen.errormsg
89        print '  ' + main + ": file '" + filen + "' does not have variable '" + vn + \
90          "' !!"
91        print '    available ones:', ncvars
92        quit(-1)
93
94oind = onc.variables[indvar]
95olon = onc.variables[lonvar]
96olat = onc.variables[latvar]
97
98indv = oind[:]
99lonv = olon[:]
100latv = olat[:]
101
102if len(olon.dimensions) == 2:
103    dimx = lonv.shape[1]
104    dimy = lonv.shape[0]
105    veclon = lonv.reshape(dimx*dimy)
106    veclat = latv.reshape(dimx*dimy)
107    Xdimn = olon.dimensions[1]
108    Ydimn = olon.dimensions[0]
109else:
110    veclon = lonv.copy()
111    veclat = latv.copy()
112    Xdimn = 'x'
113    Ydimn = 'y'
114
115vecdimn = oind.dimensions[0]
116
117if not gen.searchInlist(availProj, matProj):
118    print errormsg
119    print '  ' + fname + ": projection '" + matProj + "' not available !!"
120    print '    available ones:', availProj
121    quit(-1)
122
123ncdims = onc.dimensions
124
125if variable == 'all':
126    varns = ncvars
127else:
128    varns = variable.split(',')
129    for vn in varns:
130        if not gen.searchInlist(ncvars, vn):
131            print gen.errormsg
132            print '  ' + fname + ": file '" + filen + "' does not have " +           \
133              " variable '" + vn + "' !!"
134            print '  available ones:', ncvars
135            quit(-1)
136
137if gen.searchInlist(olon.ncattrs(), 'units'):
138    xunits = olon.getncattr('units')
139else:
140    xunits = '-'
141if gen.searchInlist(olat.ncattrs(), 'units'):
142    yunits = olat.getncattr('units')
143else:
144   yunits = '-'
145
146if type(Xmin) == type('2') and Xmin == 'all':
147    Xmin = np.min(lonv)
148    Xmax = np.max(lonv)
149    Ymin = np.min(latv)
150    Ymax = np.max(latv)
151
152# Matrix values
153if matProj == 'latlon':
154    dimx = int((Xmax - Xmin+resX)/resX)
155    dimy = int((Ymax - Ymin+resY)/resY)
156
157print 'Xmin:', Xmin, 'Xmax:', Xmax, 'Ymin:', Ymin, 'Ymax:', Ymax, 'maxdiff:', maxdiff
158
159matindt, matXt, matYt, matdifft = Sci.module_scientific.reconstruct_matrix(          \
160  vectorxpos=veclon, vectorypos=veclat, dvec=veclon.shape[0], xmin=Xmin, xmax=Xmax,  \
161  ymin=Ymin,ymax=Ymax, dmatx=dimx, dmaty=dimy, matproj=matProj, maxdiff=maxdiff)
162
163matind = matindt.transpose()
164Nfound = np.sum(matind != -1)
165Nstations = veclon.shape[0]
166print '  Nfound:', Nfound, ' number of stations:', Nstations
167
168if Nfound*1. / Nstations < 0.8:
169    print gen.errormsg
170    print '  '+main + ': only ', '{:.2f}'.format(Nfound*100./Nstations),\
171      '% of points ' + 'have been found !!'
172    print '    this is not enough. Something must went wrong!'
173    print '    Longitudes Latitudes _______'
174    for i in range(veclon.shape[0]):
175        print '    ', veclon[i], veclat[i]
176    dx2 = dimx/2
177    dy2 = dimy/2
178    print '    dx2, dy2 -/+ 5 fraction of destiny longitudes _______'
179    for j in range(-5,5):
180        print matXt[dx2-5:dx2+5,dy2+j]
181    print '    dx2, dy2 -/+ 5 fraction of destiny latitudes _______'
182    for j in range(-5,5):
183        print matYt[dx2-5:dx2+5,dy2+j]
184    print '    dx2, dy2 -/+ 5 fraction of indices equivalency _______'
185    for j in range(-5,5):
186        print matindt[dx2-5:dx2+5,dy2+j]
187    print '     min distance lon(dy2,dx2)=', matXt[dx2,dy2], ':',                    \
188      np.min(veclon - matXt[dx2,dy2])
189    print '     min distance lat(dy2,dx2)=', matYt[dx2,dy2], ':',                    \
190      np.min(veclat - matYt[dx2,dy2])
191    print '     longitude borders:', Xmin, Xmax, 'dX:', resX
192    print '     latitude borders:', Ymin, Ymax, 'dY:', resY
193    #quit(-1)
194
195# Fortran like, First 1
196matind = np.where(matind != -1, matind - 1, matind)
197
198# Creation of file
199onewnc = NetCDFFile(ofilen, 'w')
200
201# Dimensions
202newdim = onewnc.createDimension('x', dimx)
203newdim = onewnc.createDimension('y', dimy)
204
205# Variable-dimension
206newvar = onewnc.createVariable('lon', 'f8', ('y', 'x'))
207newvar[:] = matXt.transpose()
208ncvar.basicvardef(newvar, 'lon', 'Longitude', 'degrees_east')   
209
210newvar = onewnc.createVariable('lat', 'f8', ('y', 'x'))
211newvar[:] = matYt.transpose()
212ncvar.basicvardef(newvar, 'lat', 'Latitude', 'degrees_north')   
213onewnc.sync()
214
215# Variable indices
216newvar = onewnc.createVariable('vec1D_matind', 'i', ('y', 'x'), fill_value=-1)
217newvar[:] = matind
218ncvar.basicvardef(newvar, 'vec1D_matind', 'matrix with the equivalencies from 1D ' + \
219  'vector indices', '-')
220ncvar.set_attribute(newvar, 'coordinates', 'lon lat')
221
222# Variable differences
223newvar = onewnc.createVariable('vec1D_matdiff', 'i', ('y', 'x'), fill_value=-1)
224newvar[:] = matdifft.transpose()
225ncvar.basicvardef(newvar,'vec1D_matdiff', 'matrix differences respect 1D ', 'degrees')
226ncvar.set_attribute(newvar, 'coordinates', 'lon lat')
227
228# Looking for equivalencies in the 1D vector
229matlonlat = matind.copy()
230for j in range(dimy):
231    for i in range(dimx):
232        if matind[j,i] != -1:
233            matlonlat[j,i] = indv[matind[j,i]]
234
235newvar = onewnc.createVariable('lonlat_matind', 'i', ('y', 'x'), fill_value=-1)
236newvar[:] = matlonlat
237ncvar.basicvardef(newvar, 'lonlat_matind', 'matrix with the equivalencies from ' +   \
238  '2D lon, lat matrices', '-')
239ncvar.set_attribute(newvar, 'coordinates', 'lon lat')
240onewnc.sync()
241
242# Getting variables
243for vn in varns:
244    if not onewnc.variables.has_key(vn):
245        ovar = onc.variables[vn]
246        if gen.Str_Bool(opts.tvars):
247            indn0 = ovar.dimensions
248            indn = list(indn0)[::-1]
249        else:
250            indn = ovar.dimensions
251        vardims = []
252        shapevar = []
253        for dn in indn:
254            if not gen.searchInlist(onewnc.dimensions, dn) and dn != Xdimn and   \
255              dn != Ydimn:
256                if onc.dimensions[dn].isunlimited():
257                    newdim = onewnc.createDimension(dn, None)
258                else:
259                    newdim = onewnc.createDimension(dn, len(onc.dimensions[dn]))
260               
261            if dn == vecdimn: 
262                vardims.append('y')
263                vardims.append('x')
264                shapevar.append(dimy)
265                shapevar.append(dimx)
266            else: 
267                vardims.append(dn)
268                shapevar.append(len(onc.dimensions[dn]))
269
270        if ovar.dtype == type(int(2)):
271            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
272              fill_value=gen.fillValueI)
273            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI
274        elif ovar.dtype == type(np.int32(2)):
275            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
276              fill_value=gen.fillValueI)
277            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI
278        elif ovar.dtype == type(np.int64(2)):
279            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
280              fill_value=gen.fillValueI)
281            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueI
282        elif ovar.dtype == type(np.float(2.)):
283            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
284              fill_value=gen.fillValueF)
285            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueF
286        elif ovar.dtype == type(np.float32(2.)):
287            newvar= onewnc.createVariable(vn,ncvar.nctype(ovar.dtype),tuple(vardims),\
288              fill_value=gen.fillValueF)
289            varvals = np.ones(tuple(shapevar), dtype=ovar.dtype)*gen.fillValueF
290        else:
291            print gen.errormsg
292            print '  ' + fname + ': variable type:', ovar.dtype, ' not ready !!'
293            quit(-1)
294
295        print '  reconstructing:', vn, ' shape:', newvar.shape, '...'
296        # Filling variable. It would be faster if we can avoid this loop... I'm feeling lazy!
297        if not gen.searchInlist(vardims,'x') and not gen.searchInlist(vardims,'y'):
298            if gen.Str_Bool(opts.tvars):
299                newvar[:] = ovar[:].transpose() 
300            else:
301                newvar[:] = ovar[:] 
302        else:
303            if gen.Str_Bool(opts.tvars):
304                ovart = ovar[:]
305            else:
306                ovart = ovar[:].transpose()
307            print '  Lluis shapes ovart:', ovart.shape, 'newvar:', newvar.shape
308            if newvar.dtype == type(float(2.)) or newvar.dtype == type(np.float(2.)) \
309              or newvar.dtype == type(np.float32(2)) or                              \
310              newvar.dtype == type(np.float64(2)):
311                newvals = Sci.module_scientific.fill3dr_2dvec(matind=matindt,        \
312                  inmat=ovart, id1=ovart.shape[0], id2=ovart.shape[1],               \
313                  od1=newvar.shape[2], od2=newvar.shape[1], od3=newvar.shape[0])
314            else:
315                newvals = Sci.module_scientific.fill3di_2dvec(matind=matindt,        \
316                  inmat=ovart, id1=ovart.shape[0], id2=ovart.shape[1],               \
317                  od1=newvar.shape[2], od2=newvar.shape[1], od3=newvar.shape[0])
318            newvar[:] = newvals.transpose()
319
320        # Attributes
321        for atn in ovar.ncattrs():
322            if atn != '_FillValue' and atn != 'units':
323                atv = ovar.getncattr(atn)
324                ncvar.set_attribute(newvar, atn, atv)
325        ncvar.set_attribute(newvar, 'coordinates', 'lon lat')
326        onewnc.sync()
327   
328# Global attributes
329for atn in onc.ncattrs():
330    atv = onc.getncattr(atn)
331    ncvar.set_attribute(onewnc, atn, atv)
332onewnc.sync()
333ncvar.add_global_PyNCplot(onewnc, main, fname, '0.1')
334onc.close()
335
336# Reconstructing times
337#otime = onewnc.variables['time']
338#ncvar.set_attribute(otime, 'units', 'seconds since ' + opts.year + '-01-01 00:00:00')
339
340onewnc.sync()
341onewnc.close()
342
343print fname + ": Successful writing of file '" + ofilen + ".nc' !!"
Note: See TracBrowser for help on using the repository browser.