1 | # Python script to rectify ORCHIDEE netCDF files forcings |
---|
2 | # L. Fita, LMD. April 2016 |
---|
3 | # |
---|
4 | # More about: |
---|
5 | # https://forge.ipsl.jussieu.fr/orchidee/wiki/Documentation/Forcings/CfStandard |
---|
6 | # |
---|
7 | ## e.g. # file_rectifier.py -f alb_bg_jrctip.nc -v bg_albedo -d longitude,latitude -r 0.,1. -l longitude,latitude -m maskvarname@mask -o alb_bg_jrctip_XIOS.nc -T None,range -t alb_bg_jrctip_types.txt -I None,time_counter -Y monthly@means,months@since@0001-01-01@00:00:00,standard -S albedo -L background@albedo@from@MODIS -U '1' -A 'dataset:MODIS!albedo!with!the!JRC-TIP!package:U,article:Pinty et al. 2011:U' |
---|
8 | ## e.g. # file_rectifier.py -f lai2D.nc -v LAI -d longitude,latitude -r 0.,7. -l longitude,latitude -m maskabove@0. -o lai2D_XIOS.nc -T None,PFT -t PFTmap_types.txt -I None,time_counter -Y monthly@means,months@since@0001-01-01@00:00:00,standard -S lai -L leaf@area@index -U '-' -z 0. |
---|
9 | ## e.g. # file_rectifier.py -f PFTmap.nc -v maxvegetfrac -r 0.,1. -d lon,lat -l lon,lat -m masksumrange@0.9999999,1.e-8,2. -o PFTmap_XIOS.nc -T None,veget -I time_counter,time_counter -t PFTmap_types.txt -Y time,years@since@0001-01-01@00:00:00,standard -S maxvegetfrac |
---|
10 | ## e.g. # file_rectifier.py -f soils_param.nc -v soilcolor -r 1,8 -d x,y -l nav_lon,nav_lat -m maskabove@1.e-8 -o soils_param_colorXIOS.nc -t colorvalues.txt -A 'dataset:Henderson-Sellers!and!Wilson:S' -S 'soil_color' |
---|
11 | ## e.g. # file_rectifier.py -f soils_param.nc -v soiltext -r 0,7 -d x,y -l nav_lon,nav_lat -m maskabove@1.e-8 -o soils_param_Zobler_textXIOS.nc -t soilsparam_text.txt -S soil_texture -A 'article:Zobler!L.;!1986;!NASA:S' |
---|
12 | ## e.g. # file_rectifier.py -f PFTmap_IPCC_2000.nc -v maxvegetfrac -r 0.,1. -d lon,lat -l lon,lat -m masksumrange@0.9999999,1.e-8,2. -o PFTmap_IPCC_2000.XIOS.nc -T None,veget -I time_counter,time_counter -t PFTmap_types.txt -Y time,years@since@0001-01-01@00:00:00,standard -S maxvegetfrac |
---|
13 | ## e.g. # file_rectifier.py -f cartepente2d_15min.nc -v pente -r 0.,10. -d longitude,latitude -l longitude,latitude -m maskmissing -o cartepente2d_15min_XIOS.nc -S slope -A 'database:NGDC!-!etopo2:S' |
---|
14 | ## e.g. # file_rectifier.py -f orchidee_bbg_clim.nc -v flx_co2_bbg_year -r 0.,2. -d lon,lat -l lon,lat -m maskfile@soils_param_Zobler_textXIOS.nc@mask@lon@lat -o orchidee_bbg_clim_XIOS.nc -S flx_co2_bbg_year -z 0. -A 'dataset:D.!Haugustaine!and!G.!van!der!Werf!12-2003:S' |
---|
15 | ## e.g. # file_rectifier.py -f orchidee_fertilizer_1995.nc -v N_qt_WRICE_year -r 0.,6.e+06 -d lon,lat -l lon,lat -m maskabove@0. -o orchidee_fertilizer_1995_N_qt_WRICE_yearXIOS.nc -z 0. -S N_qt_WRICE_year -A 'dataset:A.!BOUWMAN!RIVM!07-2003:S' |
---|
16 | ## e.g. # file_rectifier.py -f orchidee_fertilizer_1995.nc -v N_qt_OTHER_year -r 0.,3.e+07 -d lon,lat -l lon,lat -m maskabove@0. -o orchidee_fertilizer_1995_N_qt_OTHER_yearXIOS.nc -z 0. -S N_qt_OTHER_year -A 'dataset:A.!BOUWMAN!RIVM!07-2003:S' |
---|
17 | ## e.g. # file_rectifier.py -f routing.nc -v topoind -r 1.,9900. -d x,y -l nav_lon,nav_lat -m maskmissing -o routing_topoind_XIOS.nc -S topo_index -A 'associate:nav_lat!na:S' |
---|
18 | ## e.g. # file_rectifier.py -f routing.nc -v trip -r 1.,99. -d x,y -l nav_lon,nav_lat -m maskmissing -o routing_trip_XIOS.nc -S trip -A 'associate:nav_lat!na:S' |
---|
19 | ## e.g. # file_rectifier.py -f routing.nc -v basins -r 1.,6740. -d x,y -l nav_lon,nav_lat -m maskmissing -o routing_basins_XIOS.nc -S trip -A 'associate:nav_lat!na:S' |
---|
20 | # |
---|
21 | |
---|
22 | from optparse import OptionParser |
---|
23 | import numpy as np |
---|
24 | from netCDF4 import Dataset as NetCDFFile |
---|
25 | import os |
---|
26 | import re |
---|
27 | import numpy.ma as ma |
---|
28 | |
---|
29 | main = 'file_rectifier.py' |
---|
30 | version='1.0' |
---|
31 | date='April 2016' |
---|
32 | errormsg = 'ERROR -- error -- ERROR -- error' |
---|
33 | warnmsg = 'WARNING -- warning -- WARNING -- warning' |
---|
34 | |
---|
35 | def searchInlist(listname, nameFind): |
---|
36 | """ Function to search a value within a list |
---|
37 | listname = list |
---|
38 | nameFind = value to find |
---|
39 | >>> searInlist(['1', '2', '3', '5'], '5') |
---|
40 | True |
---|
41 | """ |
---|
42 | for x in listname: |
---|
43 | if x == nameFind: |
---|
44 | return True |
---|
45 | return False |
---|
46 | |
---|
47 | def printing_class(classobj): |
---|
48 | """ Function to print all the values of a given class |
---|
49 | """ |
---|
50 | |
---|
51 | valscls = vars(classobj) |
---|
52 | for attrcls in valscls: |
---|
53 | print attrcls, ':', valscls[attrcls] |
---|
54 | |
---|
55 | def count_cond(var, val, con): |
---|
56 | """ Function to count values of a variable which attain a condition |
---|
57 | [var]= variable |
---|
58 | [val]= value |
---|
59 | [con]= statistics/condition |
---|
60 | 'eq': equal than [val] |
---|
61 | 'neq': not equal than [val] |
---|
62 | 'lt': less than [val] |
---|
63 | 'le': less and equal than [val] |
---|
64 | 'gt': greater than [val] |
---|
65 | 'ge': greater and equal than [val] |
---|
66 | 'in': inside the range [val[0]] <= [var] <= [val[1]] (where [val]=[val1, val2] |
---|
67 | 'out': inside the range [val[0]] > [var]; [val[1]] < [var] (where [val]=[val1, val2] |
---|
68 | >>> vals = np.arange(100) |
---|
69 | >>> count_cond(vals,50,'eq') |
---|
70 | 1 |
---|
71 | >>> count_cond(vals,50,'neq') |
---|
72 | 99 |
---|
73 | >>> count_cond(vals,50,'lt') |
---|
74 | 50 |
---|
75 | >>> count_cond(vals,50,'le') |
---|
76 | 51 |
---|
77 | >>> count_cond(vals,50,'gt') |
---|
78 | 49 |
---|
79 | >>> count_cond(vals,50,'ge') |
---|
80 | 50 |
---|
81 | >>> count_cond(vals,[25,75],'in') |
---|
82 | 51 |
---|
83 | >>> count_cond(vals,[25,75],'out') |
---|
84 | 49 |
---|
85 | """ |
---|
86 | fname = 'count_cond' |
---|
87 | |
---|
88 | cons = ['eq', 'neq', 'lt', 'le', 'gt', 'ge', 'in', 'out'] |
---|
89 | |
---|
90 | if con == 'eq': |
---|
91 | Nvals = np.sum(var == val) |
---|
92 | elif con == 'neq': |
---|
93 | Nvals = np.sum(var != val) |
---|
94 | elif con == 'lt': |
---|
95 | Nvals = np.sum(var < val) |
---|
96 | elif con == 'le': |
---|
97 | Nvals = np.sum(var <= val) |
---|
98 | elif con == 'gt': |
---|
99 | Nvals = np.sum(var > val) |
---|
100 | elif con == 'ge': |
---|
101 | Nvals = np.sum(var >= val) |
---|
102 | elif con == 'in': |
---|
103 | if len(val) != 2: |
---|
104 | print errormsg |
---|
105 | print ' ' + fname + ": for condition '" + con + "' we must have two " + \ |
---|
106 | " values!!" |
---|
107 | print ' we have:', val |
---|
108 | quit(-1) |
---|
109 | Nvals = np.sum((var >= val[0])*(var <= val[1])) |
---|
110 | elif con == 'out': |
---|
111 | if len(val) != 2: |
---|
112 | print errormsg |
---|
113 | print ' ' + fname + ": for condition '" + con + "' we must have two " + \ |
---|
114 | " values!!" |
---|
115 | print ' we have:', val |
---|
116 | quit(-1) |
---|
117 | Nvals = np.sum((var < val[0])+(var > val[1])) |
---|
118 | else: |
---|
119 | print errormsg |
---|
120 | print ' ' + fname + ": condition '" + con + "' not ready!!" |
---|
121 | print ' only ready:', cons |
---|
122 | quit(-1) |
---|
123 | |
---|
124 | return Nvals |
---|
125 | |
---|
126 | def index_mat(mat,val): |
---|
127 | """ Function to provide the coordinates of a given value inside a matrix |
---|
128 | index_mat(mat,val) |
---|
129 | mat= matrix with values |
---|
130 | val= value to search |
---|
131 | >>> index_mat(np.arange(27).reshape(3,3,3),22) |
---|
132 | [2 1 1] |
---|
133 | """ |
---|
134 | |
---|
135 | fname = 'index_mat' |
---|
136 | |
---|
137 | matshape = mat.shape |
---|
138 | |
---|
139 | matlist = list(mat.flatten()) |
---|
140 | ifound = matlist.index(val) |
---|
141 | |
---|
142 | Ndims = len(matshape) |
---|
143 | valpos = np.zeros((Ndims), dtype=int) |
---|
144 | baseprevdims = np.zeros((Ndims), dtype=int) |
---|
145 | |
---|
146 | for dimid in range(Ndims): |
---|
147 | baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims]) |
---|
148 | if dimid == 0: |
---|
149 | alreadyplaced = 0 |
---|
150 | else: |
---|
151 | alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid]) |
---|
152 | valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid]) |
---|
153 | |
---|
154 | return valpos |
---|
155 | |
---|
156 | def multi_index_mat(mat,val): |
---|
157 | """ Function to provide the multiple coordinates of a given value inside a matrix |
---|
158 | index_mat(mat,val) |
---|
159 | mat= matrix with values |
---|
160 | val= value to search |
---|
161 | >>> vals = np.ones((24), dtype=np.float).reshape(2,3,4) |
---|
162 | vals[:,:,2] = 0. |
---|
163 | vals[1,:,:] = np.pi |
---|
164 | vals[:,2,:] = -1. |
---|
165 | multi_index_mat(vals,1.) |
---|
166 | [array([0, 0, 0]), array([0, 0, 1]), array([0, 0, 3]), array([0, 1, 0]), array([0, 1, 1]), array([0, 1, 3])] |
---|
167 | """ |
---|
168 | fname = 'multi_index_mat' |
---|
169 | |
---|
170 | matshape = mat.shape |
---|
171 | |
---|
172 | ivalpos = [] |
---|
173 | matlist = list(mat.flatten()) |
---|
174 | Lmatlist = len(matlist) |
---|
175 | |
---|
176 | val0 = val - val |
---|
177 | if val != val0: |
---|
178 | valdiff = val0 |
---|
179 | else: |
---|
180 | valdiff = np.ones((1), dtype = type(val)) |
---|
181 | |
---|
182 | ifound = 0 |
---|
183 | while ifound < Lmatlist: |
---|
184 | if matlist.count(val) == 0: |
---|
185 | ifound = Lmatlist + 1 |
---|
186 | else: |
---|
187 | ifound = matlist.index(val) |
---|
188 | |
---|
189 | Ndims = len(matshape) |
---|
190 | valpos = np.zeros((Ndims), dtype=int) |
---|
191 | baseprevdims = np.zeros((Ndims), dtype=int) |
---|
192 | |
---|
193 | for dimid in range(Ndims): |
---|
194 | baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims]) |
---|
195 | if dimid == 0: |
---|
196 | alreadyplaced = 0 |
---|
197 | else: |
---|
198 | alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid]) |
---|
199 | valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid]) |
---|
200 | matlist[ifound] = valdiff |
---|
201 | ivalpos.append(valpos) |
---|
202 | |
---|
203 | return ivalpos |
---|
204 | |
---|
205 | class variable_inf(object): |
---|
206 | """ Class to provide the information from a given variable |
---|
207 | var = object netCDF variable |
---|
208 | self.name: name of the variable |
---|
209 | self.dtype: type of the variable |
---|
210 | self.attributes: list with the name of attributes |
---|
211 | self.FillValue: value of the missing value |
---|
212 | self.dimns: name of the dimensions of the variable |
---|
213 | self.dims: dimensions of the variable |
---|
214 | self.Ndims: number of dimensions |
---|
215 | self.dimx: length of dimension along x-axis |
---|
216 | self.dimy: length of dimension along y-axis |
---|
217 | self.sname: standard name |
---|
218 | self.lname: long name |
---|
219 | self.corr: attribute 'coordinates' |
---|
220 | self.units: units of the variable |
---|
221 | """ |
---|
222 | |
---|
223 | def __init__(self, var): |
---|
224 | |
---|
225 | if var is None: |
---|
226 | self.name = None |
---|
227 | self.dimns = None |
---|
228 | self.dims = None |
---|
229 | self.Ndims = None |
---|
230 | self.dimx = None |
---|
231 | self.dimy = None |
---|
232 | self.sname = None |
---|
233 | self.lname = None |
---|
234 | self.corr = None |
---|
235 | self.units = None |
---|
236 | self.FillValue = None |
---|
237 | self.dtype = None |
---|
238 | self.attributes = None |
---|
239 | else: |
---|
240 | self.name = var._name |
---|
241 | self.dimns = var.dimensions |
---|
242 | self.dtype = var.dtype |
---|
243 | self.attributes = var.ncattrs() |
---|
244 | self.dims = var.shape |
---|
245 | |
---|
246 | if searchInlist(self.attributes, 'standard_name'): |
---|
247 | self.sname = var.getncattr('standard_name') |
---|
248 | else: |
---|
249 | # print ' variable_inf.classpy: variable "' + self.name + '" does not have attribute "standard_name"' |
---|
250 | self.sname = None |
---|
251 | |
---|
252 | if searchInlist(self.attributes, 'long_name'): |
---|
253 | self.lname = var.getncattr('long_name') |
---|
254 | else: |
---|
255 | # print ' variable_inf.classpy: variable "' + self.name + '" does not have attribute "long_name"' |
---|
256 | self.lname = None |
---|
257 | |
---|
258 | if searchInlist(self.attributes, 'coordinates'): |
---|
259 | self.coor = var.getncattr('coordinates') |
---|
260 | else: |
---|
261 | # print ' variable_inf.classpy: variable "' + self.name + '" does not have attribute "coordinates"' |
---|
262 | self.coor = None |
---|
263 | |
---|
264 | if searchInlist(self.attributes, 'units'): |
---|
265 | self.units = var.getncattr('units') |
---|
266 | else: |
---|
267 | # print ' variable_inf.classpy: variable "' + self.name + '" does not have attribute "units"' |
---|
268 | self.units = None |
---|
269 | |
---|
270 | if searchInlist(self.attributes, '_FillValue'): |
---|
271 | self.FillValue = var.getncattr('_FillValue') |
---|
272 | else: |
---|
273 | # print ' variable_inf.classpy: variable "' + self.name + '" does not have attribute "_FillValue"' |
---|
274 | self.FillValue = None |
---|
275 | |
---|
276 | self.Ndims = len(self.dims) |
---|
277 | if self.Ndims == 1: |
---|
278 | self.dimx=self.dims[0] |
---|
279 | if self.Ndims == 2: |
---|
280 | self.dimy=self.dims[0] |
---|
281 | self.dimx=self.dims[1] |
---|
282 | if self.Ndims == 3: |
---|
283 | self.dimy=self.dims[1] |
---|
284 | self.dimx=self.dims[2] |
---|
285 | if self.Ndims == 4: |
---|
286 | self.dimy=self.dims[2] |
---|
287 | self.dimx=self.dims[3] |
---|
288 | |
---|
289 | |
---|
290 | def lonlat2D(ncfobj, lonn, latn): |
---|
291 | """ Function to create 2D matricies of longitudes and latitudes from the file |
---|
292 | ncfobj= netCDF object file |
---|
293 | lonn= longitude name |
---|
294 | latn= latitude name |
---|
295 | """ |
---|
296 | fname='lonlat2D' |
---|
297 | |
---|
298 | lonvar = ncfobj.variables[lonn] |
---|
299 | latvar = ncfobj.variables[latn] |
---|
300 | |
---|
301 | loninf = variable_inf(lonvar) |
---|
302 | latinf = variable_inf(latvar) |
---|
303 | |
---|
304 | if loninf.Ndims == 1: |
---|
305 | lonv = lonvar[:] |
---|
306 | latv = latvar[:] |
---|
307 | dx=loninf.dims[0] |
---|
308 | dy=latinf.dims[0] |
---|
309 | lonval = np.zeros((dy, dx), dtype=float) |
---|
310 | latval = np.zeros((dy, dx), dtype=float) |
---|
311 | for iy in range(dy): |
---|
312 | lonval[iy,:] = lonv |
---|
313 | for ix in range(dx): |
---|
314 | latval[:,ix] = latv |
---|
315 | elif loninf.Ndims == 2: |
---|
316 | lonval = lonvar[:] |
---|
317 | latval = latvar[:] |
---|
318 | elif loninf.Ndims == 3: |
---|
319 | lonval = lonvar[0,:,:] |
---|
320 | latval = latvar[0,:,:] |
---|
321 | else: |
---|
322 | print errormsg |
---|
323 | print ' ' + fname + ' dimension ', lonvar.Ndims, ' of the lon/lat matrices not ready !!!!' |
---|
324 | |
---|
325 | return lonval, latval |
---|
326 | |
---|
327 | def varval_ind(var, ind): |
---|
328 | """ Function to provide a value from a variable providing a value for each index |
---|
329 | var= variable |
---|
330 | ind= list of indices for each dimension of the variable |
---|
331 | """ |
---|
332 | fname = 'varval_ind' |
---|
333 | if len(var.shape) != len(ind): |
---|
334 | print errormsg |
---|
335 | print ' ' + fname + ': different shape:', var.shape,'and indices:', ind, '!!' |
---|
336 | quit(-1) |
---|
337 | |
---|
338 | varslice=[] |
---|
339 | for iv in ind: |
---|
340 | varslice.append(iv) |
---|
341 | |
---|
342 | val = varval_ind[slice(varslice)] |
---|
343 | |
---|
344 | return val |
---|
345 | |
---|
346 | def retype(val, vtype): |
---|
347 | """ Function to transform a value to a given type |
---|
348 | retype(val, vtype) |
---|
349 | [val]= value |
---|
350 | [vtype]= type to transform |
---|
351 | >>> retype(0, type(True)) |
---|
352 | False |
---|
353 | """ |
---|
354 | |
---|
355 | fname = 'retype' |
---|
356 | |
---|
357 | if val == 'h': |
---|
358 | print fname + '_____________________________________________________________' |
---|
359 | print retype.__doc__ |
---|
360 | quit() |
---|
361 | |
---|
362 | if vtype == type(int(1)): |
---|
363 | newval = int(val) |
---|
364 | elif vtype == type(float(1)): |
---|
365 | newval = float(val) |
---|
366 | # elif vtype == type(float32(1)): |
---|
367 | # newval = float32(val) |
---|
368 | elif vtype == type(np.int(1)): |
---|
369 | typeinf = np.iinfo(np.int) |
---|
370 | if (abs(np.float64(val)) > typeinf.max): |
---|
371 | print warnmsg |
---|
372 | print ' ' + fname + ': value to transform ', val,' overpasses type:', \ |
---|
373 | vtype,' limits!!' |
---|
374 | print ' Changing value to kind largest allowed value:', typeinf.max |
---|
375 | newval = abs(np.float64(val))/np.float64(val) * np.int(typeinf.max) |
---|
376 | else: |
---|
377 | newval = np.int(val) |
---|
378 | elif vtype == type(np.int16(1)): |
---|
379 | typeinf = np.iinfo(np.int16) |
---|
380 | if (abs(np.float64(val)) > typeinf.max): |
---|
381 | print warnmsg |
---|
382 | print ' ' + fname + ': value to transform ', val,' overpasses type:', \ |
---|
383 | vtype,' limits!!' |
---|
384 | print ' Changing value to kind largest allowed value:', typeinf.max |
---|
385 | newval = abs(np.float64(val))/np.float64(val) * np.int16(typeinf.max) |
---|
386 | else: |
---|
387 | newval = np.int16(val) |
---|
388 | elif vtype == type(np.int32(1)): |
---|
389 | typeinf = np.iinfo(np.int32) |
---|
390 | if (abs(np.float64(val)) > typeinf.max): |
---|
391 | print warnmsg |
---|
392 | print ' ' + fname + ': value to transform ', val,' overpasses type:', \ |
---|
393 | vtype,' limits!!' |
---|
394 | print ' Changing value to kind largest allowed value:', typeinf.max |
---|
395 | newval = abs(np.float64(val))/np.float64(val) * np.int32(typeinf.max) |
---|
396 | else: |
---|
397 | newval = np.int32(val) |
---|
398 | elif vtype == type(np.int64(1)): |
---|
399 | newval = np.int64(val) |
---|
400 | elif vtype == type(np.float(1)): |
---|
401 | newval = np.float(val) |
---|
402 | # elif vtype == type(np.float16(1)): |
---|
403 | # newval = np.float16(val) |
---|
404 | elif vtype == type(np.float32(1)): |
---|
405 | newval = np.float32(val) |
---|
406 | elif vtype == type(np.float64(1)): |
---|
407 | newval = np.float64(val) |
---|
408 | elif vtype == type(True): |
---|
409 | if val == 0: |
---|
410 | newval = False |
---|
411 | else: |
---|
412 | newval = True |
---|
413 | elif vtype == '|S1': |
---|
414 | newval = str(val) |
---|
415 | else: |
---|
416 | print errormsg |
---|
417 | print ' ' + fname + ': variable type "', vtype, '" not ready!!!' |
---|
418 | quit(-1) |
---|
419 | |
---|
420 | return newval |
---|
421 | |
---|
422 | def fillvalue_kind(vartype, fillval0): |
---|
423 | """ Function to provide the fiven fill_Value for a kind of variable |
---|
424 | [vartype]= type of the variable |
---|
425 | [fillval0]= value to take as fill value, 'std' for the standard value |
---|
426 | Float = 1.e20 |
---|
427 | Character = '-' |
---|
428 | Integer = -99999 |
---|
429 | Integer16 = -99999 |
---|
430 | Float64 = 1.e20 |
---|
431 | Integer32 = -99999 |
---|
432 | """ |
---|
433 | fname = 'fillvalue_kind' |
---|
434 | fillValue = 1.e20 |
---|
435 | fillValueF = 1.e20 |
---|
436 | fillValueC = '-' |
---|
437 | fillValueI = -99999 |
---|
438 | fillValueI16 = -99999 |
---|
439 | fillValueF64 = 1.e20 |
---|
440 | fillValueI32 = -99999 |
---|
441 | |
---|
442 | if vartype == type('c'): |
---|
443 | if fillval0 == 'std': |
---|
444 | fillval = fillvalueC |
---|
445 | else: |
---|
446 | fillval = fillval0 |
---|
447 | elif vartype == type(int(1)): |
---|
448 | if fillval0 == 'std': |
---|
449 | fillval = fillValueI |
---|
450 | else: |
---|
451 | fillval = int(fillval0) |
---|
452 | elif vartype == type(np.int16(1)): |
---|
453 | if fillval0 == 'std': |
---|
454 | fillval = np.int(9999999) |
---|
455 | else: |
---|
456 | fillval = np.int16(fillval0) |
---|
457 | elif vartype == type(np.int32(1)): |
---|
458 | if fillval0 == 'std': |
---|
459 | fillval = np.int32(fillValueI) |
---|
460 | else: |
---|
461 | fillval = np.int32(fillval0) |
---|
462 | elif vartype == type(np.float(1.)): |
---|
463 | if fillval0 == 'std': |
---|
464 | fillval = np.float(fillValueF) |
---|
465 | else: |
---|
466 | fillval = np.float(fillval0) |
---|
467 | elif vartype == type(np.float32(1.)): |
---|
468 | if fillval0 == 'std': |
---|
469 | fillval = np.float32(fillValueF) |
---|
470 | else: |
---|
471 | fillval = np.float32(fillval0) |
---|
472 | elif vartype == type(np.float64(1.)): |
---|
473 | if fillval0 == 'std': |
---|
474 | fillval = np.float64(fillValueF) |
---|
475 | else: |
---|
476 | fillval = np.float64(fillval0) |
---|
477 | else: |
---|
478 | print errormsg |
---|
479 | print ' ' + fname + ': variable type', vartype, 'not ready !!' |
---|
480 | quit(-1) |
---|
481 | |
---|
482 | return fillval |
---|
483 | |
---|
484 | def read_ASCIIlist(filen): |
---|
485 | """ Function to build a list from an ASCII lines file |
---|
486 | filen= name of the ASCII file |
---|
487 | """ |
---|
488 | import os |
---|
489 | fname = 'read_ASCIIlist' |
---|
490 | |
---|
491 | if not os.path.isfile(filen): |
---|
492 | print errormsg |
---|
493 | print ' ' + fname + ": ASCII types file '" + filen + "' does not exist !!" |
---|
494 | quit(-1) |
---|
495 | |
---|
496 | of = open(filen, 'r') |
---|
497 | filevals = [] |
---|
498 | for linef in of: |
---|
499 | filevals.append(linef.replace('\n','')) |
---|
500 | of.close() |
---|
501 | |
---|
502 | return filevals |
---|
503 | |
---|
504 | def set_attributek(ncv, attrname, attrval, attrkind): |
---|
505 | """ Sets a value of an attribute of a netCDF variable with a kind. Removes previous attribute value if exists |
---|
506 | ncvar = object netcdf variable |
---|
507 | attrname = name of the attribute |
---|
508 | attrvalue = value of the attribute |
---|
509 | atrtrkind = kind of attribute: 'S', string ('!' as spaces); 'U', unicode ('!' as spaces); 'I', integer; |
---|
510 | 'Inp32', numpy integer 32; 'R', real; ' npfloat', np.float; 'D', np.float64 |
---|
511 | """ |
---|
512 | fname = 'set_attributek' |
---|
513 | validk = {'S': 'string', 'U': 'unicode', 'I': 'integer', \ |
---|
514 | 'Inp32': 'integer long (np.int32)', 'R': 'float', 'npfloat': 'np.float', \ |
---|
515 | 'D': 'float long (np.float64)'} |
---|
516 | |
---|
517 | if type(attrkind) == type('s'): |
---|
518 | if attrkind == 'S': |
---|
519 | attrvalue = str(attrval.replace('!', ' ')) |
---|
520 | elif attrkind == 'U': |
---|
521 | attrvalue = unicode(attrval.replace('!',' ')) |
---|
522 | elif attrkind == 'I': |
---|
523 | attrvalue = np.int(attrval) |
---|
524 | elif attrkind == 'R': |
---|
525 | attrvalue = float(attrval) |
---|
526 | elif attrkind == 'npfloat': |
---|
527 | attrvalue = np.float(attrval) |
---|
528 | elif attrkind == 'D': |
---|
529 | attrvalue = np.float64(attrval) |
---|
530 | else: |
---|
531 | print errormsg |
---|
532 | print ' ' + fname + ": '" + attrkind + "' kind of attribute is not ready!" |
---|
533 | print ' valid values: _______' |
---|
534 | for key in validk.keys(): |
---|
535 | print ' ', key,':', validk[key] |
---|
536 | quit(-1) |
---|
537 | else: |
---|
538 | if attrkind == type(str('a')): |
---|
539 | attrvalue = str(attrval.replace('!', ' ')) |
---|
540 | elif attrkind == type(unicode('a')): |
---|
541 | attrvalue = unicode(attrval.replace('!',' ')) |
---|
542 | elif attrkind == type(np.int(1)): |
---|
543 | attrvalue = np.int(attrval) |
---|
544 | elif attrkind == np.dtype('i'): |
---|
545 | attrvalue = np.int32(attrval) |
---|
546 | elif attrkind == type(float(1.)): |
---|
547 | attrvalue = float(attrval) |
---|
548 | elif attrkind == type(np.float(1.)): |
---|
549 | attrvalue = np.float(attrval) |
---|
550 | elif attrkind == np.dtype('float32'): |
---|
551 | attrvalue = np.array(attrval, dtype='f') |
---|
552 | elif attrkind == type(np.float32(1.)): |
---|
553 | attrvalue = np.float32(attrval) |
---|
554 | elif attrkind == type(np.float64(1.)): |
---|
555 | attrvalue = np.float64(attrval) |
---|
556 | else: |
---|
557 | print errormsg |
---|
558 | print ' ' + fname + ": '" + attrkind + "' kind of attribute is not ready!" |
---|
559 | print ' valid values: _______' |
---|
560 | for key in validk.keys(): |
---|
561 | print ' ', key,':', validk[key] |
---|
562 | quit(-1) |
---|
563 | |
---|
564 | attvar = ncv.ncattrs() |
---|
565 | if searchInlist(attvar, attrname): |
---|
566 | attr = ncv.delncattr(attrname) |
---|
567 | |
---|
568 | attr = ncv.setncattr(attrname, attrvalue) |
---|
569 | |
---|
570 | return attr |
---|
571 | |
---|
572 | def nctype(vartype): |
---|
573 | """ Function to provide the string for the variable creation in a netCDF file |
---|
574 | [vartype]= type of the variable |
---|
575 | """ |
---|
576 | fname = 'nctype' |
---|
577 | |
---|
578 | if vartype == type('c'): |
---|
579 | ncs = 'c' |
---|
580 | elif vartype == type(int(1)): |
---|
581 | ncs = 'i' |
---|
582 | elif vartype == type(np.int16(1)): |
---|
583 | ncs = 'i' |
---|
584 | elif vartype == type(np.int32(1)): |
---|
585 | ncs = 'i8' |
---|
586 | elif vartype == type(np.float(1.)): |
---|
587 | ncs = 'f' |
---|
588 | elif vartype == type(np.float32(1.)): |
---|
589 | ncs = 'f4' |
---|
590 | elif vartype == type(np.float64(1.)): |
---|
591 | ncs = 'f8' |
---|
592 | else: |
---|
593 | print errormsg |
---|
594 | print ' ' + fname + ': variable type', vartype, 'not ready !!' |
---|
595 | quit(-1) |
---|
596 | |
---|
597 | return ncs |
---|
598 | |
---|
599 | def writing_str_nc(varo, values, Lchar): |
---|
600 | """ Function to write string values in a netCDF variable as a chain of 1char values |
---|
601 | varo= netCDF variable object |
---|
602 | values = list of values to introduce |
---|
603 | Lchar = length of the string in the netCDF file |
---|
604 | """ |
---|
605 | fname = 'writing_str_nc' |
---|
606 | |
---|
607 | Nvals = len(values) |
---|
608 | for iv in range(Nvals): |
---|
609 | if values[iv] is None: |
---|
610 | stringv = 'None' |
---|
611 | else: |
---|
612 | stringv = values[iv] |
---|
613 | charvals = np.zeros((Lchar), dtype='string_') |
---|
614 | Lstr = len(stringv) |
---|
615 | charvals[Lstr:Lchar] = '' |
---|
616 | |
---|
617 | for ich in range(Lstr): |
---|
618 | chrv = stringv[ich:ich+1] |
---|
619 | charvals[ich] = stringv[ich:ich+1] |
---|
620 | varo[iv,:] = charvals |
---|
621 | |
---|
622 | return |
---|
623 | |
---|
624 | ####### ###### ##### #### ### ## # |
---|
625 | maskhelp="(option to gnerate the mask: [maskoption]@[values]\ |
---|
626 | 'maskfile'@[filename]@[maskname]@[lonname]@[latname]: take mask values from another file [filename] with mask in variable [maskname] and [lonname], [latname] lon lat variables;\ |
---|
627 | 'maskvarname'@[varname]: reference file contents a variable [varname] with its current mask;\ |
---|
628 | 'maskbelow'@[value]: mask below a given [value];\ |
---|
629 | 'maskabove'@[value]: mask above a given [value];\ |
---|
630 | 'maskrange'@[value1],[value2]: mask within a given range [value1] <= [value] <= [value2];\ |
---|
631 | 'masksumrange'@[value1],[value2],pvalue3] masking taking the sum of all the types at each grid point (SUMij) which [value2] <= SUMij <= [value1] and [value1] < SUMij <= [value3]);\ |
---|
632 | 'maskmissing': use missing value for the mask" |
---|
633 | extrattrhelp="',' list of extra attributes to be added to the variable as [attrn]:[value]:[attrk] ([attrk], attribute kind: 'S', string ('!' as spaces); 'U', unicode ('!' as spaces); 'I', integer; 'Inp32', numpy integer 32; 'R', real; ' npfloat', np.float; 'D', np.float64)" |
---|
634 | |
---|
635 | parser = OptionParser() |
---|
636 | parser.add_option("-b", "--lonlat_bounds", dest="Llbounds", |
---|
637 | help="[OPTIONAL] when is set and it is given a 'True' value, lon/lat bounds are computed", metavar="FILE") |
---|
638 | parser.add_option("-d", "--lonlatdimname", dest="lonlatdimn", |
---|
639 | help="current names of the longitude and latitude dimensions as londimn,latdimn in the reference file", metavar="LABEL") |
---|
640 | parser.add_option("-f", "--netCDF_file", dest="ncfile", |
---|
641 | help="file to use as reference", metavar="FILE") |
---|
642 | parser.add_option("-l", "--lonlatname", dest="lonlatn", |
---|
643 | help="current names of the longitude and latitude variables as lonn,latn in the reference file", metavar="LABEL") |
---|
644 | parser.add_option("-m", "--maskname", dest="maskn", |
---|
645 | help=maskhelp, metavar="LABEL") |
---|
646 | parser.add_option("-o", "--outputfile", dest="ofile", |
---|
647 | help="output file", metavar="FILE") |
---|
648 | parser.add_option("-r", "--range", dest="range", |
---|
649 | help="range of the variable as varmin,varmax", metavar="VALUE") |
---|
650 | parser.add_option("-A", "--ExtraAttributes", dest="extrattrs", help=extrattrhelp, metavar="LABEL") |
---|
651 | parser.add_option("-T", "--typesvarn", dest="typen", |
---|
652 | help="[IF APPLICABLE] current name of the type variable and dimension in the reference file (as [typen],[typedimn], " + |
---|
653 | "'None' for no variable)", metavar="LABEL") |
---|
654 | parser.add_option("-I", "--timevarn", dest="timen", |
---|
655 | help="[IF APPLICABLE] current name of the time variable and dimension in the reference file (as [timen],[timedimn])", |
---|
656 | metavar="LABEL") |
---|
657 | parser.add_option("-t", "--types_file", dest="typesfile", |
---|
658 | help="[OPTIONAL] ASCII file with the values (one consecutive line per value, @ for spaces) for the description of " + |
---|
659 | "the types of the variable", metavar="FILE") |
---|
660 | parser.add_option("-v", "--variable", dest="varname", |
---|
661 | help="variable to check", metavar="LABEL") |
---|
662 | parser.add_option("-z", "--ocean_value", dest="oceanval", |
---|
663 | help="[OPTIONAL] current ocean value to assign as fill_value of the variable", metavar="VALUE") |
---|
664 | parser.add_option("-S", "--varstd_name", dest="svname", |
---|
665 | help="[OPTIONAL] label as the standard_name for the variable", metavar="LABEL") |
---|
666 | parser.add_option("-L", "--varlong_name", dest="lvname", |
---|
667 | help="[OPTIONAL] label as the long_name for the variable (@ for spaces)", metavar="LABEL") |
---|
668 | parser.add_option("-U", "--varunits", dest="vunits", |
---|
669 | help="[OPTIONAL] label as units for the variable", metavar="LABEL") |
---|
670 | parser.add_option("-Y", "--time_attributes", dest="timeattrs", |
---|
671 | help="[OPTIONAL] long name and units of the time variable (as [long_name],[units],[calendar] @ for spaces)", metavar="LABEL") |
---|
672 | |
---|
673 | (opts, args) = parser.parse_args() |
---|
674 | |
---|
675 | ####### ####### |
---|
676 | ## MAIN |
---|
677 | ####### |
---|
678 | ofile = opts.ofile |
---|
679 | |
---|
680 | if not os.path.isfile(opts.ncfile): |
---|
681 | print errormsg |
---|
682 | print ' ' + main + ": File '" + opts.ncfile + "' does not exist !!" |
---|
683 | quit(-1) |
---|
684 | |
---|
685 | ncobj = NetCDFFile(opts.ncfile, 'r') |
---|
686 | newncobj = NetCDFFile(opts.ofile, 'w') |
---|
687 | |
---|
688 | if opts.typen is not None: |
---|
689 | typen = opts.typen.split(',')[0] |
---|
690 | typedimn = opts.typen.split(',')[1] |
---|
691 | |
---|
692 | if opts.timen is not None: |
---|
693 | timen = opts.timen.split(',')[0] |
---|
694 | timedimn = opts.timen.split(',')[1] |
---|
695 | |
---|
696 | londimn = opts.lonlatdimn.split(',')[0] |
---|
697 | latdimn = opts.lonlatdimn.split(',')[1] |
---|
698 | # Checking presence of mandatory variables |
---|
699 | ## |
---|
700 | lonn = opts.lonlatn.split(',')[0] |
---|
701 | latn = opts.lonlatn.split(',')[1] |
---|
702 | varcheck = [lonn, latn, opts.varname] |
---|
703 | |
---|
704 | for varchk in varcheck: |
---|
705 | if not searchInlist(ncobj.variables, varchk): |
---|
706 | print errormsg |
---|
707 | print ' ' + main + ": File '" + opts.ncfile +"' does not have variable '" + \ |
---|
708 | varchk + "' !!" |
---|
709 | quit(-1) |
---|
710 | |
---|
711 | # Longitude |
---|
712 | olon = ncobj.variables[lonn] |
---|
713 | loninf = variable_inf(olon) |
---|
714 | |
---|
715 | if len(olon.shape) > 1: |
---|
716 | newdim = newncobj.createDimension('lon',loninf.dims[1]) |
---|
717 | newdim = newncobj.createDimension('lat',loninf.dims[0]) |
---|
718 | vardims = ('lat', 'lon') |
---|
719 | else: |
---|
720 | newdim = newncobj.createDimension('lon',loninf.dims[0]) |
---|
721 | vardims = ('lon') |
---|
722 | |
---|
723 | print " Creation of variable 'longitude' ..." |
---|
724 | newlon = newncobj.createVariable('lon', 'f8', vardims) |
---|
725 | newattr = set_attributek(newlon, 'standard_name', 'longitude', 'S') |
---|
726 | newattr = set_attributek(newlon, 'long_name', 'Longitude', 'S') |
---|
727 | newattr = set_attributek(newlon, 'units', 'degrees_East', 'S') |
---|
728 | newattr = set_attributek(newlon, 'axis', 'X', 'S') |
---|
729 | newattr = set_attributek(newlon, '_CoordinateAxisType', 'Lon', 'S') |
---|
730 | newlon[:] = olon[:] |
---|
731 | if opts.Llbounds is not None and opts.Llbounds == 'True': |
---|
732 | newattr = set_attributek(newlon, 'bounds', 'lon_bnds', 'S') |
---|
733 | |
---|
734 | # Latitude |
---|
735 | olat = ncobj.variables[latn] |
---|
736 | latinf = variable_inf(olat) |
---|
737 | |
---|
738 | if len(olat.shape) > 1: |
---|
739 | vardims = ('lat', 'lon') |
---|
740 | else: |
---|
741 | newdim = newncobj.createDimension('lat',latinf.dims[0]) |
---|
742 | vardims = ('lat') |
---|
743 | |
---|
744 | print " Creation of variable 'latitude' ..." |
---|
745 | newlat = newncobj.createVariable('lat', 'f8', vardims) |
---|
746 | newattr = set_attributek(newlat, 'standard_name', 'latitude', 'S') |
---|
747 | newattr = set_attributek(newlat, 'long_name', 'Latitude', 'S') |
---|
748 | newattr = set_attributek(newlat, 'units', 'degrees_North', 'S') |
---|
749 | newattr = set_attributek(newlat, 'axis', 'Y', 'S') |
---|
750 | newattr = set_attributek(newlat, '_CoordinateAxisType', 'Lat', 'S') |
---|
751 | newlat[:] = olat[:] |
---|
752 | if opts.Llbounds is not None and opts.Llbounds == 'True': |
---|
753 | newattr = set_attributek(newlat, 'bounds', 'lat_bnds', 'S') |
---|
754 | |
---|
755 | newncobj.sync() |
---|
756 | if searchInlist(ncobj.variables, 'lon_bnds') or searchInlist(ncobj.variables, \ |
---|
757 | 'bounds_lon'): |
---|
758 | print ' ' + main + ": file contens 'longitude/latitude bounds'!!" |
---|
759 | print ' writting them in output file' |
---|
760 | if searchInlist(ncobj.variables, 'lon_bnds'): |
---|
761 | olbnd = ncobj.variables['lon_bnds'] |
---|
762 | oLbnd = ncobj.variables['lat_bnds'] |
---|
763 | else: |
---|
764 | olbnd = ncobj.variables['bounds_lon'] |
---|
765 | oLbnd = ncobj.variables['bounds_lat'] |
---|
766 | |
---|
767 | lbnd = olbnd[:] |
---|
768 | Lbnd = oLbnd[:] |
---|
769 | nv = len(olbnd.shape) |
---|
770 | |
---|
771 | newdim = newncobj.createDimension('nv',nv) |
---|
772 | newvar = newncobj.createVariable('lon_bnds', 'f8', ('lon','nv')) |
---|
773 | newvar = lbnd[:] |
---|
774 | newvar = newncobj.createVariable('lat_bnds', 'f8', ('lat','nv')) |
---|
775 | newvar = Lbnd[:] |
---|
776 | newattr = set_attributek(newlon, 'bounds', 'lon_bnds', 'S') |
---|
777 | newattr = set_attributek(newlat, 'bounds', 'lat_bnds', 'S') |
---|
778 | |
---|
779 | |
---|
780 | if opts.Llbounds is not None and opts.Llbounds == 'True': |
---|
781 | # Only working for 'regular' projections! lon(i) < lon(i+1); lat(j) < lat(j+1) |
---|
782 | # Bounds |
---|
783 | print " Creation of variable 'lon_bnds', 'lat_bnds' ..." |
---|
784 | if len(olat.shape) == 1: |
---|
785 | nv = 2 |
---|
786 | newdim = newncobj.createDimension('nv',nv) |
---|
787 | lon_bounds = np.zeros((loninf.dims[0],nv), dtype=np.float32) |
---|
788 | for i in range(loninf.dims[0]-1): |
---|
789 | dlon = (olon[i] - olon[i+1])*0.5 |
---|
790 | lon_bounds[i,0] = olon[i] - dlon |
---|
791 | lon_bounds[i,1] = olon[i] + dlon |
---|
792 | i = loninf.dims[0] |
---|
793 | dlon = olon[i-1] - olon[i] |
---|
794 | lon_bounds[i,0] = olon[i] - dlon |
---|
795 | lon_bounds[i,1] = olon[i] + dlon |
---|
796 | newvar = newncobj.createVariable('lon_bnds', 'f8', ('lon','nv')) |
---|
797 | newvar[:] = lon_bounds[:] |
---|
798 | |
---|
799 | lat_bounds = np.zeros((latinf.dims[0],nv), dtype=np.float32) |
---|
800 | for i in range(latinf.dims[0]-1): |
---|
801 | dlat = (olat[i] - olat[i+1])*0.5 |
---|
802 | lat_bounds[i,0] = olat[i] - dlat |
---|
803 | lat_bounds[i,1] = olat[i] + dlat |
---|
804 | i = latinf.dims[0] |
---|
805 | dlat = olat[i-1] - olat[i] |
---|
806 | lat_bounds[i,0] = olat[i] - dlat |
---|
807 | lat_bounds[i,1] = olat[i] + dlat |
---|
808 | newvar = newncobj.createVariable('lat_bnds', 'f8', ('lat','nv')) |
---|
809 | newvar[:] = lat_bounds[:] |
---|
810 | elif len(olat.shape) == 2: |
---|
811 | # At this stage assuming contiguos projections |
---|
812 | nv = 4 |
---|
813 | newdim = newncobj.createDimension('nv',nv) |
---|
814 | lon_bounds = np.zeros((loninf.dims[0],loninf.dims[1],nv), dtype=np.float32) |
---|
815 | lat_bounds = np.zeros((loninf.dims[0],loninf.dims[1],nv), dtype=np.float32) |
---|
816 | for j in range(1,latinf.dims[0]-1): |
---|
817 | for i in range(1,loninf.dims[1]-1): |
---|
818 | lon_bounds[j,i,0] = olon[j,i] - (olon[j+1,i-1] - olon[j,i])*0.5 |
---|
819 | lon_bounds[j,i,1] = olon[j,i] - (olon[j+1,i+1] - olon[j,i])*0.5 |
---|
820 | lon_bounds[j,i,2] = olon[j,i] - (olon[j-1,i+1] - olon[j,i])*0.5 |
---|
821 | lon_bounds[j,i,3] = olon[j,i] - (olon[j-1,i-1] - olon[j,i])*0.5 |
---|
822 | lat_bounds[j,i,0] = olat[j,i] - (olat[j+1,i-1] - olat[j,i])*0.5 |
---|
823 | lat_bounds[j,i,1] = olat[j,i] - (olat[j+1,i+1] - olat[j,i])*0.5 |
---|
824 | lat_bounds[j,i,2] = olat[j,i] - (olat[j-1,i+1] - olat[j,i])*0.5 |
---|
825 | lat_bounds[j,i,3] = olat[j,i] - (olat[j-1,i-1] - olat[j,i])*0.5 |
---|
826 | |
---|
827 | i = loninf.dims[1]-1 |
---|
828 | for j in range(1,latinf.dims[0]): |
---|
829 | lon_bounds[j,i,0] = olon[j,i] - (olon[j+1,i-1] - olon[j,i])*0.5 |
---|
830 | lon_bounds[j,i,1] = olon[j,i] + (olon[j+1,i-1] - olon[j,i])*0.5 |
---|
831 | lon_bounds[j,i,2] = olon[j,i] + (olon[j-1,i-1] - olon[j,i])*0.5 |
---|
832 | lon_bounds[j,i,3] = olon[j,i] - (olon[j-1,i-1] - olon[j,i])*0.5 |
---|
833 | j = 0 |
---|
834 | lon_bounds[j,i,0] = olon[j,i] - np.abs(olon[j+1,i-1] - olon[j,i])*0.5 |
---|
835 | lon_bounds[j,i,1] = olon[j,i] + np.abs(olon[j+1,i-1] - olon[j,i])*0.5 |
---|
836 | lon_bounds[j,i,2] = olon[j,i] + np.abs(olon[j+1,i-1] - olon[j,i])*0.5 |
---|
837 | lon_bounds[j,i,3] = olon[j,i] - np.abs(olon[j+1,i-1] - olon[j,i])*0.5 |
---|
838 | j = latinf.dims[0] |
---|
839 | lon_bounds[j,i,0] = olon[j,i] - np.abs(olon[j-1,i-1] - olon[j,i])*0.5 |
---|
840 | lon_bounds[j,i,1] = olon[j,i] + np.abs(olon[j-1,i-1] - olon[j,i])*0.5 |
---|
841 | lon_bounds[j,i,2] = olon[j,i] + np.abs(olon[j-1,i-1] - olon[j,i])*0.5 |
---|
842 | lon_bounds[j,i,3] = olon[j,i] - np.abs(olon[j-1,i-1] - olon[j,i])*0.5 |
---|
843 | |
---|
844 | j = latinf.dims[0]-1 |
---|
845 | for i in range(1,loninf.dims[1]): |
---|
846 | lat_bounds[j,i,0] = olat[j,i] + (olat[j-1,i-1] - olat[j,i])*0.5 |
---|
847 | lat_bounds[j,i,1] = olat[j,i] + (olat[j-1,i+1] - olat[j,i])*0.5 |
---|
848 | lat_bounds[j,i,2] = olat[j,i] - (olat[j-1,i+1] - olat[j,i])*0.5 |
---|
849 | lat_bounds[j,i,3] = olat[j,i] - (olat[j-1,i-1] - olat[j,i])*0.5 |
---|
850 | i = 0 |
---|
851 | lat_bounds[j,i,0] = olat[j,i] + np.abs(olat[j-1,i+1] - olat[j,i])*0.5 |
---|
852 | lat_bounds[j,i,1] = olat[j,i] + np.abs(olat[j-1,i+1] - olat[j,i])*0.5 |
---|
853 | lat_bounds[j,i,2] = olat[j,i] - np.abs(olat[j-1,i+1] - olat[j,i])*0.5 |
---|
854 | lat_bounds[j,i,3] = olat[j,i] - np.abs(olat[j-1,i+1] - olat[j,i])*0.5 |
---|
855 | i = loninf.dims[1] |
---|
856 | lat_bounds[j,i,0] = olat[j,i] + np.abs(olat[j-1,i-1] - olat[j,i])*0.5 |
---|
857 | lat_bounds[j,i,1] = olat[j,i] + np.abs(olat[j-1,i-1] - olat[j,i])*0.5 |
---|
858 | lat_bounds[j,i,2] = olat[j,i] - np.abs(olat[j-1,i-1] - olat[j,i])*0.5 |
---|
859 | lat_bounds[j,i,3] = olat[j,i] - np.abs(olat[j-1,i-1] - olat[j,i])*0.5 |
---|
860 | i = 0 |
---|
861 | j = 0 |
---|
862 | lon_bounds[j,i,0] = olon[j,i] - np.abs(olon[j+1,i+1] - olon[j,i])*0.5 |
---|
863 | lon_bounds[j,i,1] = olon[j,i] + np.abs(olon[j+1,i+1] - olon[j,i])*0.5 |
---|
864 | lon_bounds[j,i,2] = olon[j,i] + np.abs(olon[j+1,i+1] - olon[j,i])*0.5 |
---|
865 | lon_bounds[j,i,3] = olon[j,i] - np.abs(olon[j+1,i+1] - olon[j,i])*0.5 |
---|
866 | lat_bounds[j,i,0] = olat[j,i] + np.abs(olat[j+1,i+1] - olat[j,i])*0.5 |
---|
867 | lat_bounds[j,i,1] = olat[j,i] + np.abs(olat[j+1,i+1] - olat[j,i])*0.5 |
---|
868 | lat_bounds[j,i,2] = olat[j,i] - np.abs(olat[j+1,i+1] - olat[j,i])*0.5 |
---|
869 | lat_bounds[j,i,3] = olat[j,i] - np.abs(olat[j+1,i+1] - olat[j,i])*0.5 |
---|
870 | |
---|
871 | newvar = newncobj.createVariable('lon_bnds', 'f8', ('lat', 'lon','nv')) |
---|
872 | newvar[:] = lon_bounds[:] |
---|
873 | newvar = newncobj.createVariable('lat_bnds', 'f8', ('lat', 'lon','nv')) |
---|
874 | newvar[:] = lat_bounds[:] |
---|
875 | |
---|
876 | # Variable |
---|
877 | ovar = ncobj.variables[opts.varname] |
---|
878 | varinf = variable_inf(ovar) |
---|
879 | varvals = ovar[:] |
---|
880 | |
---|
881 | if varinf.Ndims == 3 and (opts.typen is None or opts.timen is None): |
---|
882 | print errormsg |
---|
883 | print ' ' + main + ': curent variable is 3D !!' |
---|
884 | print ' var dim names:',varinf.dimns |
---|
885 | print ' var dims lengths:',varinf.dims |
---|
886 | print ' current variables in file:', list(ncobj.variables) |
---|
887 | print " either provide the values for the type (with '-T [typevals]' flag)" |
---|
888 | print " either provide the values for the time (with '-I [timevals]' flag)" |
---|
889 | quit(-1) |
---|
890 | elif varinf.Ndims == 4 and (opts.typen is None or opts.timen is None): |
---|
891 | print errormsg |
---|
892 | print ' ' + main + ': curent variable is 4D !!' |
---|
893 | print ' var dim names:',varinf.dimns |
---|
894 | print ' var dims lengths:',varinf.dims |
---|
895 | print ' current variables in file:', list(ncobj.variables) |
---|
896 | print " provide the values for the type (with '-T [typevals]' flag)" |
---|
897 | print " and provide the values for the time (with '-I [timevals]' flag)" |
---|
898 | quit(-1) |
---|
899 | |
---|
900 | if opts.range is not None: |
---|
901 | varmin = retype(opts.range.split(',')[0], varinf.dtype) |
---|
902 | varmax = retype(opts.range.split(',')[1], varinf.dtype) |
---|
903 | range2 = retype(0.5*(varmax + varmin), varinf.dtype) |
---|
904 | else: |
---|
905 | print warnmsg |
---|
906 | print ' ' + main + ': no range of valid values is provided!!' |
---|
907 | print ' assuming validity of all values' |
---|
908 | range2 = retype(0.5*(np.min(varvals) + np.max(varvals)), varinf.dtype) |
---|
909 | |
---|
910 | dimx=varinf.dims[varinf.Ndims-1] |
---|
911 | dimy=varinf.dims[varinf.Ndims-2] |
---|
912 | |
---|
913 | missf = False |
---|
914 | if searchInlist(varinf.attributes,'missing_value'): |
---|
915 | missf = True |
---|
916 | missv = retype(ovar.getncattr('missing_value'), varinf.dtype) |
---|
917 | print warnmsg |
---|
918 | print ' ' + main + ": variable '" + opts.varname + "' has missing value!!" |
---|
919 | if type(varvals) != type(ma.array([1, 2, 3])): |
---|
920 | print warnmsg |
---|
921 | print ' ' + main + ": variable '" + opts.varname + "' is not masked!!" |
---|
922 | print " masking it using the 'missing_value':", missv |
---|
923 | varvals = ma.masked_equal(varvals, missv) |
---|
924 | |
---|
925 | # Types |
---|
926 | typevals = None |
---|
927 | if opts.typen is not None: |
---|
928 | if typen != 'None': |
---|
929 | if not searchInlist(ncobj.variables, typen): |
---|
930 | print errormsg |
---|
931 | print ' ' + main + ": file '" + opts.filen + "' does not have type " + \ |
---|
932 | "variable '" + typen + "' !!" |
---|
933 | quit(-1) |
---|
934 | otype = ncobj.variables[typen] |
---|
935 | if len(otype.shape) != 2: |
---|
936 | print warnmsg |
---|
937 | print ' ' + main + ": file provided types variable '" + typen + \ |
---|
938 | "' does not have the expected shape (Ntypes, LString)!!" |
---|
939 | print " please, remove the vriable (as None,[typedimn]) from the " + \ |
---|
940 | "'-T [typevals] flag'" |
---|
941 | quit(-1) |
---|
942 | else: |
---|
943 | Ntypes = otype.shape[0] |
---|
944 | Lstring = otype.shape[1] |
---|
945 | |
---|
946 | newdim = newncobj.createDimension('Ntypes',Ntypes) |
---|
947 | newdim = newncobj.createDimension('Lstring',Lstring) |
---|
948 | typevalues = otype[:] |
---|
949 | typevals = [] |
---|
950 | for ityp in range(typevalues.shape[0]): |
---|
951 | string = '' |
---|
952 | for il in ityp: |
---|
953 | string = string + ityp[il:il+1] |
---|
954 | typevals.append(string) |
---|
955 | else: |
---|
956 | Ntypes = varinf.dims[varinf.Ndims-3] |
---|
957 | Lstring = 250 |
---|
958 | newdim = newncobj.createDimension('Ntypes',Ntypes) |
---|
959 | newdim = newncobj.createDimension('Lstring',Lstring) |
---|
960 | if opts.typesfile is None: |
---|
961 | print errormsg |
---|
962 | print ' ' + main + ": variable has a types dimension '" + typedimn + \ |
---|
963 | "' but there is not a variable with values !!" |
---|
964 | print ' provide an ASCII file with the required types values with ' + \ |
---|
965 | "the flag '-t [ASCII_file]' with the python call" |
---|
966 | quit(-1) |
---|
967 | else: |
---|
968 | typevals = read_ASCIIlist(opts.typesfile) |
---|
969 | if len(typevals) != Ntypes: |
---|
970 | print errormsg |
---|
971 | print ' ' + main + ': number of values:', len(typevals), \ |
---|
972 | "from the ASCII file '" + opts.typesfile + "' does not work " + \ |
---|
973 | 'for the size:', varinf.dims, 'of the current variable!!' |
---|
974 | quit(-1) |
---|
975 | else: |
---|
976 | if opts.typesfile is not None: |
---|
977 | if not os.path.isfile(opts.typesfile): |
---|
978 | print errormsg |
---|
979 | print ' ' + main + ": ASCII types file '" + opts.typesfile + \ |
---|
980 | "' does not exist !!" |
---|
981 | quit(-1) |
---|
982 | typevals = read_ASCIIlist(opts.typesfile) |
---|
983 | Ntypes = len(typevals) |
---|
984 | Lstring = 250 |
---|
985 | newdim = newncobj.createDimension('Ntypes',Ntypes) |
---|
986 | newdim = newncobj.createDimension('Lstring',Lstring) |
---|
987 | |
---|
988 | # Time |
---|
989 | timevals = None |
---|
990 | if opts.timen is not None: |
---|
991 | if timen != 'None': |
---|
992 | if not searchInlist(ncobj.variables, timen): |
---|
993 | print errormsg |
---|
994 | print ' ' + main + ": file '" + opts.ncfile + "' does not have time " + \ |
---|
995 | " variable '" + timen + "' !!" |
---|
996 | quit(-1) |
---|
997 | otime = ncobj.variables[timen] |
---|
998 | timeinf = variable_inf(otime) |
---|
999 | |
---|
1000 | dimt = otime.shape[0] |
---|
1001 | newdim = newncobj.createDimension('time', dimt) |
---|
1002 | timevals = otime[:] |
---|
1003 | if timeinf.lname is None or timeinf.units is None or \ |
---|
1004 | not searchInlist(timeinf.attribtues, 'calendar'): |
---|
1005 | print warnmsg |
---|
1006 | print ' ' + main + ": current time varibale in file '" + timen + \ |
---|
1007 | "' does not have units!!" |
---|
1008 | if opts.timeattrs is None: |
---|
1009 | print " add flag '-Y [long_name],[timeunits],[calendar]' (@ for " + \ |
---|
1010 | 'spaces) to the python call to use it as the long_name and the ' + \ |
---|
1011 | 'units for the time variable' |
---|
1012 | quit(-1) |
---|
1013 | else: |
---|
1014 | tlname = opts.timeattrs.split(',')[0].replace('@',' ') |
---|
1015 | tunits = opts.timeattrs.split(',')[1].replace('@',' ') |
---|
1016 | tcalendar = opts.timeattrs.split(',')[2].replace('@',' ') |
---|
1017 | else: |
---|
1018 | tlname = timeinf.lname |
---|
1019 | tunits = timeinf.units |
---|
1020 | tcalendar = otime.getncattr('calendar') |
---|
1021 | else: |
---|
1022 | dimt = varinf.dims[0] |
---|
1023 | newdim = newncobj.createDimension('time', dimt) |
---|
1024 | timevals = np.arange((dimt), dtype=np.float64) |
---|
1025 | if opts.timeattrs is None: |
---|
1026 | print " add flag '-Y [long_name],[timeunits]' (@ for spaces) to " + \ |
---|
1027 | 'the python call to use it as the long_name and the units for '+ \ |
---|
1028 | ' the time variable' |
---|
1029 | quit(-1) |
---|
1030 | else: |
---|
1031 | tlname = opts.timeattrs.split(',')[0].replace('@',' ') |
---|
1032 | tunits = opts.timeattrs.split(',')[1].replace('@',' ') |
---|
1033 | tcalendar = opts.timeattrs.split(',')[2].replace('@',' ') |
---|
1034 | tu = tunits.split(' ')[0] |
---|
1035 | if tu[0:5] == 'month': |
---|
1036 | print warnmsg |
---|
1037 | print ' ' + main + ": changing time units from '" + tu + \ |
---|
1038 | "' to 'days' !!" |
---|
1039 | tunits = tunits.replace('months', 'days').replace('month','days') |
---|
1040 | daysmonth = np.array([0., 31., 28., 31., 30., 31., 30., 31., 31., \ |
---|
1041 | 30., 31., 30.], dtype=np.float32) |
---|
1042 | timevals = np.zeros((dimt), dtype=np.float32) |
---|
1043 | for it in range(1,dimt): |
---|
1044 | timevals[it] = np.sum(daysmonth[0:it+1]) |
---|
1045 | |
---|
1046 | lons, lats = lonlat2D(ncobj, lonn, latn) |
---|
1047 | |
---|
1048 | # Variable |
---|
1049 | print " Creation of variable '" + opts.varname + "' ..." |
---|
1050 | |
---|
1051 | fillvalt = fillvalue_kind(varinf.dtype, 'std') |
---|
1052 | |
---|
1053 | if varinf.FillValue is None: |
---|
1054 | print warnmsg |
---|
1055 | print " variable does not have 'fill_value' attribute!!" |
---|
1056 | if opts.oceanval is None: |
---|
1057 | print warnmsg |
---|
1058 | print " you can add flag '-z [oceanval]' to the python call to use it " + \ |
---|
1059 | ' as the current value to mask the ocean values' |
---|
1060 | oceanval = None |
---|
1061 | else: |
---|
1062 | oceanval = retype(opts.oceanval, varinf.dtype) |
---|
1063 | print ' using:',oceanval,'as a complementary way to generate the ocean mask' |
---|
1064 | else: |
---|
1065 | oceanval = retype(opts.oceanval, varinf.dtype) |
---|
1066 | if oceanval != varinf.FillValue: |
---|
1067 | print warnmsg |
---|
1068 | print ' current fill_value:', varinf.FillValue,' and ocean values:', \ |
---|
1069 | oceanval,' do not coincide!!' |
---|
1070 | print ' using ocean value to create the mask' |
---|
1071 | |
---|
1072 | slicemask = [] |
---|
1073 | for idim in range(varinf.Ndims-2): |
---|
1074 | slicemask.append(0) |
---|
1075 | slicemask.append(slice(0,varinf.dims[varinf.Ndims-2])) |
---|
1076 | slicemask.append(slice(0,varinf.dims[varinf.Ndims-1])) |
---|
1077 | |
---|
1078 | # Mask |
---|
1079 | print " Creation of 'mask' ..." |
---|
1080 | newvar = newncobj.createVariable('mask', 'i4', ('lat','lon')) |
---|
1081 | maskn = opts.maskn.split('@')[0] |
---|
1082 | if maskn == 'maskfile': |
---|
1083 | maskf = opts.maskn.split('@')[1] |
---|
1084 | maskv = opts.maskn.split('@')[2] |
---|
1085 | maskl = opts.maskn.split('@')[3] |
---|
1086 | maskL = opts.maskn.split('@')[4] |
---|
1087 | |
---|
1088 | print warnmsg |
---|
1089 | print ' ' + main + ": taking mask values from file '" + maskf + "' variable '" +\ |
---|
1090 | maskv + "' !!" |
---|
1091 | print " make sure that it is on the same projection and resolution as '" + \ |
---|
1092 | opts.ncfile + "' !!" |
---|
1093 | |
---|
1094 | compncobj = NetCDFFile(maskf,'r') |
---|
1095 | if not searchInlist(compncobj.variables,maskv): |
---|
1096 | print errormsg |
---|
1097 | print ' ' + fname + ": file '" + maskf + '" does not have mask variable "' +\ |
---|
1098 | maskn + "' !!" |
---|
1099 | quit(-1) |
---|
1100 | if not searchInlist(compncobj.variables,maskl): |
---|
1101 | print errormsg |
---|
1102 | print ' ' + fname + ": file '" + maskf + '" does not have lon variable "' + \ |
---|
1103 | maskl + "' !!" |
---|
1104 | quit(-1) |
---|
1105 | if not searchInlist(compncobj.variables,maskL): |
---|
1106 | print errormsg |
---|
1107 | print ' ' + fname + ": file '" + maskf + '" does not have lat variable "' + \ |
---|
1108 | maskL + "' !!" |
---|
1109 | quit(-1) |
---|
1110 | |
---|
1111 | omask = compncobj.variables[maskv] |
---|
1112 | maskvals = omask[:] |
---|
1113 | complons, complats = lonlat2D(compncobj, maskl, maskL) |
---|
1114 | |
---|
1115 | # Checking a bit |
---|
1116 | if maskvals.shape[1] != dimx: |
---|
1117 | print errormsg |
---|
1118 | print ' ' + main + ": complementary mask file '" + maskf + "' x-shape ", \ |
---|
1119 | maskvals.shape[1], "is different than '" + opts.ncfile + "': ", dimx,'!!' |
---|
1120 | quit(-1) |
---|
1121 | if maskvals.shape[0] != dimy: |
---|
1122 | print errormsg |
---|
1123 | print ' ' + main + ": complementary mask file '" + maskf + "' y-shape ", \ |
---|
1124 | maskvals.shape[0], "is different than '" + opts.ncfile + "': ", dimy,'!!' |
---|
1125 | quit(-1) |
---|
1126 | if not np.all(complons - lons == 0.): |
---|
1127 | print errormsg |
---|
1128 | print ' ' + main + ": complementary mask file '" + maskf + "' and '" + \ |
---|
1129 | opts.ncfile + "' have different longitude projection !!" |
---|
1130 | quit(-1) |
---|
1131 | if not np.all(complats - lats == 0.): |
---|
1132 | print errormsg |
---|
1133 | print ' ' + main + ": complementary mask file '" + maskf + "' and '" + \ |
---|
1134 | opts.ncfile + "' have different latitude projection !!" |
---|
1135 | quit(-1) |
---|
1136 | compncobj.close() |
---|
1137 | |
---|
1138 | elif maskn == 'maskvarname': |
---|
1139 | maskv = opts.maskn.split('@')[1] |
---|
1140 | if not searchInlist(ncobj.variables,maskv): |
---|
1141 | print errormsg |
---|
1142 | print ' ' + fname + ": file '" + opts.ncfile + '" does not have variable "'+\ |
---|
1143 | maskn + "' !!" |
---|
1144 | quit(-1) |
---|
1145 | else: |
---|
1146 | omask = ncobj.variables[maskv] |
---|
1147 | maskvals = omask[:] |
---|
1148 | elif maskn == 'maskbelow': |
---|
1149 | maskv = retype(opts.maskn.split('@')[1], varinf.dtype) |
---|
1150 | if missf: |
---|
1151 | varvalsmask = varvals.filled(range2).reshape(np.prod(varvals.shape[0: \ |
---|
1152 | varinf.Ndims-2]),dimy,dimx) |
---|
1153 | else: |
---|
1154 | varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \ |
---|
1155 | dimx) |
---|
1156 | maskvals0 = np.where(np.all(varvalsmask <= maskv, axis=0), 1, 0) |
---|
1157 | if oceanval is not None: |
---|
1158 | maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1) |
---|
1159 | maskvals0 = maskvals0 * maskvals00 |
---|
1160 | # maskvals = maskvals0[tuple(slicemask)] |
---|
1161 | maskvals = maskvals0 |
---|
1162 | elif maskn == 'maskabove': |
---|
1163 | maskv = retype(opts.maskn.split('@')[1], varinf.dtype) |
---|
1164 | if missf: |
---|
1165 | varvalsmask = varvals.filled(range2).reshape(np.prod(varvals.shape[0: \ |
---|
1166 | varinf.Ndims-2]),dimy,dimx) |
---|
1167 | else: |
---|
1168 | varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \ |
---|
1169 | dimx) |
---|
1170 | maskvals0 = np.where(np.all(varvalsmask >= maskv, axis=0), 1, 0) |
---|
1171 | if oceanval is not None: |
---|
1172 | maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1) |
---|
1173 | maskvals0 = maskvals0 * maskvals00 |
---|
1174 | # maskvals = maskvals0[tuple(slicemask)] |
---|
1175 | maskvals = maskvals0 |
---|
1176 | elif maskn == 'maskrange': |
---|
1177 | maskv1 = retype(opts.maskn.split('@')[1].split(',')[0], varinf.dtype) |
---|
1178 | maskv2 = retype(opts.maskn.split('@')[1].split(',')[1], varinf.dtype) |
---|
1179 | if missf: |
---|
1180 | varvalsmask = varvals.filled(range2).reshape(np.prod( \ |
---|
1181 | varvals.shape[0:varinf.Ndims-2]),dimy,dimx) |
---|
1182 | else: |
---|
1183 | varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \ |
---|
1184 | dimx) |
---|
1185 | maskvals01 = np.where(np.all(varvalsmask >= maskv1, axis=0), 1, 0) |
---|
1186 | maskvals02 = np.where(np.all(varvalsmask <= maskv2, axis=0), 1, 0) |
---|
1187 | maskvals0 = maskvals01 * maskvals02 |
---|
1188 | if oceanval is not None: |
---|
1189 | maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1) |
---|
1190 | maskvals0 = maskvals0 * maskvals00 |
---|
1191 | # maskvals = maskvals0[tuple(slicemask)] |
---|
1192 | maskvals = maskvals0 |
---|
1193 | elif maskn == 'masksumrange': |
---|
1194 | maskv1 = retype(opts.maskn.split('@')[1].split(',')[0], varinf.dtype) |
---|
1195 | maskv2 = retype(opts.maskn.split('@')[1].split(',')[1], varinf.dtype) |
---|
1196 | maskv3 = retype(opts.maskn.split('@')[1].split(',')[2], varinf.dtype) |
---|
1197 | if missf: |
---|
1198 | varvalsmask = varvals.filled(range2).reshape(np.prod(varvals.shape[0: \ |
---|
1199 | varinf.Ndims-2]),dimy,dimx) |
---|
1200 | else: |
---|
1201 | varvalsmask = varvals.reshape(np.prod(varvals.shape[0:varinf.Ndims-2]),dimy, \ |
---|
1202 | dimx) |
---|
1203 | print ' masksumrange: _______' |
---|
1204 | print ' ', maskv2,'<= varsum <=',maskv1 |
---|
1205 | print ' ', maskv1,'< varsum <=',maskv3 |
---|
1206 | varsum = np.sum(varvals, axis=varinf.Ndims-3) |
---|
1207 | maskvals1a = np.where(varsum >= maskv2, 1, 0) |
---|
1208 | maskvals1b = np.where(varsum <= maskv1, 1, 0) |
---|
1209 | maskvals2a = np.where(varsum > maskv1, 1, 0) |
---|
1210 | maskvals2b = np.where(varsum <= maskv3, 1, 0) |
---|
1211 | if varinf.Ndims == 4: |
---|
1212 | maskvals0 = maskvals1a[0,:,:]*maskvals1b[0,:,:]+maskvals2a[0,:,:]* \ |
---|
1213 | maskvals2b[0,:,:] |
---|
1214 | elif varinf.Ndims == 3: |
---|
1215 | maskvals0 = maskvals1a*maskvals1b+maskvals2a*maskvals2b |
---|
1216 | else: |
---|
1217 | print errormsg |
---|
1218 | print ' ' + main + ": masking by '" + maskn + "' with shape:", varinf.dims, \ |
---|
1219 | 'not ready!!' |
---|
1220 | quit(-1) |
---|
1221 | if oceanval is not None: |
---|
1222 | maskvals00 = np.where(np.all(varvalsmask == oceanval, axis=0), 0, 1) |
---|
1223 | maskvals0 = maskvals0 * maskvals00 |
---|
1224 | if type(maskvals0) == type(ma.masked_array(np.arange(1))): |
---|
1225 | msk = maskvals0.mask |
---|
1226 | maskvals = np.where(msk, 0, 1) |
---|
1227 | else: |
---|
1228 | maskvals = maskvals0 |
---|
1229 | elif maskn == 'maskmissing': |
---|
1230 | if missf: |
---|
1231 | maskvals = np.where(varvals.mask[tuple(slicemask)], 0, 1) |
---|
1232 | else: |
---|
1233 | print errormsg |
---|
1234 | print ' ' + main + ": it should be masked following the missing value '" + \ |
---|
1235 | mask + "' but variable does not have attribute 'missing_value'!!" |
---|
1236 | quit(-1) |
---|
1237 | else: |
---|
1238 | print errormsg |
---|
1239 | print ' ' + main + ": mask type '" + maskn + "' not ready!!" |
---|
1240 | quit(-1) |
---|
1241 | |
---|
1242 | # This should not be needed, isn't? |
---|
1243 | #if varinf.FillValue is not None: |
---|
1244 | # msk = msk * varvals.msk |
---|
1245 | # maskvals = np.where(msk, 0, 1) |
---|
1246 | newvar[:] = maskvals[:] |
---|
1247 | |
---|
1248 | # Dimensions new varibale |
---|
1249 | vardims0 = [] |
---|
1250 | for dimn in varinf.dimns: |
---|
1251 | if dimn == londimn: vardims0.append('lon') |
---|
1252 | if dimn == latdimn: vardims0.append('lat') |
---|
1253 | if opts.typen is not None and dimn == typedimn: vardims0.append('Ntypes') |
---|
1254 | if opts.timen is not None and dimn == timedimn: vardims0.append('time') |
---|
1255 | if not searchInlist(vardims0, 'lon') or not searchInlist(vardims0, 'lat'): |
---|
1256 | print errormsg |
---|
1257 | print ' ' + main + ": new variable will not have mandatory 'lon', 'lat' " + \ |
---|
1258 | "dimensions!!" |
---|
1259 | print ' check given names of the dimensions:', londimn, ',', latdimn |
---|
1260 | quit(-1) |
---|
1261 | if len(vardims0) != varinf.Ndims: |
---|
1262 | print errormsg |
---|
1263 | print ' ' + main + ': new variable will not have same number of dimensions:', \ |
---|
1264 | len(vardims0), 'as original one', varinf.Ndims,'!!' |
---|
1265 | print ' check given names of dimensions lon, lat:', londimn, ',', latdimn |
---|
1266 | if opts.typen is not None: print ' types:', typedimn |
---|
1267 | if opts.timen is not None: print ' time:', timedimn |
---|
1268 | quit(-1) |
---|
1269 | |
---|
1270 | vardims = tuple(vardims0) |
---|
1271 | |
---|
1272 | # Masking variable values according to the mask and the range |
---|
1273 | varmask0 = np.where(maskvals == 0, True, False) |
---|
1274 | varmask = np.zeros(varinf.dims, dtype=bool) |
---|
1275 | varmask[...,:,:] = varmask0 |
---|
1276 | if varmin is not None: |
---|
1277 | if np.min(ovar[:]) < varmin: |
---|
1278 | Nvals = np.sum(varvals < varmin) |
---|
1279 | print warnmsg |
---|
1280 | print ' ' + main + ': masking ', Nvals,'values of variable lower than ', \ |
---|
1281 | 'varmin:', varmin |
---|
1282 | if np.max(ovar[:]) > varmax: |
---|
1283 | Nvals = np.sum(varvals > varmax) |
---|
1284 | print warnmsg |
---|
1285 | print ' ' + main + ': masking ', Nvals,'values of variable greater than ', \ |
---|
1286 | 'varmax:', varmax |
---|
1287 | |
---|
1288 | if missf: |
---|
1289 | varvalsfill = varvals.filled(varmax + range2) |
---|
1290 | minrangevarval = ma.masked_less(varvalsfill, varmin) |
---|
1291 | minrangevarval = np.where(minrangevarval == varmax+range2, missv, \ |
---|
1292 | minrangevarval) |
---|
1293 | else: |
---|
1294 | minrangevarval = ma.masked_less(varvals, varmin) |
---|
1295 | |
---|
1296 | if missf: |
---|
1297 | varvalsfill = varvals.filled(varmin - range2) |
---|
1298 | maxrangevarval = ma.masked_greater(varvalsfill, varmax) |
---|
1299 | maxrangevarval = np.where(maxrangevarval == varmax-range2, missv, \ |
---|
1300 | maxrangevarval) |
---|
1301 | else: |
---|
1302 | maxrangevarval = ma.masked_greater(varvals, varmax) |
---|
1303 | |
---|
1304 | rangevarmask = minrangevarval.mask*maxrangevarval.mask |
---|
1305 | |
---|
1306 | if missf: |
---|
1307 | fillvarval = ma.array(varvals.filled(missv), mask=rangevarmask) |
---|
1308 | else: |
---|
1309 | fillvarval = ma.array(varvals, mask=rangevarmask) |
---|
1310 | else: |
---|
1311 | fillvarval = varvals |
---|
1312 | |
---|
1313 | fillvarval = ma.array(fillvarval, mask=varmask) |
---|
1314 | |
---|
1315 | newvar = newncobj.createVariable(opts.varname, nctype(varinf.dtype), vardims, \ |
---|
1316 | fill_value=fillvalt) |
---|
1317 | newvar[:] = fillvarval[:] |
---|
1318 | |
---|
1319 | if varinf.sname is None: |
---|
1320 | print warnmsg |
---|
1321 | print " variable does not have 'standard_name' attribute!!" |
---|
1322 | if opts.svname is None: |
---|
1323 | print errormsg |
---|
1324 | print " add flag '-S [standard_name]' to the python call to assign the" + \ |
---|
1325 | ' attribute' |
---|
1326 | quit(-1) |
---|
1327 | else: |
---|
1328 | newattr = set_attributek(newvar, 'standard_name', opts.svname, 'S') |
---|
1329 | else: |
---|
1330 | newattr = set_attributek(newvar, 'standard_name', varinf.sname, 'S') |
---|
1331 | |
---|
1332 | if varinf.lname is None: |
---|
1333 | print warnmsg |
---|
1334 | print " variable does not have 'long_name' attribute!!" |
---|
1335 | if opts.lvname is None: |
---|
1336 | print errormsg |
---|
1337 | print " add flag '-L [long_name]' (@ for spaces) to the python call to" + \ |
---|
1338 | 'assign the attribute' |
---|
1339 | quit(-1) |
---|
1340 | else: |
---|
1341 | newattr = set_attributek(newvar, 'long_name', opts.lvname.replace('@',' '), \ |
---|
1342 | 'S') |
---|
1343 | else: |
---|
1344 | newattr = set_attributek(newvar, 'long_name', varinf.lname, 'S') |
---|
1345 | |
---|
1346 | if varinf.units is None: |
---|
1347 | print warnmsg |
---|
1348 | print " variable does not have 'units' attribute!!" |
---|
1349 | if opts.vunits is None: |
---|
1350 | print errormsg |
---|
1351 | print " add flag '-U [units]' to the python call to assign the attribute" |
---|
1352 | quit(-1) |
---|
1353 | else: |
---|
1354 | newattr = set_attributek(newvar, 'units', opts.vunits, 'S') |
---|
1355 | else: |
---|
1356 | newattr = set_attributek(newvar, 'units', varinf.units, 'S') |
---|
1357 | |
---|
1358 | newattr = set_attributek(newvar, 'coordinates', 'lon lat', 'S') |
---|
1359 | if opts.range is not None: |
---|
1360 | newattr = set_attributek(newvar, 'valid_min', varmin, varinf.dtype) |
---|
1361 | newattr = set_attributek(newvar, 'valid_max', varmax, varinf.dtype) |
---|
1362 | |
---|
1363 | if opts.extrattrs is not None: |
---|
1364 | if opts.extrattrs.find(',') != -1: |
---|
1365 | attrs = opts.extrattrs.split(',') |
---|
1366 | else: |
---|
1367 | attrs = [opts.extrattrs] |
---|
1368 | for attr in attrs: |
---|
1369 | if len(attr.split(':')) != 3: |
---|
1370 | print errormsg |
---|
1371 | print ' ' + main + ": wrong number of values:", len(attr.split(':')), \ |
---|
1372 | 'for extra attribute !!' |
---|
1373 | print " given values:'" + attr + "' 3 ':' separated are expected" |
---|
1374 | print ' which gives:', attr.split(':') |
---|
1375 | quit(-1) |
---|
1376 | attrn = attr.split(':')[0] |
---|
1377 | attrv = attr.split(':')[1] |
---|
1378 | attrk = attr.split(':')[2] |
---|
1379 | newattr = set_attributek(newvar, attrn, attrv, attrk) |
---|
1380 | |
---|
1381 | newncobj.sync() |
---|
1382 | # Types variable |
---|
1383 | if typevals is not None: |
---|
1384 | print " Creation of 'typenames' ..." |
---|
1385 | newvar = newncobj.createVariable('typenames', 'c', ('Ntypes', 'Lstring')) |
---|
1386 | writing_str_nc(newvar, typevals, Lstring) |
---|
1387 | |
---|
1388 | # Time variable |
---|
1389 | if timevals is not None: |
---|
1390 | print " Creation of 'time' ..." |
---|
1391 | newvar = newncobj.createVariable('time', 'f8', ('time')) |
---|
1392 | newattr = set_attributek(newvar, 'standard_name', 'time', 'S') |
---|
1393 | newattr = set_attributek(newvar, 'long_name', tlname, 'S') |
---|
1394 | newattr = set_attributek(newvar, 'units', tunits, 'S') |
---|
1395 | newattr = set_attributek(newvar, 'calendar', tcalendar, 'S') |
---|
1396 | newvar[:] = timevals[:] |
---|
1397 | newattr = set_attributek(newvar, 'axis', 'T', 'S') |
---|
1398 | newattr = set_attributek(newvar, '_CoordinateAxisType', 'time', 'S') |
---|
1399 | |
---|
1400 | |
---|
1401 | newncobj.sync() |
---|
1402 | |
---|
1403 | newattr = set_attributek(newncobj, 'creation',"file created with " + main, 'S') |
---|
1404 | newattr = set_attributek(newncobj, 'py_script_version',version, 'S') |
---|
1405 | newattr = set_attributek(newncobj, 'py_script_date',date, 'S') |
---|
1406 | newattr = set_attributek(newncobj, 'py_script_author',"L. Fita", 'S') |
---|
1407 | newattr = set_attributek(newncobj, 'py_script_author_institution', \ |
---|
1408 | unicode('Laboratoire de M' + unichr(233)+'t'+unichr(233)+'orologie Dynamique'), 'U') |
---|
1409 | |
---|
1410 | # Getting global attributes |
---|
1411 | ncgattrs = ncobj.ncattrs() |
---|
1412 | for attr in ncgattrs: |
---|
1413 | attrv = ncobj.getncattr(attr) |
---|
1414 | attrk = type(attrv) |
---|
1415 | newattr = set_attributek(newncobj, attr, attrv, attrk) |
---|
1416 | |
---|
1417 | ncobj.close() |
---|
1418 | newncobj.close() |
---|
1419 | print main + ": successfull written of file '" + ofile + "' !!" |
---|