1 | # -*- coding: iso-8859-15 -*- |
---|
2 | # Generic program to transfrom ASCII observational data in columns to a netcdf |
---|
3 | # L. Fita, LMD February 2015 |
---|
4 | ## e.g. # create_OBSnetcdf.py -c '#' -d ACAR/description.dat -e space -f ACAR/2012/10/ACAR_121018.asc -g true -t 19491201000000,seconds -k trajectory |
---|
5 | |
---|
6 | import numpy as np |
---|
7 | from netCDF4 import Dataset as NetCDFFile |
---|
8 | import os |
---|
9 | import re |
---|
10 | from optparse import OptionParser |
---|
11 | |
---|
12 | # version |
---|
13 | version=1.2 |
---|
14 | |
---|
15 | # Filling values for floats, integer and string |
---|
16 | fillValueF = 1.e20 |
---|
17 | fillValueI = -99999 |
---|
18 | fillValueS = '---' |
---|
19 | |
---|
20 | # Length of the string variables |
---|
21 | StringLength = 200 |
---|
22 | |
---|
23 | # Size of the map for the complementary variables/maps |
---|
24 | Ndim2D = 100 |
---|
25 | |
---|
26 | main = 'create_OBSnetcdf.py' |
---|
27 | errormsg = 'ERROR -- error -- ERROR -- error' |
---|
28 | warnmsg = 'WARNING -- warning -- WARNING -- warning' |
---|
29 | |
---|
30 | fillValue = 1.e20 |
---|
31 | |
---|
32 | def searchInlist(listname, nameFind): |
---|
33 | """ Function to search a value within a list |
---|
34 | listname = list |
---|
35 | nameFind = value to find |
---|
36 | >>> searInlist(['1', '2', '3', '5'], '5') |
---|
37 | True |
---|
38 | """ |
---|
39 | for x in listname: |
---|
40 | if x == nameFind: |
---|
41 | return True |
---|
42 | return False |
---|
43 | |
---|
44 | |
---|
45 | def typemod(value, typeval): |
---|
46 | """ Function to give back a value in a given dtype |
---|
47 | >>> print(typemod(8.2223, 'np.float64')) |
---|
48 | <type 'numpy.float64'> |
---|
49 | >>> print(typemod(8.2223, 'tuple')) |
---|
50 | <type 'tuple'> |
---|
51 | """ |
---|
52 | fname='typemod' |
---|
53 | |
---|
54 | if typeval == 'int' or typeval == 'I': |
---|
55 | return int(value) |
---|
56 | elif typeval == 'long': |
---|
57 | return long(value) |
---|
58 | elif typeval == 'float' or typeval == 'F' or typeval == 'R': |
---|
59 | return float(value) |
---|
60 | elif typeval == 'complex': |
---|
61 | return complex(value) |
---|
62 | elif typeval == 'str' or typeval == 'S': |
---|
63 | return str(value) |
---|
64 | elif typeval == 'bool': |
---|
65 | return bool(value) |
---|
66 | elif typeval == 'B': |
---|
67 | return Str_Bool(value) |
---|
68 | elif typeval == 'list': |
---|
69 | newval = [] |
---|
70 | newval.append(value) |
---|
71 | return newval |
---|
72 | elif typeval == 'dic': |
---|
73 | newval = {} |
---|
74 | newval[value] = value |
---|
75 | return newval |
---|
76 | elif typeval == 'tuple': |
---|
77 | newv = [] |
---|
78 | newv.append(value) |
---|
79 | newval = tuple(newv) |
---|
80 | return newval |
---|
81 | elif typeval == 'np.int8': |
---|
82 | return np.int8(value) |
---|
83 | elif typeval == 'np.int16': |
---|
84 | return np.int16(value) |
---|
85 | elif typeval == 'np.int32': |
---|
86 | return np.int32(value) |
---|
87 | elif typeval == 'np.int64': |
---|
88 | return np.int64(value) |
---|
89 | elif typeval == 'np.uint8': |
---|
90 | return np.uint8(value) |
---|
91 | elif typeval == 'np.uint16': |
---|
92 | return np.uint16(value) |
---|
93 | elif typeval == 'np.np.uint32': |
---|
94 | return np.uint32(value) |
---|
95 | elif typeval == 'np.uint64': |
---|
96 | return np.uint64(value) |
---|
97 | elif typeval == 'np.float' or typeval == 'R': |
---|
98 | return np.float(value) |
---|
99 | elif typeval == 'np.float16': |
---|
100 | return np.float16(value) |
---|
101 | elif typeval == 'np.float32': |
---|
102 | return np.float32(value) |
---|
103 | elif typeval == 'float32': |
---|
104 | return np.float32(value) |
---|
105 | elif typeval == 'np.float64' or typeval == 'D': |
---|
106 | return np.float64(value) |
---|
107 | elif typeval == 'float64': |
---|
108 | return np.float64(value) |
---|
109 | elif typeval == 'np.complex': |
---|
110 | return np.complex(value) |
---|
111 | elif typeval == 'np.complex64': |
---|
112 | return np.complex64(value) |
---|
113 | elif typeval == 'np.complex128': |
---|
114 | return np.complex128(value) |
---|
115 | else: |
---|
116 | print errormsg |
---|
117 | print fname + ': data type "' + typeval + '" is not ready !' |
---|
118 | print errormsg |
---|
119 | quit(-1) |
---|
120 | |
---|
121 | return |
---|
122 | |
---|
123 | def str_list_k(string, cdiv, kind): |
---|
124 | """ Function to obtain a list of types of values from a string giving a split character |
---|
125 | string= String from which to obtain a list ('None' for None) |
---|
126 | cdiv= character to use to split the string |
---|
127 | kind= kind of desired value (as string like: 'np.float', 'int', 'np.float64', ....) |
---|
128 | >>> str_list_k('1:@#:$:56', ':', 'S') |
---|
129 | ['1', '@#', '$', '56'] |
---|
130 | >>> str_list_k('1:3.4:12.3', ':', 'np.float64') |
---|
131 | [1.0, 3.3999999999999999, 12.300000000000001] |
---|
132 | """ |
---|
133 | fname = 'str_list' |
---|
134 | |
---|
135 | if string.find(cdiv) != -1: |
---|
136 | listv = string.split(cdiv) |
---|
137 | else: |
---|
138 | if string == 'None': |
---|
139 | listv = None |
---|
140 | else: |
---|
141 | listv = [string] |
---|
142 | |
---|
143 | if listv is not None: |
---|
144 | finalist = [] |
---|
145 | for lv in listv: |
---|
146 | finalist.append(typemod(lv, kind)) |
---|
147 | else: |
---|
148 | finalist = None |
---|
149 | |
---|
150 | return finalist |
---|
151 | |
---|
152 | |
---|
153 | def set_attribute(ncvar, attrname, attrvalue): |
---|
154 | """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists |
---|
155 | ncvar = object netcdf variable |
---|
156 | attrname = name of the attribute |
---|
157 | attrvalue = value of the attribute |
---|
158 | """ |
---|
159 | import numpy as np |
---|
160 | from netCDF4 import Dataset as NetCDFFile |
---|
161 | |
---|
162 | attvar = ncvar.ncattrs() |
---|
163 | if searchInlist(attvar, attrname): |
---|
164 | attr = ncvar.delncattr(attrname) |
---|
165 | |
---|
166 | attr = ncvar.setncattr(attrname, attrvalue) |
---|
167 | |
---|
168 | return ncvar |
---|
169 | |
---|
170 | def set_attributek(ncv, attrname, attrval, attrkind): |
---|
171 | """ Sets a value of an attribute of a netCDF variable with a kind. Removes previous attribute value if exists |
---|
172 | ncvar = object netcdf variable |
---|
173 | attrname = name of the attribute |
---|
174 | attrvalue = value of the attribute |
---|
175 | atrtrkind = kind of attribute: 'S', string ('!' as spaces); 'U', unicode ('!' as spaces); 'I', integer; |
---|
176 | 'Inp32', numpy integer 32; 'R', ot 'F' real; ' npfloat', np.float; 'D', np.float64 |
---|
177 | """ |
---|
178 | fname = 'set_attributek' |
---|
179 | validk = {'S': 'string', 'U': 'unicode', 'I': 'integer', \ |
---|
180 | 'Inp32': 'integer long (np.int32)', 'F': 'float', 'R': 'float', \ |
---|
181 | 'npfloat': 'np.float', 'D': 'float long (np.float64)'} |
---|
182 | |
---|
183 | if type(attrkind) == type('s'): |
---|
184 | if attrkind == 'S': |
---|
185 | attrvalue = str(attrval.replace('!', ' ')) |
---|
186 | elif attrkind == 'U': |
---|
187 | attrvalue = unicode(attrval.replace('!',' ')) |
---|
188 | elif attrkind == 'I': |
---|
189 | attrvalue = np.int(attrval) |
---|
190 | elif attrkind == 'R' or attrkind == 'F' : |
---|
191 | attrvalue = float(attrval) |
---|
192 | elif attrkind == 'npfloat': |
---|
193 | attrvalue = np.float(attrval) |
---|
194 | elif attrkind == 'D': |
---|
195 | attrvalue = np.float64(attrval) |
---|
196 | else: |
---|
197 | print errormsg |
---|
198 | print ' ' + fname + ": '" + attrkind + "' kind of attribute is not ready!" |
---|
199 | print ' valid values: _______' |
---|
200 | for key in validk.keys(): |
---|
201 | print ' ', key,':', validk[key] |
---|
202 | quit(-1) |
---|
203 | else: |
---|
204 | if attrkind == type(str('a')): |
---|
205 | attrvalue = str(attrval.replace('!', ' ')) |
---|
206 | elif attrkind == type(unicode('a')): |
---|
207 | attrvalue = unicode(attrval.replace('!',' ')) |
---|
208 | elif attrkind == type(np.int(1)): |
---|
209 | attrvalue = np.int(attrval) |
---|
210 | elif attrkind == np.dtype('i'): |
---|
211 | attrvalue = np.int32(attrval) |
---|
212 | elif attrkind == type(float(1.)): |
---|
213 | attrvalue = float(attrval) |
---|
214 | elif attrkind == type(np.float(1.)): |
---|
215 | attrvalue = np.float(attrval) |
---|
216 | elif attrkind == np.dtype('float32'): |
---|
217 | attrvalue = np.array(attrval, dtype='f') |
---|
218 | elif attrkind == type(np.float32(1.)): |
---|
219 | attrvalue = np.float32(attrval) |
---|
220 | elif attrkind == type(np.float64(1.)): |
---|
221 | attrvalue = np.float64(attrval) |
---|
222 | elif attrkind == type(np.array([1.,2.])): |
---|
223 | attrvalue = attrval |
---|
224 | else: |
---|
225 | print errormsg |
---|
226 | print ' ' + fname + ": '" + attrkind + "' kind of attribute is not ready!" |
---|
227 | print ' valid values: _______' |
---|
228 | for key in validk.keys(): |
---|
229 | print ' ', key,':', validk[key] |
---|
230 | quit(-1) |
---|
231 | |
---|
232 | attvar = ncv.ncattrs() |
---|
233 | if searchInlist(attvar, attrname): |
---|
234 | attr = ncv.delncattr(attrname) |
---|
235 | |
---|
236 | if attrname == 'original_subroutines_author': |
---|
237 | attrvalue = 'Cindy Bruyere' |
---|
238 | attr = ncv.setncattr(attrname, attrvalue) |
---|
239 | |
---|
240 | return attr |
---|
241 | |
---|
242 | |
---|
243 | def basicvardef(varobj, vstname, vlname, vunits): |
---|
244 | """ Function to give the basic attributes to a variable |
---|
245 | varobj= netCDF variable object |
---|
246 | vstname= standard name of the variable |
---|
247 | vlname= long name of the variable |
---|
248 | vunits= units of the variable |
---|
249 | """ |
---|
250 | attr = varobj.setncattr('standard_name', vstname) |
---|
251 | attr = varobj.setncattr('long_name', vlname) |
---|
252 | attr = varobj.setncattr('units', vunits) |
---|
253 | |
---|
254 | return |
---|
255 | |
---|
256 | def remove_NONascii(string): |
---|
257 | """ Function to remove that characters which are not in the standard 127 ASCII |
---|
258 | string= string to transform |
---|
259 | >>> remove_NONascii('LluÃs') |
---|
260 | Lluis |
---|
261 | """ |
---|
262 | fname = 'remove_NONascii' |
---|
263 | |
---|
264 | newstring = string |
---|
265 | |
---|
266 | RTFchar= ['á', 'é', 'Ã', 'ó', 'ú', 'à ', 'Ú', 'ì', 'ò', 'ù', 'â', 'ê', 'î', 'ÃŽ', \ |
---|
267 | 'û', 'À', 'ë', 'ï', 'ö', 'ÃŒ', 'ç', 'ñ','Ê', 'Å', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', \ |
---|
268 | 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã', 'Ã',\ |
---|
269 | 'Ã', 'Å', '\n', '\t'] |
---|
270 | ASCchar= ['a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', 'u', 'a', 'e', 'i', 'o', \ |
---|
271 | 'u', 'a', 'e', 'i', 'o', 'u', 'c', 'n','ae', 'oe', 'A', 'E', 'I', 'O', 'U', 'A', \ |
---|
272 | 'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'A', 'E', 'I', 'O', 'U', 'C', 'N',\ |
---|
273 | 'AE', 'OE', '', ' '] |
---|
274 | |
---|
275 | Nchars = len(RTFchar) |
---|
276 | for ichar in range(Nchars): |
---|
277 | foundchar = string.find(RTFchar[ichar]) |
---|
278 | if foundchar != -1: |
---|
279 | newstring = newstring.replace(RTFchar[ichar], ASCchar[ichar]) |
---|
280 | |
---|
281 | return newstring |
---|
282 | |
---|
283 | def read_description(fdobs, dbg): |
---|
284 | """ reads the description file of the observational data-set |
---|
285 | fdobs= descriptive observational data-set |
---|
286 | dbg= boolean argument for debugging |
---|
287 | * Each station should have a 'description.dat' file with: |
---|
288 | institution=Institution who creates the data |
---|
289 | department=Department within the institution |
---|
290 | scientists=names of the data producers |
---|
291 | contact=contact of the data producers |
---|
292 | description=description of the observations |
---|
293 | acknowledgement=sentence of acknowlegement |
---|
294 | comment=comment for the measurements |
---|
295 | |
---|
296 | MissingValue='|' list of ASCII values for missing values within the data |
---|
297 | (as they appear!, 'empty' for no value at all) |
---|
298 | comment=comments |
---|
299 | |
---|
300 | varN='|' list of variable names |
---|
301 | varLN='|' list of long variable names |
---|
302 | varU='|' list units of the variables |
---|
303 | varBUFR='|' list BUFR code of the variables |
---|
304 | varTYPE='|' list of variable types ('D', 'F', 'I', 'I64', 'S') |
---|
305 | varOPER='|' list of operations to do to the variables to meet their units ([oper],[val]) |
---|
306 | [oper]: |
---|
307 | -, nothing |
---|
308 | sumc, add [val] |
---|
309 | subc, rest [val] |
---|
310 | mulc, multiply by [val] |
---|
311 | divc, divide by [val] |
---|
312 | rmchar,[val],[pos], remove [val] characters from [pos]='B', beginning, 'E', end |
---|
313 | NAMElon=name of the variable with the longitude (x position) |
---|
314 | NAMElat=name of the variable with the latitude (y position) |
---|
315 | NAMEheight=ind_alt |
---|
316 | NAMEtime=name of the varibale with the time |
---|
317 | FMTtime=format of the time (as in 'C', 'CFtime' for already CF-like time) |
---|
318 | |
---|
319 | """ |
---|
320 | fname = 'read_description' |
---|
321 | |
---|
322 | descobj = open(fdobs, 'r') |
---|
323 | desc = {} |
---|
324 | |
---|
325 | namevalues = [] |
---|
326 | |
---|
327 | for line in descobj: |
---|
328 | if line[0:1] != '#' and len(line) > 1 : |
---|
329 | descn = remove_NONascii(line.split('=')[0]) |
---|
330 | descv = remove_NONascii(line.split('=')[1]) |
---|
331 | namevalues.append(descn) |
---|
332 | if descn[0:3] != 'var': |
---|
333 | if descn != 'MissingValue': |
---|
334 | desc[descn] = descv |
---|
335 | else: |
---|
336 | desc[descn] = [] |
---|
337 | for dn in descv.split('|'): desc[descn].append(dn) |
---|
338 | print ' ' + fname + ': missing values found:',desc[descn] |
---|
339 | elif descn[0:3] == 'var': |
---|
340 | desc[descn] = descv.split('|') |
---|
341 | elif descn[0:4] == 'NAME': |
---|
342 | desc[descn] = descv |
---|
343 | elif descn[0:3] == 'FMT': |
---|
344 | desc[descn] = descv |
---|
345 | if not desc.has_key('varOPER'): desc['varOPER'] = None |
---|
346 | |
---|
347 | if dbg: |
---|
348 | Nvars = len(desc['varN']) |
---|
349 | print ' ' + fname + ": description content of '" + fdobs + "'______________" |
---|
350 | for varn in namevalues: |
---|
351 | if varn[0:3] != 'var': |
---|
352 | print ' ' + varn + ':',desc[varn] |
---|
353 | elif varn == 'varN': |
---|
354 | print ' * Variables:' |
---|
355 | for ivar in range(Nvars): |
---|
356 | varname = desc['varN'][ivar] |
---|
357 | varLname = desc['varLN'][ivar] |
---|
358 | varunits = desc['varU'][ivar] |
---|
359 | if desc.has_key('varBUFR'): |
---|
360 | varbufr = desc['varBUFR'][ivar] |
---|
361 | else: |
---|
362 | varbufr = None |
---|
363 | if desc['varOPER'] is not None: |
---|
364 | opv = desc['varOPER'][ivar] |
---|
365 | else: |
---|
366 | opv = None |
---|
367 | print ' ', ivar, varname+':',varLname,'[',varunits, \ |
---|
368 | ']','bufr code:',varbufr,'oper:',opv |
---|
369 | |
---|
370 | descobj.close() |
---|
371 | |
---|
372 | return desc |
---|
373 | |
---|
374 | def value_fmt(val, miss, op, fmt): |
---|
375 | """ Function to transform an ASCII value to a given format |
---|
376 | val= value to transform |
---|
377 | miss= list of possible missing values |
---|
378 | op= operation to perform to the value |
---|
379 | fmt= format to take: |
---|
380 | 'D': float double precission |
---|
381 | 'F': float |
---|
382 | 'I': integer |
---|
383 | 'I64': 64-bits integer |
---|
384 | 'S': string |
---|
385 | >>> value_fmt('9876.12', '-999', 'F') |
---|
386 | 9876.12 |
---|
387 | """ |
---|
388 | |
---|
389 | fname = 'value_fmt' |
---|
390 | |
---|
391 | aopers = ['sumc','subc','mulc','divc', 'rmchar'] |
---|
392 | |
---|
393 | fmts = ['D', 'F', 'I', 'I64', 'S'] |
---|
394 | Nfmts = len(fmts) |
---|
395 | |
---|
396 | if not searchInlist(miss,val): |
---|
397 | if searchInlist(miss,'empty') and len(val) == 0: |
---|
398 | newval = None |
---|
399 | else: |
---|
400 | if op != '-': |
---|
401 | opern = op.split(',')[0] |
---|
402 | operv = op.split(',')[1] |
---|
403 | |
---|
404 | if not searchInlist(aopers,opern): |
---|
405 | print errormsg |
---|
406 | print ' ' + fname + ": operation '" + opern + "' not ready!!" |
---|
407 | print ' availables:',aopers |
---|
408 | quit(-1) |
---|
409 | else: |
---|
410 | opern = 'sumc' |
---|
411 | operv = '0' |
---|
412 | |
---|
413 | if not searchInlist(fmts, fmt): |
---|
414 | print errormsg |
---|
415 | print ' ' + fname + ": format '" + fmt + "' not ready !!" |
---|
416 | quit(-1) |
---|
417 | else: |
---|
418 | if fmt == 'D': |
---|
419 | opv = np.float32(operv) |
---|
420 | if opern == 'sumc': newval = np.float32(val) + opv |
---|
421 | elif opern == 'subc': newval = np.float32(val) - opv |
---|
422 | elif opern == 'mulc': newval = np.float32(val) * opv |
---|
423 | elif opern == 'divc': newval = np.float32(val) / opv |
---|
424 | elif opern == 'rmchar': |
---|
425 | opv = int(operv) |
---|
426 | Lval = len(val) |
---|
427 | if op.split(',')[2] == 'B': |
---|
428 | newval = np.float32(val[opv:Lval+1]) |
---|
429 | elif op.split(',')[2] == 'E': |
---|
430 | newval = np.float32(val[Lval-opv:Lval]) |
---|
431 | else: |
---|
432 | print errormsg |
---|
433 | print ' ' + fname + ": operation '" + opern + "' not " +\ |
---|
434 | " work with '" + op.split(',')[2] + "' !!" |
---|
435 | quit(-1) |
---|
436 | elif fmt == 'F': |
---|
437 | opv = np.float(operv) |
---|
438 | if opern == 'sumc': newval = np.float(val) + opv |
---|
439 | elif opern == 'subc': newval = np.float(val) - opv |
---|
440 | elif opern == 'mulc': newval = np.float(val) * opv |
---|
441 | elif opern == 'divc': newval = np.float(val) / opv |
---|
442 | elif opern == 'rmchar': |
---|
443 | opv = int(operv) |
---|
444 | Lval = len(val) |
---|
445 | if op.split(',')[2] == 'B': |
---|
446 | newval = np.float(val[opv:Lval+1]) |
---|
447 | elif op.split(',')[2] == 'E': |
---|
448 | newval = np.float(val[0:Lval-opv]) |
---|
449 | else: |
---|
450 | print errormsg |
---|
451 | print ' ' + fname + ": operation '" + opern + "' not " +\ |
---|
452 | " work with '" + op.split(',')[2] + "' !!" |
---|
453 | quit(-1) |
---|
454 | |
---|
455 | elif fmt == 'I': |
---|
456 | opv = int(operv) |
---|
457 | if opern == 'sumc': newval = int(val) + opv |
---|
458 | elif opern == 'subc': newval = int(val) - opv |
---|
459 | elif opern == 'mulc': newval = int(val) * opv |
---|
460 | elif opern == 'divc': newval = int(val) / opv |
---|
461 | elif opern == 'rmchar': |
---|
462 | opv = int(operv) |
---|
463 | Lval = len(val) |
---|
464 | if op.split(',')[2] == 'B': |
---|
465 | newval = int(val[opv:Lval+1]) |
---|
466 | elif op.split(',')[2] == 'E': |
---|
467 | newval = int(val[0:Lval-opv]) |
---|
468 | else: |
---|
469 | print errormsg |
---|
470 | print ' ' + fname + ": operation '" + opern + "' not " +\ |
---|
471 | " work with '" + op.split(',')[2] + "' !!" |
---|
472 | quit(-1) |
---|
473 | elif fmt == 'I64': |
---|
474 | opv = np.int64(operv) |
---|
475 | if opern == 'sumc': newval = np.int64(val) + opv |
---|
476 | elif opern == 'subc': newval = np.int64(val) - opv |
---|
477 | elif opern == 'mulc': newval = np.int64(val) * opv |
---|
478 | elif opern == 'divc': newval = np.int64(val) / opv |
---|
479 | elif opern == 'rmchar': |
---|
480 | opv = int(operv) |
---|
481 | Lval = len(val) |
---|
482 | if op.split(',')[2] == 'B': |
---|
483 | newval = np.int64(val[opv:Lval+1]) |
---|
484 | elif op.split(',')[2] == 'E': |
---|
485 | newval = np.int64(val[0:Lval-opv]) |
---|
486 | else: |
---|
487 | print errormsg |
---|
488 | print ' ' + fname + ": operation '" + opern + "' not " +\ |
---|
489 | " work with '" + op.split(',')[2] + "' !!" |
---|
490 | quit(-1) |
---|
491 | elif fmt == 'S': |
---|
492 | if opern == 'rmchar': |
---|
493 | opv = int(operv) |
---|
494 | Lval = len(val) |
---|
495 | if op.split(',')[2] == 'B': |
---|
496 | newval = val[opv:Lval+1] |
---|
497 | elif op.split(',')[2] == 'E': |
---|
498 | newval = val[0:Lval-opv] |
---|
499 | else: |
---|
500 | print errormsg |
---|
501 | print ' ' + fname + ": operation '" + opern + "' not " +\ |
---|
502 | " work with '" + op.split(',')[2] + "' !!" |
---|
503 | quit(-1) |
---|
504 | else: |
---|
505 | newval = val |
---|
506 | |
---|
507 | else: |
---|
508 | newval = None |
---|
509 | |
---|
510 | return newval |
---|
511 | |
---|
512 | def getting_fixedline(line, cuts, types, dbg=False): |
---|
513 | """ Function to get the values from a line of text with fixed lenght of different values |
---|
514 | line: line with values |
---|
515 | cuts: character number where a value ends |
---|
516 | types: consecutive type of values |
---|
517 | 'I': integer |
---|
518 | 'R': real |
---|
519 | 'D': float64 |
---|
520 | 'S': string |
---|
521 | 'B': boolean |
---|
522 | dbg: debug mode (default False) |
---|
523 | >>> Sline=' 87007 03012015 25.6 6.4 9.4 5 15' |
---|
524 | >>> getting_fixedline(Sline, [8, 17, 23, 29, 36, 40, 45, 50], ['I', 'R', 'R', 'R', 'I', 'R', 'R', 'R']) |
---|
525 | [87007, 3012015.0, 25.6, 6.4, -99999, -99999, 9.4, 5.0, 15.0] |
---|
526 | """ |
---|
527 | fname = 'getting_fixedline' |
---|
528 | |
---|
529 | if len(cuts) + 1 != len(types): |
---|
530 | print errormsg |
---|
531 | print ' ' + fname + ': The number of types :', len(types), 'must be +1', \ |
---|
532 | 'number of cuts:', len(cuts) |
---|
533 | quit(-1) |
---|
534 | |
---|
535 | values = [] |
---|
536 | val = line[0:cuts[0]] |
---|
537 | if len(val.replace(' ','')) >= 1: |
---|
538 | values.append(typemod(val, types[0])) |
---|
539 | else: |
---|
540 | if types[0] == 'I': values.append(fillValueI) |
---|
541 | elif types[0] == 'R': values.append(fillValueF) |
---|
542 | elif types[0] == 'D': values.append(fillValueF) |
---|
543 | elif types[0] == 'S': values.append(fillValueS) |
---|
544 | elif types[0] == 'B': values.append(fillValueB) |
---|
545 | |
---|
546 | Ncuts = len(cuts) |
---|
547 | for ic in range(1,Ncuts): |
---|
548 | val = line[cuts[ic-1]:cuts[ic]] |
---|
549 | if dbg: print ic, ':', val, '-->', types[ic] |
---|
550 | if len(val.replace(' ','')) >= 1: |
---|
551 | values.append(typemod(val, types[ic])) |
---|
552 | else: |
---|
553 | if types[ic] == 'I': values.append(fillValueI) |
---|
554 | elif types[ic] == 'R': values.append(fillValueF) |
---|
555 | elif types[ic] == 'D': values.append(fillValueF) |
---|
556 | elif types[ic] == 'S': values.append(fillValueS) |
---|
557 | elif types[ic] == 'B': values.append(fillValueB) |
---|
558 | |
---|
559 | # Last value |
---|
560 | Lline = len(line) |
---|
561 | val = line[cuts[Ncuts-1]:Lline] |
---|
562 | if len(val.replace(' ','')) >= 1: |
---|
563 | values.append(typemod(val, types[Ncuts])) |
---|
564 | else: |
---|
565 | if types[Ncuts] == 'I': values.append(fillValueI) |
---|
566 | elif types[Ncuts] == 'R': values.append(fillValueF) |
---|
567 | elif types[Ncuts] == 'D': values.append(fillValueF) |
---|
568 | elif types[Ncuts] == 'S': values.append(fillValueS) |
---|
569 | elif types[Ncuts] == 'B': values.append(fillValueB) |
---|
570 | |
---|
571 | return values |
---|
572 | |
---|
573 | |
---|
574 | def getting_fixedline_NOk(line, cuts, miss, dbg=False): |
---|
575 | """ Function to get the values from a line of text with fixed lenght of |
---|
576 | different values without types |
---|
577 | line: line with values |
---|
578 | cuts: character number where a value ends |
---|
579 | miss: character form issed values |
---|
580 | dbg: debug mode (default False) |
---|
581 | >>> Sline=' 87007 03012015 25.6 6.4 9.4 5 15' |
---|
582 | >>> getting_fixedline_NOk(Sline, [8, 17, 23, 29, 36, 40, 45, 50], '#') |
---|
583 | ['87007', '03012015', '25.6', '6.4', '#', '#', '9.4', '5', '15'] |
---|
584 | """ |
---|
585 | fname = 'getting_fixedline_NOk' |
---|
586 | |
---|
587 | values = [] |
---|
588 | val = line[0:cuts[0]] |
---|
589 | if len(val.replace(' ','')) >= 1: |
---|
590 | values.append(val.replace(' ','')) |
---|
591 | else: |
---|
592 | values.append(miss) |
---|
593 | |
---|
594 | Ncuts = len(cuts) |
---|
595 | for ic in range(1,Ncuts): |
---|
596 | val = line[cuts[ic-1]:cuts[ic]] |
---|
597 | if dbg: print ic, ':', val |
---|
598 | if len(val.replace(' ','')) >= 1: |
---|
599 | values.append(val.replace(' ','')) |
---|
600 | else: |
---|
601 | values.append(miss) |
---|
602 | |
---|
603 | # Last value |
---|
604 | Lline = len(line) |
---|
605 | val = line[cuts[Ncuts-1]:Lline] |
---|
606 | if len(val.replace(' ','')) >= 1: |
---|
607 | values.append(val.replace(' ','')) |
---|
608 | else: |
---|
609 | values.append(miss) |
---|
610 | |
---|
611 | return values |
---|
612 | |
---|
613 | def read_datavalues(dataf, comchar, colchar, fmt, jBl, oper, miss, varns, dbg): |
---|
614 | """ Function to read from an ASCII file values in column |
---|
615 | dataf= data file |
---|
616 | comchar= list of the characters indicating comment in the file |
---|
617 | colchar= character which indicate end of value in the column |
---|
618 | dbg= debug mode or not |
---|
619 | fmt= list of kind of values to be found |
---|
620 | jBl= number of lines to jump from the beginning of file |
---|
621 | oper= list of operations to perform |
---|
622 | miss= missing value |
---|
623 | varns= list of name of the variables to find |
---|
624 | """ |
---|
625 | |
---|
626 | fname = 'read_datavalues' |
---|
627 | |
---|
628 | ofile = open(dataf, 'r') |
---|
629 | Nvals = len(fmt) |
---|
630 | |
---|
631 | if oper is None: |
---|
632 | opers = [] |
---|
633 | for ioper in range(Nvals): |
---|
634 | opers.append('-') |
---|
635 | else: |
---|
636 | opers = oper |
---|
637 | |
---|
638 | finalvalues = {} |
---|
639 | |
---|
640 | iline = 0 |
---|
641 | for line in ofile: |
---|
642 | line = line.replace('\n','').replace(chr(13),'') |
---|
643 | if not searchInlist(comchar,line[0:1]) and len(line) > 1 and iline > jBl-1: |
---|
644 | # Getting values |
---|
645 | if colchar[0:4] != 'None': |
---|
646 | values0 = line.split(colchar) |
---|
647 | else: |
---|
648 | ltypes=str_list_k(colchar[5:len(colchar)+1],',','I') |
---|
649 | values0 = getting_fixedline_NOk(line, ltypes, miss[0], dbg=dbg) |
---|
650 | # Removing no-value columns |
---|
651 | values = [] |
---|
652 | for iv in values0: |
---|
653 | if len(iv) > 0: values.append(iv) |
---|
654 | |
---|
655 | Nvalues = len(values) |
---|
656 | # Checkings for wierd characters at the end of lines (use it to check) |
---|
657 | # if values[Nvalues-1][0:4] == '-999': |
---|
658 | # print line,'last v:',values[Nvalues-1],'len:',len(values[Nvalues-1]) |
---|
659 | # for ic in range(len(values[Nvalues-1])): |
---|
660 | # print ic,ord(values[Nvalues-1][ic:ic+1]) |
---|
661 | # quit() |
---|
662 | |
---|
663 | if len(values[Nvalues-1]) == 0: |
---|
664 | Nvalues = Nvalues - 1 |
---|
665 | |
---|
666 | if Nvalues != Nvals: |
---|
667 | print warnmsg |
---|
668 | print ' ' + fname + ': number of formats:',Nvals,' and number of ', \ |
---|
669 | 'values:',Nvalues,' with split character *' + colchar + \ |
---|
670 | '* does not coincide!!' |
---|
671 | print ' * what is found is ________' |
---|
672 | if Nvalues > Nvals: |
---|
673 | Nshow = Nvals |
---|
674 | for ivar in range(Nshow): |
---|
675 | print ' ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar] |
---|
676 | print ' missing formats for:',values[Nshow:Nvalues+1] |
---|
677 | print ' values not considered, continue' |
---|
678 | else: |
---|
679 | Nshow = Nvalues |
---|
680 | for ivar in range(Nshow): |
---|
681 | print ' ',varns[ivar],'fmt:',fmt[ivar],'value:',values[ivar] |
---|
682 | print ' excess of formats:',fmt[Nshow:Nvals+1] |
---|
683 | quit(-1) |
---|
684 | |
---|
685 | # Reading and transforming values |
---|
686 | if dbg: print ' ' + fname + ': values found _________' |
---|
687 | |
---|
688 | for ivar in range(Nvals): |
---|
689 | if dbg: |
---|
690 | print iline, ',', ivar, ' ', varns[ivar],'value:',values[ivar], \ |
---|
691 | miss,opers[ivar], fmt[ivar] |
---|
692 | |
---|
693 | if iline == 0: |
---|
694 | listvals = [] |
---|
695 | listvals.append(value_fmt(values[ivar], miss, opers[ivar], \ |
---|
696 | fmt[ivar])) |
---|
697 | finalvalues[varns[ivar]] = listvals |
---|
698 | else: |
---|
699 | listvals = finalvalues[varns[ivar]] |
---|
700 | listvals.append(value_fmt(values[ivar], miss, opers[ivar], \ |
---|
701 | fmt[ivar])) |
---|
702 | finalvalues[varns[ivar]] = listvals |
---|
703 | else: |
---|
704 | # First line without values |
---|
705 | if iline == 0: iline = -1 |
---|
706 | |
---|
707 | iline = iline + 1 |
---|
708 | |
---|
709 | ofile.close() |
---|
710 | |
---|
711 | return finalvalues |
---|
712 | |
---|
713 | def writing_str_nc(varo, values, Lchar): |
---|
714 | """ Function to write string values in a netCDF variable as a chain of 1char values |
---|
715 | varo= netCDF variable object |
---|
716 | values = list of values to introduce |
---|
717 | Lchar = length of the string in the netCDF file |
---|
718 | """ |
---|
719 | |
---|
720 | Nvals = len(values) |
---|
721 | |
---|
722 | for iv in range(Nvals): |
---|
723 | stringv=values[iv] |
---|
724 | charvals = np.chararray(Lchar) |
---|
725 | Lstr = len(stringv) |
---|
726 | charvals[Lstr:Lchar] = '' |
---|
727 | |
---|
728 | for ich in range(Lstr): |
---|
729 | charvals[ich] = stringv[ich:ich+1] |
---|
730 | |
---|
731 | varo[iv,:] = charvals |
---|
732 | |
---|
733 | return |
---|
734 | |
---|
735 | def Stringtimes_CF(tvals, fmt, Srefdate, tunits, dbg): |
---|
736 | """ Function to transform a given data in String formats to a CF date |
---|
737 | tvals= string temporal values |
---|
738 | fmt= format of the the time values |
---|
739 | Srefdate= reference date in [YYYY][MM][DD][HH][MI][SS] format |
---|
740 | tunits= units to use ('weeks', 'days', 'hours', 'minutes', 'seconds') |
---|
741 | dbg= debug |
---|
742 | >>> Stringtimes_CF(['19760217082712','20150213101837'], '%Y%m%d%H%M%S', |
---|
743 | '19491201000000', 'hours', False) |
---|
744 | [229784.45333333 571570.31027778] |
---|
745 | """ |
---|
746 | import datetime as dt |
---|
747 | |
---|
748 | fname = 'Stringtimes' |
---|
749 | |
---|
750 | dimt = len(tvals) |
---|
751 | |
---|
752 | yrref = int(Srefdate[0:4]) |
---|
753 | monref = int(Srefdate[4:6]) |
---|
754 | dayref = int(Srefdate[6:8]) |
---|
755 | horref = int(Srefdate[8:10]) |
---|
756 | minref = int(Srefdate[10:12]) |
---|
757 | secref = int(Srefdate[12:14]) |
---|
758 | refdate = dt.datetime( yrref, monref, dayref, horref, minref, secref) |
---|
759 | |
---|
760 | cftimes = np.zeros((dimt), dtype=np.float) |
---|
761 | |
---|
762 | Nfmt=len(fmt.split('%')) |
---|
763 | |
---|
764 | if dbg: print ' ' + fname + ': fmt=',fmt,'refdate:',Srefdate,'uits:',tunits, \ |
---|
765 | 'date dt_days dt_time deltaseconds CFtime _______' |
---|
766 | for it in range(dimt): |
---|
767 | |
---|
768 | # Removing excess of mili-seconds (up to 6 decimals) |
---|
769 | if fmt.split('%')[Nfmt-1] == 'f': |
---|
770 | tpoints = tvals[it].split('.') |
---|
771 | if len(tpoints[len(tpoints)-1]) > 6: |
---|
772 | milisec = '{0:.6f}'.format(np.float('0.'+tpoints[len(tpoints)-1]))[0:7] |
---|
773 | newtval = '' |
---|
774 | for ipt in range(len(tpoints)-1): |
---|
775 | newtval = newtval + tpoints[ipt] + '.' |
---|
776 | newtval = newtval + str(milisec)[2:len(milisec)+1] |
---|
777 | else: |
---|
778 | newtval = tvals[it] |
---|
779 | tval = dt.datetime.strptime(newtval, fmt) |
---|
780 | else: |
---|
781 | tval = dt.datetime.strptime(tvals[it], fmt) |
---|
782 | |
---|
783 | deltatime = tval - refdate |
---|
784 | deltaseconds = deltatime.days*24.*3600. + deltatime.seconds + \ |
---|
785 | deltatime.microseconds/100000. |
---|
786 | if tunits == 'weeks': |
---|
787 | deltat = 7.*24.*3600. |
---|
788 | elif tunits == 'days': |
---|
789 | deltat = 24.*3600. |
---|
790 | elif tunits == 'hours': |
---|
791 | deltat = 3600. |
---|
792 | elif tunits == 'minutes': |
---|
793 | deltat = 60. |
---|
794 | elif tunits == 'seconds': |
---|
795 | deltat = 1. |
---|
796 | else: |
---|
797 | print errormsg |
---|
798 | print ' ' + fname + ": time units '" + tunits + "' not ready !!" |
---|
799 | quit(-1) |
---|
800 | |
---|
801 | cftimes[it] = deltaseconds / deltat |
---|
802 | if dbg: |
---|
803 | print ' ' + tvals[it], deltatime, deltaseconds, cftimes[it] |
---|
804 | |
---|
805 | return cftimes |
---|
806 | |
---|
807 | def adding_complementary(onc, dscn, okind, dvalues, tvals, refCFt, CFtu, Nd, dbg): |
---|
808 | """ Function to add complementary variables as function of the observational type |
---|
809 | onc= netcdf objecto file to add the variables |
---|
810 | dscn= description dictionary |
---|
811 | okind= observational kind |
---|
812 | dvalues= values |
---|
813 | tvals= CF time values |
---|
814 | refCFt= reference time of CF time (in [YYYY][MM][DD][HH][MI][SS]) |
---|
815 | CFtu= CF time units |
---|
816 | Nd= size of the domain |
---|
817 | dbg= debugging flag |
---|
818 | """ |
---|
819 | import numpy.ma as ma |
---|
820 | |
---|
821 | fname = 'adding_complementary' |
---|
822 | |
---|
823 | # Kind of observations which require de integer lon/lat (for the 2D map) |
---|
824 | map2D=['multi-points', 'trajectory'] |
---|
825 | |
---|
826 | SrefCFt = refCFt[0:4] +'-'+ refCFt[4:6] +'-'+ refCFt[6:8] + ' ' + refCFt[8:10] + \ |
---|
827 | ':'+ refCFt[10:12] +':'+ refCFt[12:14] |
---|
828 | |
---|
829 | if dscn['NAMElon'] == '-' or dscn['NAMElat'] == '-': |
---|
830 | print errormsg |
---|
831 | print ' ' + fname + ": to complement a '" + okind + "' observation kind " + \ |
---|
832 | " a given longitude ('NAMElon':",dscn['NAMElon'],") and latitude ('" + \ |
---|
833 | "'NAMElat:'", dscn['NAMElat'],') from the data has to be provided and ' + \ |
---|
834 | 'any are given !!' |
---|
835 | quit(-1) |
---|
836 | |
---|
837 | if not dvalues.has_key(dscn['NAMElon']) or not dvalues.has_key(dscn['NAMElat']): |
---|
838 | print errormsg |
---|
839 | print ' ' + fname + ": observations do not have 'NAMElon':", \ |
---|
840 | dscn['NAMElon'],"and/or 'NAMElat:'", dscn['NAMElat'],' !!' |
---|
841 | print ' available data:',dvalues.keys() |
---|
842 | quit(-1) |
---|
843 | |
---|
844 | if okind == 'trajectory': |
---|
845 | if dscn['NAMEheight'] == '-': |
---|
846 | print warnmsg |
---|
847 | print ' ' + fname + ": to complement a '" + okind + "' observation " + \ |
---|
848 | "kind a given height ('NAMEheight':",dscn['NAMEheight'],"') might " + \ |
---|
849 | 'be provided and any is given !!' |
---|
850 | quit(-1) |
---|
851 | |
---|
852 | if not dvalues.has_key(dscn['NAMEheight']): |
---|
853 | print errormsg |
---|
854 | print ' ' + fname + ": observations do not have 'NAMEtime':", \ |
---|
855 | dscn['NAMEtime'],' !!' |
---|
856 | print ' available data:',dvalues.keys() |
---|
857 | quit(-1) |
---|
858 | |
---|
859 | if searchInlist(map2D, okind): |
---|
860 | # A new 2D map with the number of observation will be added for that 'NAMElon' |
---|
861 | # and 'NAMElat' are necessary. A NdxNd domain space size will be used. |
---|
862 | objfile.createDimension('lon2D',Nd) |
---|
863 | objfile.createDimension('lat2D',Nd) |
---|
864 | lonvals = ma.masked_equal(dvalues[dscn['NAMElon']], None) |
---|
865 | latvals = ma.masked_equal(dvalues[dscn['NAMElat']], None) |
---|
866 | |
---|
867 | minlon = min(lonvals) |
---|
868 | maxlon = max(lonvals) |
---|
869 | minlat = min(latvals) |
---|
870 | maxlat = max(latvals) |
---|
871 | |
---|
872 | blon = (maxlon - minlon)/(Nd-1) |
---|
873 | blat = (maxlat - minlat)/(Nd-1) |
---|
874 | |
---|
875 | newvar = onc.createVariable( 'lon2D', 'f8', ('lon2D')) |
---|
876 | basicvardef(newvar, 'longitude', 'longitude map observations','degrees_East') |
---|
877 | newvar[:] = minlon + np.arange(Nd)*blon |
---|
878 | newattr = set_attribute(newvar, 'axis', 'X') |
---|
879 | |
---|
880 | newvar = onc.createVariable( 'lat2D', 'f8', ('lat2D')) |
---|
881 | basicvardef(newvar, 'latitude', 'latitude map observations', 'degrees_North') |
---|
882 | newvar[:] = minlat + np.arange(Nd)*blat |
---|
883 | newattr = set_attribute(newvar, 'axis', 'Y') |
---|
884 | |
---|
885 | if dbg: |
---|
886 | print ' ' + fname + ': minlon=',minlon,'maxlon=',maxlon |
---|
887 | print ' ' + fname + ': minlat=',minlat,'maxlat=',maxlat |
---|
888 | print ' ' + fname + ': precission on x-axis=', blon*(Nd-1), 'y-axis=', \ |
---|
889 | blat*(Nd-1) |
---|
890 | |
---|
891 | if okind == 'multi-points': |
---|
892 | map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI |
---|
893 | |
---|
894 | dimt = len(tvals) |
---|
895 | Nlost = 0 |
---|
896 | for it in range(dimt): |
---|
897 | lon = dvalues[dscn['NAMElon']][it] |
---|
898 | lat = dvalues[dscn['NAMElat']][it] |
---|
899 | if lon is not None and lat is not None: |
---|
900 | ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon)) |
---|
901 | ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat)) |
---|
902 | |
---|
903 | if map2D[ilat,ilon] == fillValueI: |
---|
904 | map2D[ilat,ilon] = 1 |
---|
905 | else: |
---|
906 | map2D[ilat,ilon] = map2D[ilat,ilon] + 1 |
---|
907 | if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon] |
---|
908 | |
---|
909 | newvar = onc.createVariable( 'mapobs', 'f4', ('lat2D', 'lon2D'), \ |
---|
910 | fill_value = fillValueI) |
---|
911 | basicvardef(newvar, 'map_observations', 'number of observations', '-') |
---|
912 | newvar[:] = map2D |
---|
913 | newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D') |
---|
914 | |
---|
915 | elif okind == 'trajectory': |
---|
916 | # A new 2D map with the trajectory 'NAMElon' and 'NAMElat' and maybe 'NAMEheight' |
---|
917 | # are necessary. A NdxNdxNd domain space size will be used. Using time as |
---|
918 | # reference variable |
---|
919 | if dscn['NAMEheight'] == '-': |
---|
920 | # No height |
---|
921 | map2D = np.ones((Nd, Nd), dtype=np.float)*fillValueI |
---|
922 | |
---|
923 | dimt = len(tvals) |
---|
924 | Nlost = 0 |
---|
925 | if dbg: print ' time-step lon lat ix iy passes _______' |
---|
926 | for it in range(dimt): |
---|
927 | lon = dvalues[dscn['NAMElon']][it] |
---|
928 | lat = dvalues[dscn['NAMElat']][it] |
---|
929 | if lon is not None and lat is not None: |
---|
930 | ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon)) |
---|
931 | ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat)) |
---|
932 | |
---|
933 | if map2D[ilat,ilon] == fillValueI: |
---|
934 | map2D[ilat,ilon] = 1 |
---|
935 | else: |
---|
936 | map2D[ilat,ilon] = map2D[ilat,ilon] + 1 |
---|
937 | if dbg: print it, lon, lat, ilon, ilat, map2D[ilat,ilon] |
---|
938 | |
---|
939 | newvar = onc.createVariable( 'trjobs', 'i', ('lat2D', 'lon2D'), \ |
---|
940 | fill_value = fillValueI) |
---|
941 | basicvardef(newvar, 'trajectory_observations', 'number of passes', '-' ) |
---|
942 | newvar[:] = map2D |
---|
943 | newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D') |
---|
944 | |
---|
945 | else: |
---|
946 | ivn = 0 |
---|
947 | for vn in dscn['varN']: |
---|
948 | if vn == dscn['NAMEheight']: |
---|
949 | zu = dscn['varU'][ivn] |
---|
950 | break |
---|
951 | ivn = ivn + 1 |
---|
952 | |
---|
953 | objfile.createDimension('z2D',Nd) |
---|
954 | zvals = ma.masked_equal(dvalues[dscn['NAMEheight']], None) |
---|
955 | minz = min(zvals) |
---|
956 | maxz = max(zvals) |
---|
957 | |
---|
958 | bz = (maxz - minz)/(Nd-1) |
---|
959 | |
---|
960 | newvar = onc.createVariable( 'z2D', 'f8', ('z2D')) |
---|
961 | basicvardef(newvar, 'z2D', 'z-coordinate map observations', zu) |
---|
962 | newvar[:] = minz + np.arange(Nd)*bz |
---|
963 | newattr = set_attribute(newvar, 'axis', 'Z') |
---|
964 | |
---|
965 | if dbg: |
---|
966 | print ' ' + fname + ': zmin=',minz,zu,'zmax=',maxz,zu |
---|
967 | print ' ' + fname + ': precission on z-axis=', bz*(Nd-1), zu |
---|
968 | |
---|
969 | map3D = np.ones((Nd, Nd, Nd), dtype=int)*fillValueI |
---|
970 | dimt = len(tvals) |
---|
971 | Nlost = 0 |
---|
972 | if dbg: print ' time-step lon lat z ix iy iz passes _______' |
---|
973 | for it in range(dimt): |
---|
974 | lon = dvalues[dscn['NAMElon']][it] |
---|
975 | lat = dvalues[dscn['NAMElat']][it] |
---|
976 | z = dvalues[dscn['NAMEheight']][it] |
---|
977 | if lon is not None and lat is not None and z is not None: |
---|
978 | ilon = int((Nd-1)*(lon - minlon)/(maxlon - minlon)) |
---|
979 | ilat = int((Nd-1)*(lat - minlat)/(maxlat - minlat)) |
---|
980 | iz = int((Nd-1)*(z - minz)/(maxz - minz)) |
---|
981 | |
---|
982 | if map3D[iz,ilat,ilon] == fillValueI: |
---|
983 | map3D[iz,ilat,ilon] = 1 |
---|
984 | else: |
---|
985 | map3D[iz,ilat,ilon] = map3D[iz,ilat,ilon] + 1 |
---|
986 | if dbg: print it, lon, lat, z, ilon, ilat, iz, map3D[iz,ilat,ilon] |
---|
987 | |
---|
988 | newvar = onc.createVariable( 'trjobs', 'i', ('z2D', 'lat2D', 'lon2D'), \ |
---|
989 | fill_value = fillValueI) |
---|
990 | basicvardef(newvar, 'trajectory_observations', 'number of passes', '-') |
---|
991 | newvar[:] = map3D |
---|
992 | newattr = set_attribute(newvar, 'coordinates', 'lon2D lat2D z2D') |
---|
993 | |
---|
994 | onc.sync() |
---|
995 | return |
---|
996 | |
---|
997 | def adding_station_desc(onc,stdesc): |
---|
998 | """ Function to add a station description in a netCDF file |
---|
999 | onc= netCDF object |
---|
1000 | stdesc= station description name, lon, lat, height |
---|
1001 | """ |
---|
1002 | fname = 'adding_station_desc' |
---|
1003 | |
---|
1004 | newdim = onc.createDimension('nst',1) |
---|
1005 | |
---|
1006 | newvar = objfile.createVariable( 'station', 'c', ('nst','StrLength')) |
---|
1007 | writing_str_nc(newvar, [stdesc[0].replace('!', ' ')], StringLength) |
---|
1008 | |
---|
1009 | newvar = objfile.createVariable( 'lonstGDM', 'c', ('nst','StrLength')) |
---|
1010 | Gv = int(stdesc[1]) |
---|
1011 | Dv = int((stdesc[1] - Gv)*60.) |
---|
1012 | Mv = int((stdesc[1] - Gv - Dv/60.)*3600.) |
---|
1013 | writing_str_nc(newvar, [str(Gv)+"d" + str(Dv)+"m" + str(Mv)+'s'], StringLength) |
---|
1014 | |
---|
1015 | if onc.variables.has_key('lon'): |
---|
1016 | print warnmsg |
---|
1017 | print ' ' + fname + ": variable 'lon' already exist !!" |
---|
1018 | print " renaming it as 'lonst'" |
---|
1019 | lonname = 'lonst' |
---|
1020 | else: |
---|
1021 | lonname = 'lon' |
---|
1022 | |
---|
1023 | newvar = objfile.createVariable( lonname, 'f4', ('nst')) |
---|
1024 | basicvardef(newvar, lonname, 'longitude', 'degrees_West' ) |
---|
1025 | newvar[:] = stdesc[1] |
---|
1026 | |
---|
1027 | newvar = objfile.createVariable( 'latstGDM', 'c', ('nst','StrLength')) |
---|
1028 | Gv = int(stdesc[2]) |
---|
1029 | Dv = int((stdesc[2] - Gv)*60.) |
---|
1030 | Mv = int((stdesc[2] - Gv - Dv/60.)*3600.) |
---|
1031 | writing_str_nc(newvar, [str(Gv)+"d" + str(Dv)+"m" + str(Mv)+'s'], StringLength) |
---|
1032 | |
---|
1033 | if onc.variables.has_key('lat'): |
---|
1034 | print warnmsg |
---|
1035 | print ' ' + fname + ": variable 'lat' already exist !!" |
---|
1036 | print " renaming it as 'latst'" |
---|
1037 | latname = 'latst' |
---|
1038 | else: |
---|
1039 | latname = 'lat' |
---|
1040 | |
---|
1041 | newvar = objfile.createVariable( latname, 'f4', ('nst')) |
---|
1042 | basicvardef(newvar, lonname, 'latitude', 'degrees_North' ) |
---|
1043 | newvar[:] = stdesc[2] |
---|
1044 | |
---|
1045 | if onc.variables.has_key('height'): |
---|
1046 | print warnmsg |
---|
1047 | print ' ' + fname + ": variable 'height' already exist !!" |
---|
1048 | print " renaming it as 'heightst'" |
---|
1049 | heightname = 'heightst' |
---|
1050 | else: |
---|
1051 | heightname = 'height' |
---|
1052 | |
---|
1053 | newvar = objfile.createVariable( heightname, 'f4', ('nst')) |
---|
1054 | basicvardef(newvar, heightname, 'height above sea level', 'm' ) |
---|
1055 | newvar[:] = stdesc[3] |
---|
1056 | |
---|
1057 | return |
---|
1058 | |
---|
1059 | def oper_values(dvals, opers): |
---|
1060 | """ Function to operate the values according to given parameters |
---|
1061 | dvals= datavalues |
---|
1062 | opers= operations |
---|
1063 | """ |
---|
1064 | fname = 'oper_values' |
---|
1065 | |
---|
1066 | newdvals = {} |
---|
1067 | varnames = dvals.keys() |
---|
1068 | |
---|
1069 | aopers = ['sumc','subc','mulc','divc'] |
---|
1070 | |
---|
1071 | Nopers = len(opers) |
---|
1072 | for iop in range(Nopers): |
---|
1073 | vn = varnames[iop] |
---|
1074 | print vn,'oper:',opers[iop] |
---|
1075 | if opers[iop] != '-': |
---|
1076 | opern = opers[iop].split(',')[0] |
---|
1077 | operv = np.float(opers[iop].split(',')[1]) |
---|
1078 | |
---|
1079 | vvals = np.array(dvals[vn]) |
---|
1080 | |
---|
1081 | if opern == 'sumc': |
---|
1082 | newvals = np.where(vvals is None, None, vvals+operv) |
---|
1083 | elif opern == 'subc': |
---|
1084 | newvals = np.where(vvals is None, None, vvals-operv) |
---|
1085 | elif opern == 'mulc': |
---|
1086 | newvals = np.where(vvals is None, None, vvals*operv) |
---|
1087 | elif opern == 'divc': |
---|
1088 | newvals = np.where(vvals is None, None, vvals/operv) |
---|
1089 | else: |
---|
1090 | print errormsg |
---|
1091 | print ' ' + fname + ": operation '" + opern + "' not ready!!" |
---|
1092 | print ' availables:',aopers |
---|
1093 | quit(-1) |
---|
1094 | |
---|
1095 | newdvals[vn] = list(newvals) |
---|
1096 | else: |
---|
1097 | newdvals[vn] = dvals[vn] |
---|
1098 | |
---|
1099 | return newdvals |
---|
1100 | |
---|
1101 | def WMOcodevar(unitsn, onc, Lstrcode=1024): |
---|
1102 | """ Function to add a variabe providing description of a units based on WMO code |
---|
1103 | unitsn= name of the units (all WMO codes derived units must be labelled as |
---|
1104 | 'wmo_code_[num]' associated to a file valled wmo_[num].code) |
---|
1105 | onc= netCDF object file to add the description |
---|
1106 | Lstrcode= length of description of codes |
---|
1107 | |
---|
1108 | wmo_[num].code must have the structure: ('#' for comments) |
---|
1109 | reference|web page, document reference with the code |
---|
1110 | short_description|main description of the code |
---|
1111 | long_description|long description of the code |
---|
1112 | codeTYPE|type of the code (andy of 'D', 'F', 'I', 'S') |
---|
1113 | @| Values line giving the start of the values |
---|
1114 | [val1]|[meaning of first value] |
---|
1115 | (...) |
---|
1116 | [valN]|[meaning of last value] |
---|
1117 | """ |
---|
1118 | fname = 'WMOcodevar' |
---|
1119 | |
---|
1120 | # From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python |
---|
1121 | folder = os.path.dirname(os.path.realpath(__file__)) |
---|
1122 | |
---|
1123 | code = unitsn.split('_')[2] |
---|
1124 | |
---|
1125 | infile = folder + '/wmo_' + code +'.code' |
---|
1126 | |
---|
1127 | if not os.path.isfile(infile): |
---|
1128 | print warnmsg |
---|
1129 | print ' ' + fname + ": WMO code file '" + infile + "' does not exist !!" |
---|
1130 | return |
---|
1131 | # main expected values |
---|
1132 | descvals = ['wmo_code', 'reference', 'short_description', 'long_description', \ |
---|
1133 | 'codeTYPE', '@'] |
---|
1134 | |
---|
1135 | availcodetype = ['D', 'F', 'I', 'S'] |
---|
1136 | |
---|
1137 | codevalues = {} |
---|
1138 | ocode = open(infile, 'r') |
---|
1139 | inivals = False |
---|
1140 | codvals = [] |
---|
1141 | codmeanings = [] |
---|
1142 | for line in ocode: |
---|
1143 | if len(line) > 1 and line[0:1] != '#': |
---|
1144 | linev = line.replace('\n','').replace('\t',' ').replace('\r','') |
---|
1145 | if not inivals: |
---|
1146 | Sv = linev.split('|')[0] |
---|
1147 | Vv = linev.split('|')[1] |
---|
1148 | if searchInlist(descvals, Sv): |
---|
1149 | if Sv != '@': codevalues[Sv] = Vv |
---|
1150 | else: inivals = True |
---|
1151 | else: |
---|
1152 | Svv = linev.split('|')[0] |
---|
1153 | Vvv = linev.split('|')[1] |
---|
1154 | codvals.append(Svv) |
---|
1155 | codmeanings.append(Vvv) |
---|
1156 | |
---|
1157 | # Creating variable |
---|
1158 | if not searchInlist(onc.dimensions, 'Lstringcode'): |
---|
1159 | print ' ' + fname + ": Adding string length dimension 'Lstringcode' for " + \ |
---|
1160 | " code descriptions" |
---|
1161 | newdim = onc.createDimension('Lstringcode', Lstrcode) |
---|
1162 | Ncodes = len(codvals) |
---|
1163 | codedimn = 'wmo_code_' + str(code) |
---|
1164 | if not searchInlist(onc.dimensions, codedimn): |
---|
1165 | print ' ' + fname + ": Adding '" + codedimn + "' dimension for code " + \ |
---|
1166 | "description" |
---|
1167 | newdim = onc.createDimension(codedimn, Ncodes) |
---|
1168 | onc.sync() |
---|
1169 | |
---|
1170 | # Variable with the value of the code |
---|
1171 | if not onc.variables.has_key(codedimn): |
---|
1172 | if codevalues['codeTYPE'] == 'D': |
---|
1173 | newvar = onc.createVariable(codedimn, 'f8', (codedimn)) |
---|
1174 | for iv in range(Ncodes): newvar[iv] = np.float64(codvals[iv]) |
---|
1175 | elif codevalues['codeTYPE'] == 'F': |
---|
1176 | newvar = onc.createVariable(codedimn, 'f', (codedimn)) |
---|
1177 | for iv in range(Ncodes): newvar[iv] = np.float(codvals[iv]) |
---|
1178 | elif codevalues['codeTYPE'] == 'I': |
---|
1179 | newvar = onc.createVariable(codedimn, 'i', (codedimn)) |
---|
1180 | for iv in range(Ncodes): newvar[iv] = int(codvals[iv]) |
---|
1181 | elif codevalues['codeTYPE'] == 'S': |
---|
1182 | newvar = onc.createVariable(codedimn, 'c', (codedimn, 'Lstringcode')) |
---|
1183 | writing_str_nc(newvar, codvals, Lstrcode) |
---|
1184 | else: |
---|
1185 | print errormsg |
---|
1186 | print ' ' + fname + ": codeTYPE '" + codevalues['codeTYPE'] + "' not" + \ |
---|
1187 | " ready !!" |
---|
1188 | print ' available ones:', availcodetype |
---|
1189 | quit(-1) |
---|
1190 | |
---|
1191 | for descv in descvals: |
---|
1192 | if descv != '@' and descv != 'codeTYPE': |
---|
1193 | newvar.setncattr(descv, codevalues[descv]) |
---|
1194 | # Variable with the meaning of the code |
---|
1195 | if not onc.variables.has_key(codedimn+'_meaning'): |
---|
1196 | print ' '+fname+": Adding '" + codedimn + "_meaning' variable for code " + \ |
---|
1197 | "description" |
---|
1198 | newvar = onc.createVariable(codedimn+'_meaning','c',(codedimn,'Lstringcode')) |
---|
1199 | writing_str_nc(newvar, codmeanings, Lstrcode) |
---|
1200 | newvar.setncattr('description', 'meaning of WMO code ' + str(code)) |
---|
1201 | |
---|
1202 | onc.sync() |
---|
1203 | |
---|
1204 | return |
---|
1205 | |
---|
1206 | def EXTRAcodevar(unitsn, onc, Lstrcode=1024): |
---|
1207 | """ Function to add a variabe providing description of a units based on an extra |
---|
1208 | code |
---|
1209 | unitsn= name of the units (all codes derived units must be labelled as |
---|
1210 | 'extra_code_[ref]' associated to a file valled extra_[ref].code) |
---|
1211 | onc= netCDF object file to add the description |
---|
1212 | Lstrcode= length of description of codes |
---|
1213 | |
---|
1214 | extra_[ref].code must have the structure: ('#' for comments) |
---|
1215 | reference|web page, document reference with the code |
---|
1216 | short_description|main description of the code |
---|
1217 | long_description|long description of the code |
---|
1218 | codeTYPE|type of the code (andy of 'D', 'F', 'I', 'S') |
---|
1219 | @| Values line giving the start of the values |
---|
1220 | [val1]|[meaning of first value] |
---|
1221 | (...) |
---|
1222 | [valN]|[meaning of last value] |
---|
1223 | """ |
---|
1224 | fname = 'EXTRAcodevar' |
---|
1225 | |
---|
1226 | # From http://stackoverflow.com/questions/4934806/how-can-i-find-scripts-directory-with-python |
---|
1227 | folder = os.path.dirname(os.path.realpath(__file__)) |
---|
1228 | |
---|
1229 | code = unitsn.split('_')[2] |
---|
1230 | |
---|
1231 | infile = folder + '/extra_' + code +'.code' |
---|
1232 | |
---|
1233 | if not os.path.isfile(infile): |
---|
1234 | print warnmsg |
---|
1235 | print ' ' + fname + ": EXTRA code file '" + infile + "' does not exist !!" |
---|
1236 | return |
---|
1237 | # main expected values |
---|
1238 | descvals = ['reference', 'short_description', 'long_description', \ |
---|
1239 | 'codeTYPE', '@'] |
---|
1240 | |
---|
1241 | availcodetype = ['D', 'F', 'I', 'S'] |
---|
1242 | |
---|
1243 | codevalues = {} |
---|
1244 | ocode = open(infile, 'r') |
---|
1245 | inivals = False |
---|
1246 | codvals = [] |
---|
1247 | codmeanings = [] |
---|
1248 | for line in ocode: |
---|
1249 | if len(line) > 1 and line[0:1] != '#': |
---|
1250 | linev = line.replace('\n','').replace('\t',' ').replace('\r','') |
---|
1251 | if not inivals: |
---|
1252 | Sv = linev.split('|')[0] |
---|
1253 | Vv = linev.split('|')[1] |
---|
1254 | if searchInlist(descvals, Sv): |
---|
1255 | if Sv != '@': codevalues[Sv] = Vv |
---|
1256 | else: inivals = True |
---|
1257 | else: |
---|
1258 | Svv = linev.split('|')[0] |
---|
1259 | Vvv = linev.split('|')[1] |
---|
1260 | codvals.append(Svv) |
---|
1261 | codmeanings.append(Vvv) |
---|
1262 | |
---|
1263 | # Creating variable |
---|
1264 | if not searchInlist(onc.dimensions, 'Lstringcode'): |
---|
1265 | print ' ' + fname + ": Adding string length dimension 'Lstringcode' for " + \ |
---|
1266 | " code descriptions" |
---|
1267 | newdim = onc.createDimension('Lstringcode', Lstrcode) |
---|
1268 | Ncodes = len(codvals) |
---|
1269 | codedimn = 'extra_code_' + str(code) |
---|
1270 | if not searchInlist(onc.dimensions, codedimn): |
---|
1271 | print ' ' + fname + ": Adding '" + codedimn + "' dimension for code " + \ |
---|
1272 | "description" |
---|
1273 | newdim = onc.createDimension(codedimn, Ncodes) |
---|
1274 | onc.sync() |
---|
1275 | |
---|
1276 | # Variable with the value of the code |
---|
1277 | if not onc.variables.has_key(codedimn): |
---|
1278 | if codevalues['codeTYPE'] == 'D': |
---|
1279 | newvar = onc.createVariable(codedimn, 'f8', (codedimn)) |
---|
1280 | for iv in range(Ncodes): newvar[iv] = np.float64(codvals[iv]) |
---|
1281 | elif codevalues['codeTYPE'] == 'F': |
---|
1282 | newvar = onc.createVariable(codedimn, 'f', (codedimn)) |
---|
1283 | for iv in range(Ncodes): newvar[iv] = np.float(codvals[iv]) |
---|
1284 | elif codevalues['codeTYPE'] == 'I': |
---|
1285 | newvar = onc.createVariable(codedimn, 'i', (codedimn)) |
---|
1286 | for iv in range(Ncodes): newvar[iv] = int(codvals[iv]) |
---|
1287 | elif codevalues['codeTYPE'] == 'S': |
---|
1288 | newvar = onc.createVariable(codedimn, 'c', (codedimn, 'Lstringcode')) |
---|
1289 | writing_str_nc(newvar, codvals, Lstrcode) |
---|
1290 | else: |
---|
1291 | print errormsg |
---|
1292 | print ' ' + fname + ": codeTYPE '" + codevalues['codeTYPE'] + "' not" + \ |
---|
1293 | " ready !!" |
---|
1294 | print ' available ones:', availcodetype |
---|
1295 | quit(-1) |
---|
1296 | |
---|
1297 | for descv in descvals: |
---|
1298 | if descv != '@' and descv != 'codeTYPE': |
---|
1299 | newvar.setncattr(descv, codevalues[descv]) |
---|
1300 | # Variable with the meaning of the code |
---|
1301 | if not onc.variables.has_key(codedimn+'_meaning'): |
---|
1302 | print ' '+fname+": Adding '" + codedimn + "_meaning' variable for code " + \ |
---|
1303 | "description" |
---|
1304 | newvar = onc.createVariable(codedimn+'_meaning','c',(codedimn,'Lstringcode')) |
---|
1305 | writing_str_nc(newvar, codmeanings, Lstrcode) |
---|
1306 | newvar.setncattr('description', 'meaning of EXTRA code ' + str(code)) |
---|
1307 | |
---|
1308 | onc.sync() |
---|
1309 | |
---|
1310 | return |
---|
1311 | |
---|
1312 | def add_global_PyNCplot(ObjFile, pyscript, funcname, version): |
---|
1313 | """ Function to add global attributes from 'PyNCplot' to a given netCDF |
---|
1314 | ObjFile= netCDF file object to which add the global attributes |
---|
1315 | funcname= name of the function from which file comes from |
---|
1316 | version= version of the function |
---|
1317 | """ |
---|
1318 | fname = 'add_global_PyNCplot' |
---|
1319 | |
---|
1320 | # Global values |
---|
1321 | ObjFile.setncattr('author', 'L. Fita') |
---|
1322 | newattr = set_attributek(ObjFile, 'institution', unicode('Centro de ' + \ |
---|
1323 | 'Investigaciones del Mar y la Atm' + unichr(243) + 'sfera (CIMA)'), 'U') |
---|
1324 | newattr = set_attributek(ObjFile, 'institution2', unicode('Instituto Franco-' + \ |
---|
1325 | 'Argentino sobre Estudios de Clima y sus Impactos (CNRS, UMI-3351-IFAECI'), 'U') |
---|
1326 | newattr = set_attributek(ObjFile, 'center', unicode('Consejo Nacional de ' + \ |
---|
1327 | 'Investigaciones Cient' + unichr(237) + 'ficas y T' + unichr(233) + \ |
---|
1328 | 'cnicas (CONICET)'), 'U') |
---|
1329 | ObjFile.setncattr('university', 'Universidad de Buenos Aires (UBA)') |
---|
1330 | ObjFile.setncattr('city', 'Buenos Aires') |
---|
1331 | ObjFile.setncattr('country', 'Argentina') |
---|
1332 | ObjFile.setncattr('tool', 'PyNCplot') |
---|
1333 | ObjFile.setncattr('url', 'http://www.xn--llusfb-5va.cat/python/PyNCplot') |
---|
1334 | ObjFile.setncattr('script', pyscript) |
---|
1335 | if funcname is not None: |
---|
1336 | ObjFile.setncattr('function', funcname) |
---|
1337 | ObjFile.setncattr('version', version) |
---|
1338 | |
---|
1339 | ObjFile.sync() |
---|
1340 | |
---|
1341 | return |
---|
1342 | |
---|
1343 | ####### ###### ##### #### ### ## # |
---|
1344 | |
---|
1345 | strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \ |
---|
1346 | " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')" |
---|
1347 | |
---|
1348 | kindobs=['stations-map','multi-points', 'single-station', 'trajectory'] |
---|
1349 | strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \ |
---|
1350 | "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\ |
---|
1351 | "'trajectory': following a trajectory" |
---|
1352 | |
---|
1353 | parser = OptionParser() |
---|
1354 | parser.add_option("-c", "--comments", dest="charcom", |
---|
1355 | help="':', list of characters used for comments", metavar="VALUES") |
---|
1356 | parser.add_option("-d", "--descriptionfile", dest="fdesc", |
---|
1357 | help="description file to use" + read_description.__doc__, metavar="FILE") |
---|
1358 | parser.add_option("-e", "--end_column", dest="endcol", |
---|
1359 | help="character to indicate end of the column ('space', for ' ', or 'None,[fixcoljumps]' for fixed size columns with [fixcoljumps]: ',' separated list of positions of column ends within line)", metavar="VALUE") |
---|
1360 | parser.add_option("-f", "--file", dest="obsfile", |
---|
1361 | help="observational file to use", metavar="FILE") |
---|
1362 | parser.add_option("-j", "--jumpBlines", dest="jumpBlines", |
---|
1363 | help="number of lines to jump from the beggining of file", metavar="VALUE") |
---|
1364 | parser.add_option("-g", "--debug", dest="debug", |
---|
1365 | help="whether debug is required ('false', 'true')", metavar="VALUE") |
---|
1366 | parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, |
---|
1367 | help=strkObs, metavar="VALUE") |
---|
1368 | parser.add_option("-s", "--stationLocation", dest="stloc", |
---|
1369 | help="name ('!' for spaces), longitude, latitude and height of the station (only for 'single-station')", |
---|
1370 | metavar="FILE") |
---|
1371 | parser.add_option("-t", "--CFtime", dest="CFtime", help=strCFt, metavar="VALUE") |
---|
1372 | |
---|
1373 | (opts, args) = parser.parse_args() |
---|
1374 | |
---|
1375 | ####### ####### |
---|
1376 | ## MAIN |
---|
1377 | ####### |
---|
1378 | |
---|
1379 | ofile='OBSnetcdf.nc' |
---|
1380 | |
---|
1381 | if opts.charcom is None: |
---|
1382 | print warnmsg |
---|
1383 | print ' ' + main + ': No list of comment characters provided!!' |
---|
1384 | print ' assuming no need!' |
---|
1385 | charcomments = [] |
---|
1386 | else: |
---|
1387 | charcomments = opts.charcom.split(':') |
---|
1388 | |
---|
1389 | if opts.jumpBlines is None: |
---|
1390 | print warnmsg |
---|
1391 | print ' ' + main + ': No number of lines to jump from beggining of file provided!!' |
---|
1392 | print ' assuming no need!' |
---|
1393 | jumpBlines = 0 |
---|
1394 | else: |
---|
1395 | jumpBlines = int(opts.jumpBlines) |
---|
1396 | |
---|
1397 | |
---|
1398 | if opts.fdesc is None: |
---|
1399 | print errormsg |
---|
1400 | print ' ' + main + ': No description file for the observtional data provided!!' |
---|
1401 | quit(-1) |
---|
1402 | |
---|
1403 | if opts.endcol is None: |
---|
1404 | print warnmsg |
---|
1405 | print ' ' + main + ': No list of comment characters provided!!' |
---|
1406 | print " assuming 'space'" |
---|
1407 | endcol = ' ' |
---|
1408 | else: |
---|
1409 | if opts.endcol == 'space': |
---|
1410 | endcol = ' ' |
---|
1411 | else: |
---|
1412 | endcol = opts.endcol |
---|
1413 | |
---|
1414 | if opts.obsfile is None: |
---|
1415 | print errormsg |
---|
1416 | print ' ' + main + ': No observations file provided!!' |
---|
1417 | quit(-1) |
---|
1418 | |
---|
1419 | if opts.debug is None: |
---|
1420 | print warnmsg |
---|
1421 | print ' ' + main + ': No debug provided!!' |
---|
1422 | print " assuming 'False'" |
---|
1423 | debug = False |
---|
1424 | else: |
---|
1425 | if opts.debug == 'true': |
---|
1426 | debug = True |
---|
1427 | else: |
---|
1428 | debug = False |
---|
1429 | |
---|
1430 | if not os.path.isfile(opts.fdesc): |
---|
1431 | print errormsg |
---|
1432 | print ' ' + main + ": description file '" + opts.fdesc + "' does not exist !!" |
---|
1433 | quit(-1) |
---|
1434 | |
---|
1435 | if not os.path.isfile(opts.obsfile): |
---|
1436 | print errormsg |
---|
1437 | print ' ' + main + ": observational file '" + opts.obsfile + "' does not exist !!" |
---|
1438 | quit(-1) |
---|
1439 | |
---|
1440 | if opts.CFtime is None: |
---|
1441 | print warnmsg |
---|
1442 | print ' ' + main + ': No CFtime criteria are provided !!' |
---|
1443 | print " either time is already in CF-format ('timeFMT=CFtime') in '" + \ |
---|
1444 | opts.fdesc + "'" |
---|
1445 | print " or assuming refdate: '19491201000000' and time units: 'hours'" |
---|
1446 | referencedate = '19491201000000' |
---|
1447 | timeunits = 'hours' |
---|
1448 | else: |
---|
1449 | referencedate = opts.CFtime.split(',')[0] |
---|
1450 | timeunits = opts.CFtime.split(',')[1] |
---|
1451 | |
---|
1452 | if opts.obskind is None: |
---|
1453 | print warnmsg |
---|
1454 | print ' ' + main + ': No kind of observations provided !!' |
---|
1455 | print " assuming 'single-station': single station on a fixed position at 0,0,0" |
---|
1456 | obskind = 'single-station' |
---|
1457 | stationdesc = [0.0, 0.0, 0.0, 0.0] |
---|
1458 | else: |
---|
1459 | obskind = opts.obskind |
---|
1460 | if obskind == 'single-station': |
---|
1461 | if opts.stloc is None: |
---|
1462 | print errormsg |
---|
1463 | print ' ' + main + ': No station location provided !!' |
---|
1464 | quit(-1) |
---|
1465 | else: |
---|
1466 | stvals = opts.stloc.split(',') |
---|
1467 | stationdesc = [stvals[0], np.float(stvals[1]), np.float(stvals[2]), \ |
---|
1468 | np.float(stvals[3])] |
---|
1469 | else: |
---|
1470 | obskind = opts.obskind |
---|
1471 | |
---|
1472 | # Reading description file |
---|
1473 | ## |
---|
1474 | description = read_description(opts.fdesc, debug) |
---|
1475 | |
---|
1476 | Nvariables=len(description['varN']) |
---|
1477 | NlongN=len(description['varLN']) |
---|
1478 | NvarU=len(description['varU']) |
---|
1479 | formats = description['varTYPE'] |
---|
1480 | |
---|
1481 | if Nvariables != NlongN: |
---|
1482 | print errormsg |
---|
1483 | print ' ' + main + ': number of variables:', Nvariables,' and number of ' + \ |
---|
1484 | 'long names', NlongN,' does not coincide!!' |
---|
1485 | print ' * what is found is _______' |
---|
1486 | if Nvariables > NlongN: |
---|
1487 | Nshow = NlongN |
---|
1488 | for ivar in range(Nshow): |
---|
1489 | print ivar,' ',description['varN'][ivar],':', description['varLN'][ivar],\ |
---|
1490 | '[', description['varU'][ivar], '] fmt:', formats[ivar] |
---|
1491 | print ' missing values for:', description['varN'][Nshow:Nvariables+1] |
---|
1492 | else: |
---|
1493 | Nshow = Nvariables |
---|
1494 | for ivar in range(Nshow): |
---|
1495 | print ivar,' ',description['varN'][ivar],':', description['varLN'][ivar],\ |
---|
1496 | '[', description['varU'][ivar], '] fmt:', formats[ivar] |
---|
1497 | print ' excess of long names for :', description['varLN'][Nshow:NlongN+1] |
---|
1498 | |
---|
1499 | quit(-1) |
---|
1500 | |
---|
1501 | if Nvariables != NvarU: |
---|
1502 | print errormsg |
---|
1503 | print ' ' + main + ': number of variables:', Nvariables,' and number of ' + \ |
---|
1504 | 'units', NvarU,' does not coincide!!' |
---|
1505 | print ' * what is found is _______' |
---|
1506 | if Nvariables > NvarU: |
---|
1507 | Nshow = NvarU |
---|
1508 | for ivar in range(Nshow): |
---|
1509 | print ivar,' ',description['varN'][ivar],':', description['varLN'][ivar],\ |
---|
1510 | '[', description['varU'][ivar], '] fmt:', formats[ivar] |
---|
1511 | print ' missing values for:', description['varN'][Nshow:Nvariables+1] |
---|
1512 | else: |
---|
1513 | Nshow = Nvariables |
---|
1514 | for ivar in range(Nshow): |
---|
1515 | print ivar,' ',description['varN'][ivar],':', description['varLN'][ivar],\ |
---|
1516 | '[', description['varU'][ivar], '] fmt:', formats[ivar] |
---|
1517 | print ' excess of units for :', description['varU'][Nshow:NvarU+1] |
---|
1518 | |
---|
1519 | quit(-1) |
---|
1520 | |
---|
1521 | print 'Number of variables', Nvariables |
---|
1522 | print 'Number of long names', len(description['varLN']) |
---|
1523 | |
---|
1524 | if len(formats) != Nvariables: |
---|
1525 | print errormsg |
---|
1526 | print ' ' + main + ': number of formats:',len(formats),' and number of ' + \ |
---|
1527 | 'variables', Nvariables,' does not coincide!!' |
---|
1528 | print ' * what is found is _______' |
---|
1529 | if Nvariables > len(formats): |
---|
1530 | Nshow = len(formats) |
---|
1531 | for ivar in range(Nshow): |
---|
1532 | print ivar,' ',description['varN'][ivar],':', description['varLN'][ivar],\ |
---|
1533 | '[', description['varU'][ivar], '] fmt:', formats[ivar] |
---|
1534 | print ' missing values for:', description['varN'][Nshow:Nvariables+1] |
---|
1535 | else: |
---|
1536 | Nshow = Nvariables |
---|
1537 | for ivar in range(Nshow): |
---|
1538 | print ivar,' ',description['varN'][ivar],':', description['varLN'][ivar],\ |
---|
1539 | '[', description['varU'][ivar], '] fmt:', formats[ivar] |
---|
1540 | print ' excess of formats for:', formats[Nshow:len(formats)+1] |
---|
1541 | |
---|
1542 | # quit(-1) |
---|
1543 | |
---|
1544 | # Reading values |
---|
1545 | ## |
---|
1546 | datavalues = read_datavalues(opts.obsfile, charcomments, endcol, formats, jumpBlines, |
---|
1547 | description['varOPER'], description['MissingValue'], description['varN'], debug) |
---|
1548 | |
---|
1549 | # Total number of values |
---|
1550 | Ntvalues = len(datavalues[description['varN'][0]]) |
---|
1551 | if obskind == 'stations-map': |
---|
1552 | print main + ': total values found:',Ntvalues |
---|
1553 | else: |
---|
1554 | print main + ': total temporal values found:',Ntvalues |
---|
1555 | |
---|
1556 | objfile = NetCDFFile(ofile, 'w') |
---|
1557 | |
---|
1558 | # Creation of dimensions |
---|
1559 | ## |
---|
1560 | if obskind == 'stations-map': |
---|
1561 | rowsdim = 'Npoints' |
---|
1562 | dimlength = Ntvalues |
---|
1563 | else: |
---|
1564 | rowsdim = 'time' |
---|
1565 | dimlength = None |
---|
1566 | |
---|
1567 | objfile.createDimension(rowsdim,dimlength) |
---|
1568 | objfile.createDimension('StrLength',StringLength) |
---|
1569 | |
---|
1570 | # Creation of variables |
---|
1571 | ## |
---|
1572 | for ivar in range(Nvariables): |
---|
1573 | varn = description['varN'][ivar] |
---|
1574 | print " including: '" + varn + "' ... .. ." |
---|
1575 | |
---|
1576 | if formats[ivar] == 'D': |
---|
1577 | newvar = objfile.createVariable(varn, 'f32', (rowsdim), fill_value=fillValueF) |
---|
1578 | basicvardef(newvar, varn, description['varLN'][ivar], \ |
---|
1579 | description['varU'][ivar]) |
---|
1580 | newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn]) |
---|
1581 | elif formats[ivar] == 'F': |
---|
1582 | newvar = objfile.createVariable(varn, 'f', (rowsdim), fill_value=fillValueF) |
---|
1583 | basicvardef(newvar, varn, description['varLN'][ivar], \ |
---|
1584 | description['varU'][ivar]) |
---|
1585 | newvar[:] = np.where(datavalues[varn] is None, fillValueF, datavalues[varn]) |
---|
1586 | elif formats[ivar] == 'I': |
---|
1587 | newvar = objfile.createVariable(varn, 'i', (rowsdim), fill_value=fillValueI) |
---|
1588 | basicvardef(newvar, varn, description['varLN'][ivar], \ |
---|
1589 | description['varU'][ivar]) |
---|
1590 | # Why is not wotking with integers? |
---|
1591 | vals = np.array(datavalues[varn]) |
---|
1592 | for iv in range(Ntvalues): |
---|
1593 | if vals[iv] is None: vals[iv] = fillValueI |
---|
1594 | newvar[:] = vals |
---|
1595 | elif formats[ivar] == 'S': |
---|
1596 | newvar = objfile.createVariable(varn, 'c', (rowsdim,'StrLength')) |
---|
1597 | basicvardef(newvar, varn, description['varLN'][ivar], \ |
---|
1598 | description['varU'][ivar]) |
---|
1599 | vals = datavalues[varn] |
---|
1600 | for iv in range(Ntvalues): |
---|
1601 | if vals[iv] is None: vals[iv] = fillValueS |
---|
1602 | writing_str_nc(newvar, vals, StringLength) |
---|
1603 | |
---|
1604 | # Getting new variables to describe certain units as codeWMO_[num] from an |
---|
1605 | # external file |
---|
1606 | if description['varU'][ivar][0:8] == 'wmo_code': |
---|
1607 | WMOcodevar(description['varU'][ivar], objfile) |
---|
1608 | # Getting new variables to describe certain units as codeEXTRA_[ref] from an |
---|
1609 | # external file |
---|
1610 | if description['varU'][ivar][0:10] == 'extra_code': |
---|
1611 | EXTRAcodevar(description['varU'][ivar], objfile) |
---|
1612 | |
---|
1613 | # Extra variable descriptions/attributes |
---|
1614 | if description.has_key('varBUFR'): |
---|
1615 | set_attribute(newvar,'bufr_code',description['varBUFR'][ivar]) |
---|
1616 | |
---|
1617 | objfile.sync() |
---|
1618 | |
---|
1619 | # Time variable in CF format |
---|
1620 | ## |
---|
1621 | if description['FMTtime'] == 'CFtime': |
---|
1622 | timevals = datavalues[description['NAMEtime']] |
---|
1623 | iv = 0 |
---|
1624 | for ivn in description['varN']: |
---|
1625 | if ivn == description['NAMEtime']: |
---|
1626 | tunits = description['varU'][iv] |
---|
1627 | break |
---|
1628 | iv = iv + 1 |
---|
1629 | else: |
---|
1630 | # Time as a composition of different columns |
---|
1631 | tcomposite = description['NAMEtime'].find('@') |
---|
1632 | if tcomposite != -1: |
---|
1633 | timevars = description['NAMEtime'].split('@') |
---|
1634 | |
---|
1635 | print warnmsg |
---|
1636 | print ' ' + main + ': time values as combination of different columns!' |
---|
1637 | print ' combining:',timevars,' with a final format: ',description['FMTtime'] |
---|
1638 | |
---|
1639 | timeSvals = [] |
---|
1640 | if debug: print ' ' + main + ': creating times _______' |
---|
1641 | for it in range(Ntvalues): |
---|
1642 | tSvals = '' |
---|
1643 | for tvar in timevars: |
---|
1644 | tSvals = tSvals + datavalues[tvar][it] + ' ' |
---|
1645 | |
---|
1646 | timeSvals.append(tSvals[0:len(tSvals)-1]) |
---|
1647 | if debug: print it, '*' + timeSvals[it] + '*' |
---|
1648 | |
---|
1649 | timevals = Stringtimes_CF(timeSvals, description['FMTtime'].replace('@',' '),\ |
---|
1650 | referencedate, timeunits, debug) |
---|
1651 | else: |
---|
1652 | timevals = Stringtimes_CF(datavalues[description['NAMEtime']], \ |
---|
1653 | description['FMTtime'], referencedate, timeunits, debug) |
---|
1654 | |
---|
1655 | CFtimeRef = referencedate[0:4] +'-'+ referencedate[4:6] +'-'+ referencedate[6:8] + \ |
---|
1656 | ' ' + referencedate[8:10] +':'+ referencedate[10:12] +':'+ referencedate[12:14] |
---|
1657 | tunits = timeunits + ' since ' + CFtimeRef |
---|
1658 | |
---|
1659 | if objfile.variables.has_key('time'): |
---|
1660 | print warnmsg |
---|
1661 | print ' ' + main + ": variable 'time' already exist !!" |
---|
1662 | print " renaming it as 'CFtime'" |
---|
1663 | timeCFname = 'CFtime' |
---|
1664 | newdim = objfile.renameDimension('time','CFtime') |
---|
1665 | newvar = objfile.createVariable( timeCFname, 'f8', ('CFtime')) |
---|
1666 | basicvardef(newvar, timeCFname, 'time', tunits ) |
---|
1667 | else: |
---|
1668 | if not searchInlist(objfile.dimensions, 'time'): |
---|
1669 | newdim = objfile.createDimension('time',None) |
---|
1670 | timeCFname = 'time' |
---|
1671 | newvar = objfile.createVariable( timeCFname, 'f8', ('time')) |
---|
1672 | newvar[:] = np.zeros(timevals.shape[0]) |
---|
1673 | basicvardef(newvar, timeCFname, 'time', tunits ) |
---|
1674 | |
---|
1675 | set_attribute(newvar, 'calendar', 'standard') |
---|
1676 | if obskind == 'stations-map': |
---|
1677 | newvar[:] = timevals[0] |
---|
1678 | else: |
---|
1679 | newvar[:] = timevals |
---|
1680 | |
---|
1681 | # Global attributes |
---|
1682 | ## |
---|
1683 | for descn in description.keys(): |
---|
1684 | if descn[0:3] != 'var' and descn[0:4] != 'NAME' and descn[0:3] != 'FMT': |
---|
1685 | set_attribute(objfile, descn, description[descn]) |
---|
1686 | |
---|
1687 | add_global_PyNCplot(objfile, main, 'main', version) |
---|
1688 | |
---|
1689 | objfile.sync() |
---|
1690 | |
---|
1691 | # Adding new variables as function of the observational type |
---|
1692 | ## 'multi-points', 'single-station', 'trajectory' |
---|
1693 | |
---|
1694 | if obskind != 'single-station': |
---|
1695 | adding_complementary(objfile, description, obskind, datavalues, timevals, \ |
---|
1696 | referencedate, timeunits, Ndim2D, debug) |
---|
1697 | else: |
---|
1698 | # Adding three variables with the station location, longitude, latitude and height |
---|
1699 | adding_station_desc(objfile,stationdesc) |
---|
1700 | |
---|
1701 | objfile.sync() |
---|
1702 | objfile.close() |
---|
1703 | |
---|
1704 | print main + ": Successfull generation of netcdf observational file '" + ofile + "' !!" |
---|