1 __version__ = '1.9.3'
2 import g2clib
3 import struct
4 import string
5 import math
6 import warnings
7 import operator
8 from datetime import datetime
9 try:
10 from StringIO import StringIO
11 except ImportError:
12 from io import BytesIO as StringIO
13
14 import numpy as np
15 from numpy import ma
16 try:
17 import pyproj
18 except ImportError:
19 try:
20 from mpl_toolkits.basemap import pyproj
21 except:
22 raise ImportError("either pyproj or basemap required")
23
24
25 _earthparams={0:6367470.0,
26 1:'Spherical - radius specified in m by data producer',
27 2:(6378160.0,6356775.0),
28 3:'OblateSpheroid - major and minor axes specified in km by data producer',
29 4:(6378137.0,6356752.314),
30 5:'WGS84',
31 6:6371229.0,
32 7:'OblateSpheroid - major and minor axes specified in m by data producer',
33 8:6371200.0,
34 255:'Missing'}
35 for _n in range(192):
36 if not _n in _earthparams: _earthparams[_n]='Reserved'
37 for _n in range(192,255):
38 _earthparams[_n] = 'Reserved for local use'
39
40 _table0={1:('Melbourne (WMC)','ammc'),
41 2:('Melbourne - BMRC (WMC)',None),
42 3:('Melbourne (WMC)',None),
43 4:('Moscow (WMC)',None),
44 5:('Moscow (WMC)',None),
45 6:('Moscow (WMC)',None),
46 7:('US National Weather Service - NCEP (WMC)','kwbc'),
47 8:('US National Weather Service - NWSTG (WMC)',None),
48 9:('US National Weather Service - Other (WMC)',None),
49 10:('Cairo (RSMC/RAFC)',None),
50 11:('Cairo (RSMC/RAFC)',None),
51 12:('Dakar (RSMC/RAFC)',None),
52 13:('Dakar (RSMC/RAFC)',None),
53 14:('Nairobi (RSMC/RAFC)',None),
54 15:('Nairobi (RSMC/RAFC)',None),
55 16:('Casablanca',None),
56 17:('Tunis (RSMC)',None),
57 18:('Tunis-Casablanca (RSMC)',None),
58 19:('Tunis-Casablanca (RSMC)',None),
59 20:('Las Palmas (RAFC)',None),
60 21:('Algiers (RSMC)',None),
61 22:('ACMAD',None),
62 23:('Mozambique (NMC)',None),
63 24:('Pretoria (RSMC)',None),
64 25:('La Reunion (RSMC)',None),
65 26:('Khabarovsk (RSMC)',None),
66 27:('Khabarovsk (RSMC)',None),
67 28:('New Delhi (RSMC/RAFC)',None),
68 29:('New Delhi (RSMC/RAFC)',None),
69 30:('Novosibirsk (RSMC)',None),
70 31:('Novosibirsk (RSMC)',None),
71 32:('Tashkent (RSMC)',None),
72 33:('Jeddah (RSMC)',None),
73 34:('Japanese Meteorological Agency - Tokyo (RSMC)','rjtd'),
74 35:('Japanese Meteorological Agency - Tokyo (RSMC)',None),
75 36:('Bankok',None),
76 37:('Ulan Bator',None),
77 38:('Beijing (RSMC)','babj'),
78 39:('Beijing (RSMC)',None),
79 40:('Korean Meteorological Administration - Seoul','rksl'),
80 41:('Buenos Aires (RSMC/RAFC)',None),
81 42:('Buenos Aires (RSMC/RAFC)',None),
82 43:('Brasilia (RSMC/RAFC)',None),
83 44:('Brasilia (RSMC/RAFC)',None),
84 45:('Santiago',None),
85 46:('Brazilian Space Agency - INPE','sbsj'),
86 47:('Columbia (NMC)',None),
87 48:('Ecuador (NMC)',None),
88 49:('Peru (NMC)',None),
89 50:('Venezuela (NMC)',None),
90 51:('Miami (RSMC/RAFC)',None),
91 52:('Tropical Prediction Center (NHC), Miami (RSMC)',None),
92 53:('Canadian Meteorological Service - Montreal (RSMC)',None),
93 54:('Canadian Meteorological Service - Montreal (RSMC)','cwao'),
94 55:('San Francisco',None),
95 56:('ARINC Center',None),
96 57:('U.S. Air Force - Global Weather Center',None),
97 58:('US Navy - Fleet Numerical Oceanography Center','fnmo'),
98 59:('NOAA Forecast Systems Lab, Boulder CO',None),
99 60:('National Center for Atmospheric Research (NCAR), Boulder, CO',None),
100 61:('Service ARGOS - Landover, MD, USA',None),
101 62:('US Naval Oceanographic Office',None),
102 63:('Reserved',None),
103 64:('Honolulu',None),
104 65:('Darwin (RSMC)',None),
105 66:('Darwin (RSMC)',None),
106 67:('Melbourne (RSMC)',None),
107 68:('Reserved',None),
108 69:('Wellington (RSMC/RAFC)',None),
109 70:('Wellington (RSMC/RAFC)',None),
110 71:('Nadi (RSMC)',None),
111 72:('Singapore',None),
112 73:('Malaysia (NMC)',None),
113 74:('U.K. Met Office - Exeter (RSMC)','egrr'),
114 75:('U.K. Met Office - Exeter (RSMC)',None),
115 76:('Moscow (RSMC/RAFC)',None),
116 77:('Reserved',None),
117 78:('Offenbach (RSMC)','edzw'),
118 79:('Offenbach (RSMC)',None),
119 80:('Rome (RSMC)','cnmc'),
120 81:('Rome (RSMC)',None),
121 82:('Norrkoping',None),
122 83:('Norrkoping',None),
123 84:('French Weather Service - Toulouse','lfpw'),
124 85:('French Weather Service - Toulouse','lfpw'),
125 86:('Helsinki',None),
126 87:('Belgrade',None),
127 88:('Oslo',None),
128 89:('Prague',None),
129 90:('Episkopi',None),
130 91:('Ankara',None),
131 92:('Frankfurt/Main (RAFC)',None),
132 93:('London (WAFC)',None),
133 94:('Copenhagen',None),
134 95:('Rota',None),
135 96:('Athens',None),
136 97:('European Space Agency (ESA)',None),
137 98:('European Center for Medium-Range Weather Forecasts (RSMC)','ecmf'),
138 99:('De BiltNone), Netherlands',None),
139 100:('Brazzaville',None),
140 101:('Abidjan',None),
141 102:('Libyan Arab Jamahiriya (NMC)',None),
142 103:('Madagascar (NMC)',None),
143 104:('Mauritius (NMC)',None),
144 105:('Niger (NMC)',None),
145 106:('Seychelles (NMC)',None),
146 107:('Uganda (NMC)',None),
147 108:('Tanzania (NMC)',None),
148 109:('Zimbabwe (NMC)',None),
149 110:('Hong-Kong',None),
150 111:('Afghanistan (NMC)',None),
151 112:('Bahrain (NMC)',None),
152 113:('Bangladesh (NMC)',None),
153 114:('Bhutan (NMC)',None),
154 115:('Cambodia (NMC)',None),
155 116:("Democratic People's Republic of Korea (NMC)",None),
156 117:('Islamic Republic of Iran (NMC)',None),
157 118:('Iraq (NMC)',None),
158 119:('Kazakhstan (NMC)',None),
159 120:('Kuwait (NMC)',None),
160 121:('Kyrgyz Republic (NMC)',None),
161 122:("Lao People's Democratic Republic (NMC)",None),
162 123:('MacaoNone), China',None),
163 124:('Maldives (NMC)',None),
164 125:('Myanmar (NMC)',None),
165 126:('Nepal (NMC)',None),
166 127:('Oman (NMC)',None),
167 128:('Pakistan (NMC)',None),
168 129:('Qatar (NMC)',None),
169 130:('Republic of Yemen (NMC)',None),
170 131:('Sri Lanka (NMC)',None),
171 132:('Tajikistan (NMC)',None),
172 133:('Turkmenistan (NMC)',None),
173 134:('United Arab Emirates (NMC)',None),
174 135:('Uzbekistan (NMC)',None),
175 136:('Socialist Republic of Viet Nam (NMC)',None),
176 137:('Reserved',None),
177 138:('Reserved',None),
178 139:('Reserved',None),
179 140:('Bolivia (NMC)',None),
180 141:('Guyana (NMC)',None),
181 142:('Paraguay (NMC)',None),
182 143:('Suriname (NMC)',None),
183 144:('Uruguay (NMC)',None),
184 145:('French Guyana',None),
185 146:('Brazilian Navy Hydrographic Center',None),
186 147:('Reserved',None),
187 148:('Reserved',None),
188 149:('Reserved',None),
189 150:('Antigua and Barbuda (NMC)',None),
190 151:('Bahamas (NMC)',None),
191 152:('Barbados (NMC)',None),
192 153:('Belize (NMC)',None),
193 154:('British Caribbean Territories Center',None),
194 155:('San Jose',None),
195 156:('Cuba (NMC)',None),
196 157:('Dominica (NMC)',None),
197 158:('Dominican Republic (NMC)',None),
198 159:('El Salvador (NMC)',None),
199 160:('US NOAA/NESDIS',None),
200 161:('US NOAA Office of Oceanic and Atmospheric Research',None),
201 162:('Guatemala (NMC)',None),
202 163:('Haiti (NMC)',None),
203 164:('Honduras (NMC)',None),
204 165:('Jamaica (NMC)',None),
205 166:('Mexico',None),
206 167:('Netherlands Antilles and Aruba (NMC)',None),
207 168:('Nicaragua (NMC)',None),
208 169:('Panama (NMC)',None),
209 170:('Saint Lucia (NMC)',None),
210 171:('Trinidad and Tobago (NMC)',None),
211 172:('French Departments',None),
212 173:('Reserved',None),
213 174:('Reserved',None),
214 175:('Reserved',None),
215 176:('Reserved',None),
216 177:('Reserved',None),
217 178:('Reserved',None),
218 179:('Reserved',None),
219 180:('Reserved',None),
220 181:('Reserved',None),
221 182:('Reserved',None),
222 183:('Reserved',None),
223 184:('Reserved',None),
224 185:('Reserved',None),
225 186:('Reserved',None),
226 187:('Reserved',None),
227 188:('Reserved',None),
228 189:('Reserved',None),
229 190:('Cook Islands (NMC)',None),
230 191:('French Polynesia (NMC)',None),
231 192:('Tonga (NMC)',None),
232 193:('Vanuatu (NMC)',None),
233 194:('Brunei (NMC)',None),
234 195:('Indonesia (NMC)',None),
235 196:('Kiribati (NMC)',None),
236 197:('Federated States of Micronesia (NMC)',None),
237 198:('New Caledonia (NMC)',None),
238 199:('Niue',None),
239 200:('Papua New Guinea (NMC)',None),
240 201:('Philippines (NMC)',None),
241 202:('Samoa (NMC)',None),
242 203:('Solomon Islands (NMC)',None),
243 204:('Reserved',None),
244 205:('Reserved',None),
245 206:('Reserved',None),
246 207:('Reserved',None),
247 208:('Reserved',None),
248 209:('Reserved',None),
249 210:('Frascati',None),
250 211:('Lanion',None),
251 212:('Lisboa',None),
252 213:('Reykjavik',None),
253 214:('Madrid','lemm'),
254 215:('Zurich',None),
255 216:('Service ARGOS - ToulouseNone), FR',None),
256 217:('Bratislava',None),
257 218:('Budapest',None),
258 219:('Ljubljana',None),
259 220:('Warsaw',None),
260 221:('Zagreb',None),
261 222:('Albania (NMC)',None),
262 223:('Armenia (NMC)',None),
263 224:('Austria (NMC)',None),
264 225:('Azerbaijan (NMC)',None),
265 226:('Belarus (NMC)',None),
266 227:('Belgium (NMC)',None),
267 228:('Bosnia and Herzegovina (NMC)',None),
268 229:('Bulgaria (NMC)',None),
269 230:('Cyprus (NMC)',None),
270 231:('Estonia (NMC)',None),
271 232:('Georgia (NMC)',None),
272 233:('Dublin',None),
273 234:('Israel (NMC)',None),
274 235:('Jordan (NMC)',None),
275 236:('Latvia (NMC)',None),
276 237:('Lebanon (NMC)',None),
277 238:('Lithuania (NMC)',None),
278 239:('Luxembourg',None),
279 240:('Malta (NMC)',None),
280 241:('Monaco',None),
281 242:('Romania (NMC)',None),
282 243:('Syrian Arab Republic (NMC)',None),
283 244:('The former Yugoslav Republic of Macedonia (NMC)',None),
284 245:('Ukraine (NMC)',None),
285 246:('Republic of Moldova',None),
286 247:('Reserved',None),
287 248:('Reserved',None),
288 249:('Reserved',None),
289 250:('Reserved',None),
290 251:('Reserved',None),
291 252:('Reserved',None),
292 253:('Reserved',None),
293 254:('EUMETSAT Operations Center',None),
294 255:('Missing Value',None)}
295
297 """
298 A decimal to binary converter. Returns bits in a list.
299 """
300 retval = []
301 for i in range(maxbits - 1, -1, -1):
302 bit = int(val / (2 ** i))
303 val = (val % (2 ** i))
304 retval.append(bit)
305 return retval
306
308 """convert a float to a IEEE format 32 bit integer"""
309 ra = np.array([r],'f')
310 ia = np.empty(1,'i')
311 g2clib.rtoi_ieee(ra,ia)
312 return ia[0]
313
315 """convert an IEEE format 32 bit integer to a float"""
316 ia = np.array([i],'i')
317 ra = np.empty(1,'f')
318 g2clib.itor_ieee(ia,ra)
319 return ra[0]
320
322 """Test if string is a string like object if not return 0 """
323 try: string + ''
324 except: return 0
325 else: return 1
326
328 """
329 Class for accessing data in a GRIB Edition 2 message.
330
331 The L{Grib2Decode} function returns a list of these class instances,
332 one for each grib message in the file.
333
334 When a class instance is created, metadata in the GRIB2 file
335 is decoded and used to set various instance variables.
336
337 @ivar bitmap_indicator_flag: flag to indicate whether a bit-map is used (0 for yes, 255 for no).
338 @ivar data_representation_template: data representation template from section 5.
339 @ivar data_representation_template_number: data representation template number
340 from section 5
341 (U{Table 5.0
342 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table5-0.shtml>})
343 @ivar has_local_use_section: True if grib message contains a local use
344 section. If True the actual local use section is contained in the
345 C{_local_use_section} instance variable, as a raw byte string.
346 @ivar discipline_code: product discipline code for grib message
347 (U{Table 0.0
348 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table0-0.shtml>}).
349 @ivar earthRmajor: major (equatorial) earth radius.
350 @ivar earthRminor: minor (polar) earth radius.
351 @ivar grid_definition_info: grid definition section information from section 3.
352 See L{Grib2Encode.addgrid} for details.
353 @ivar grid_definition_template: grid definition template from section 3.
354 @ivar grid_definition_template_number: grid definition template number from section 3
355 (U{Table 3.1
356 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-1.shtml>}).
357 @ivar gridlength_in_x_direction: x (or longitudinal) direction grid length.
358 @ivar gridlength_in_y_direction: y (or latitudinal) direction grid length.
359 @ivar identification_section: data from identification section (section 1).
360 See L{Grib2Encode.__init__} for details.
361 @ivar latitude_first_gridpoint: latitude of first grid point on grid.
362 @ivar latitude_last_gridpoint: latitude of last grid point on grid.
363 @ivar longitude_first_gridpoint: longitude of first grid point on grid.
364 @ivar longitude_last_gridpoint: longitude of last grid point on grid.
365 @ivar originating_center: name of national/international originating center.
366 @ivar center_wmo_code: 4 character wmo code for originating center.
367 @ivar scanmodeflags: scanning mode flags from Table 3.4
368 (U{Table 3.4
369 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-4.shtml>}).
370
371 - bit 1:
372
373 0 - Points in the first row or column scan in the +i (+x) direction
374
375 1 - Points in the first row or column scan in the -i (-x) direction
376
377 - bit 2:
378
379 0 - Points in the first row or column scan in the -j (-y) direction
380
381 1 - Points in the first row or column scan in the +j (+y) direction
382
383 - bit 3:
384
385 0 - Adjacent points in the i (x) direction are consecutive (row-major order).
386
387 1 - Adjacent points in the j (y) direction are consecutive (column-major order).
388
389 - bit 4:
390
391 0 - All rows scan in the same direction
392
393 1 - Adjacent rows scan in the opposite direction
394
395 @ivar number_of_data_points_to_unpack: total number of data points in grib message.
396 @ivar points_in_x_direction: number of points in the x (longitudinal) direction.
397 @ivar points_in_y_direction: number of points in the y (latitudinal) direction.
398 @ivar product_definition_template: product definition template from section 4.
399 @ivar product_definition_template_number: product definition template number from section 4
400 (U{Table 4.0
401 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-0.shtml>}).
402 @ivar shape_of_earth: string describing the shape of the earth (e.g. 'Oblate Spheroid', 'Spheroid').
403 @ivar spectral_truncation_parameters: pentagonal truncation parameters that describe the
404 spherical harmonic truncation (only relevant for grid_definition_template_numbers 50-52).
405 For triangular truncation, all three of these numbers are the same.
406 @ivar latitude_of_southern_pole: the geographic latitude in degrees of the southern
407 pole of the coordinate system (for rotated lat/lon or gaussian grids).
408 @ivar longitude_of_southern_pole: the geographic longitude in degrees of the southern
409 pole of the coordinate system (for rotated lat/lon or gaussian grids).
410 @ivar angle_of_pole_rotation: The angle of rotation in degrees about the new
411 polar axis (measured clockwise when looking from the southern to the northern pole)
412 of the coordinate system. For rotated lat/lon or gaussian grids.
413 @ivar missing_value: primary missing value (for data_representation_template_numbers
414 2 and 3).
415 @ivar missing_value2: secondary missing value (for data_representation_template_numbers
416 2 and 3).
417 @ivar proj4_: instance variables with this prefix are used to set the map projection
418 parameters for U{PROJ.4<http://proj.maptools.org>}.
419 """
421 """
422 create a Grib2Decode class instance given a GRIB Edition 2 filename.
423
424 (used by L{Grib2Decode} function - not directly called by user)
425 """
426 for k,v in kwargs.items():
427 setattr(self,k,v)
428
429 gdsinfo = self.grid_definition_info
430 gdtnum = self.grid_definition_template_number
431 gdtmpl = self.grid_definition_template
432 reggrid = gdsinfo[2] == 0
433
434 if gdtnum not in [50,51,52,1200]:
435 earthR = _earthparams[gdtmpl[0]]
436 if earthR == 'Reserved': earthR = None
437 else:
438 earthR = None
439 if _isString(earthR) and (earthR.startswith('Reserved') or earthR=='Missing'):
440 self.shape_of_earth = earthR
441 self.earthRminor = None
442 self.earthRmajor = None
443 elif _isString(earthR) and earthR.startswith('Spherical'):
444 self.shape_of_earth = earthR
445 scaledearthR = gdtmpl[2]
446 earthRscale = gdtmpl[1]
447 self.earthRmajor = math.pow(10,-earthRscale)*scaledearthR
448 self.earthRminor = self.earthRmajor
449 elif _isString(earthR) and earthR.startswith('OblateSpheroid'):
450 self.shape_of_earth = earthR
451 scaledearthRmajor = gdtmpl[4]
452 earthRmajorscale = gdtmpl[3]
453 self.earthRmajor = math.pow(10,-earthRmajorscale)*scaledearthRmajor
454 self.earthRmajor = self.earthRmajor*1000.
455 scaledearthRminor = gdtmpl[6]
456 earthRminorscale = gdtmpl[5]
457 self.earthRminor = math.pow(10,-earthRminorscale)*scaledearthRminor
458 self.earthRminor = self.earthRminor*1000.
459 elif _isString(earthR) and earthR.startswith('WGS84'):
460 self.shape_of_earth = earthR
461 self.earthRmajor = 6378137.0
462 self.earthRminor = 6356752.3142
463 elif isinstance(earthR,tuple):
464 self.shape_of_earth = 'OblateSpheroid'
465 self.earthRmajor = earthR[0]
466 self.earthRminor = earthR[1]
467 else:
468 if earthR is not None:
469 self.shape_of_earth = 'Spherical'
470 self.earthRmajor = earthR
471 self.earthRminor = self.earthRmajor
472 if reggrid and gdtnum not in [50,51,52,53,100,120,1000,1200]:
473 self.points_in_x_direction = gdtmpl[7]
474 self.points_in_y_direction = gdtmpl[8]
475 if not reggrid and gdtnum == 40:
476 self.points_in_y_direction = gdtmpl[8]
477 if gdtnum in [0,1,203,205,32768]:
478 scalefact = float(gdtmpl[9])
479 divisor = float(gdtmpl[10])
480 if scalefact == 0: scalefact = 1.
481 if divisor <= 0: divisor = 1.e6
482 self.latitude_first_gridpoint = scalefact*gdtmpl[11]/divisor
483 self.longitude_first_gridpoint = scalefact*gdtmpl[12]/divisor
484 self.latitude_last_gridpoint = scalefact*gdtmpl[14]/divisor
485 self.longitude_last_gridpoint = scalefact*gdtmpl[15]/divisor
486 self.gridlength_in_x_direction = scalefact*gdtmpl[16]/divisor
487 self.gridlength_in_y_direction = scalefact*gdtmpl[17]/divisor
488 if self.latitude_first_gridpoint > self.latitude_last_gridpoint:
489 self.gridlength_in_y_direction = -self.gridlength_in_y_direction
490 if self.longitude_first_gridpoint > self.longitude_last_gridpoint:
491 self.gridlength_in_x_direction = -self.gridlength_in_x_direction
492 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4]
493 if gdtnum == 1:
494 self.latitude_of_southern_pole = scalefact*gdtmpl[19]/divisor
495 self.longitude_of_southern_pole = scalefact*gdtmpl[20]/divisor
496 self.angle_of_pole_rotation = gdtmpl[21]
497 elif gdtnum == 10:
498 self.latitude_first_gridpoint = gdtmpl[9]/1.e6
499 self.longitude_first_gridpoint = gdtmpl[10]/1.e6
500 self.latitude_last_gridpoint = gdtmpl[13]/1.e6
501 self.longitude_last_gridpoint = gdtmpl[14]/1.e6
502 self.gridlength_in_x_direction = gdtmpl[17]/1.e3
503 self.gridlength_in_y_direction= gdtmpl[18]/1.e3
504 self.proj4_lat_ts = gdtmpl[12]/1.e6
505 self.proj4_lon_0 = 0.5*(self.longitude_first_gridpoint+self.longitude_last_gridpoint)
506 self.proj4_proj = 'merc'
507 self.scanmodeflags = _dec2bin(gdtmpl[15])[0:4]
508 elif gdtnum == 20:
509 projflag = _dec2bin(gdtmpl[16])[0]
510 self.latitude_first_gridpoint = gdtmpl[9]/1.e6
511 self.longitude_first_gridpoint = gdtmpl[10]/1.e6
512 self.proj4_lat_ts = gdtmpl[12]/1.e6
513 if projflag == 0:
514 self.proj4_lat_0 = 90
515 elif projflag == 1:
516 self.proj4_lat_0 = -90
517 else:
518 raise ValueError('Invalid projection center flag = %s'%projflag)
519 self.proj4_lon_0 = gdtmpl[13]/1.e6
520 self.gridlength_in_x_direction = gdtmpl[14]/1000.
521 self.gridlength_in_y_direction = gdtmpl[15]/1000.
522 self.proj4_proj = 'stere'
523 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
524 elif gdtnum == 30:
525 self.latitude_first_gridpoint = gdtmpl[9]/1.e6
526 self.longitude_first_gridpoint = gdtmpl[10]/1.e6
527 self.gridlength_in_x_direction = gdtmpl[14]/1000.
528 self.gridlength_in_y_direction = gdtmpl[15]/1000.
529 self.proj4_lat_1 = gdtmpl[18]/1.e6
530 self.proj4_lat_2 = gdtmpl[19]/1.e6
531 self.proj4_lat_0 = gdtmpl[12]/1.e6
532 self.proj4_lon_0 = gdtmpl[13]/1.e6
533 self.proj4_proj = 'lcc'
534 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
535 elif gdtnum == 31:
536 self.latitude_first_gridpoint = gdtmpl[9]/1.e6
537 self.longitude_first_gridpoint = gdtmpl[10]/1.e6
538 self.gridlength_in_x_direction = gdtmpl[14]/1000.
539 self.gridlength_in_y_direction = gdtmpl[15]/1000.
540 self.proj4_lat_1 = gdtmpl[18]/1.e6
541 self.proj4_lat_2 = gdtmpl[19]/1.e6
542 self.proj4_lat_0 = gdtmpl[12]/1.e6
543 self.proj4_lon_0 = gdtmpl[13]/1.e6
544 self.proj4_proj = 'aea'
545 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
546 elif gdtnum == 40 or gdtnum == 41:
547 scalefact = float(gdtmpl[9])
548 divisor = float(gdtmpl[10])
549 if scalefact == 0: scalefact = 1.
550 if divisor <= 0: divisor = 1.e6
551 self.points_between_pole_and_equator = gdtmpl[17]
552 self.latitude_first_gridpoint = scalefact*gdtmpl[11]/divisor
553 self.longitude_first_gridpoint = scalefact*gdtmpl[12]/divisor
554 self.latitude_last_gridpoint = scalefact*gdtmpl[14]/divisor
555 self.longitude_last_gridpoint = scalefact*gdtmpl[15]/divisor
556 if reggrid:
557 self.gridlength_in_x_direction = scalefact*gdtmpl[16]/divisor
558 if self.longitude_first_gridpoint > self.longitude_last_gridpoint:
559 self.gridlength_in_x_direction = -self.gridlength_in_x_direction
560 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4]
561 if gdtnum == 41:
562 self.latitude_of_southern_pole = scalefact*gdtmpl[19]/divisor
563 self.longitude_of_southern_pole = scalefact*gdtmpl[20]/divisor
564 self.angle_of_pole_rotation = gdtmpl[21]
565 elif gdtnum == 50:
566 self.spectral_truncation_parameters = (gdtmpl[0],gdtmpl[1],gdtmpl[2])
567 self.scanmodeflags = [None,None,None,None]
568 elif gdtnum == 90:
569 self.proj4_lat_0 = gdtmpl[9]/1.e6
570 self.proj4_lon_0 = gdtmpl[10]/1.e6
571 self.proj4_h = self.earthRmajor * (gdtmpl[18]/1.e6)
572 dx = gdtmpl[12]
573 dy = gdtmpl[13]
574
575 if self.proj4_lat_0 == 0.:
576 self.proj4_proj = 'geos'
577
578 else:
579 self.proj4_proj = 'nsper'
580 msg = """
581 only geostationary perspective is supported.
582 lat/lon values returned by grid method may be incorrect."""
583 warnings.warn(msg)
584
585 lonmax = 90.-(180./np.pi)*np.arcsin(self.earthRmajor/self.proj4_h)
586
587 latmax = 90.-(180./np.pi)*np.arcsin(self.earthRminor/self.proj4_h)
588
589
590 latmax = int(1000*latmax)/1000.
591 lonmax = int(1000*lonmax)/1000.
592
593 self.proj4_h = self.proj4_h - self.earthRmajor
594
595 P = pyproj.Proj(proj=self.proj4_proj,\
596 a=self.earthRmajor,b=self.earthRminor,\
597 lat_0=0,lon_0=0,h=self.proj4_h)
598 x1,y1 = P(0.,latmax); x2,y2 = P(lonmax,0.)
599 width = 2*x2; height = 2*y1
600 self.gridlength_in_x_direction = width/dx
601 self.gridlength_in_y_direction = height/dy
602 self.scanmodeflags = _dec2bin(gdtmpl[16])[0:4]
603 elif gdtnum == 110:
604 self.proj4_lat_0 = gdtmpl[9]/1.e6
605 self.proj4_lon_0 = gdtmpl[10]/1.e6
606 self.gridlength_in_x_direction = gdtmpl[12]/1000.
607 self.gridlength_in_y_direction = gdtmpl[13]/1000.
608 self.proj4_proj = 'aeqd'
609 self.scanmodeflags = _dec2bin(gdtmpl[15])[0:4]
610 elif gdtnum == 204:
611 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4]
612
613 drtnum = self.data_representation_template_number
614 drtmpl = self.data_representation_template
615 if (drtnum == 2 or drtnum == 3) and drtmpl[6] != 0:
616 self.missing_value = _getieeeint(drtmpl[7])
617 if drtmpl[6] == 2:
618 self.missing_value2 = _getieeeint(drtmpl[8])
619
621 strings = []
622 keys = self.__dict__.keys()
623 keys.sort()
624 for k in keys:
625 if not k.startswith('_'):
626 strings.append('%s = %s\n'%(k,self.__dict__[k]))
627 return ''.join(strings)
628
629 - def data(self,fill_value=9.9692099683868690e+36,masked_array=True,expand=True,order=None):
630 """
631 returns an unpacked data grid. Can also be accomplished with L{values}
632 property.
633
634 @keyword fill_value: missing or masked data is filled with this value
635 (default 9.9692099683868690e+36).
636
637 @keyword masked_array: if True, return masked array if there is bitmap
638 for missing or masked data (default True).
639
640 @keyword expand: if True (default), ECMWF 'reduced' gaussian grids are
641 expanded to regular gaussian grids.
642
643 @keyword order: if 1, linear interpolation is used for expanding reduced
644 gaussian grids. if 0, nearest neighbor interpolation is used. Default
645 is 0 if grid has missing or bitmapped values, 1 otherwise.
646
647 @return: C{B{data}}, a float32 numpy regular or masked array
648 with shape (nlats,lons) containing the requested grid.
649 """
650
651
652 from redtoreg import _redtoreg
653 if not hasattr(self,'scanmodeflags'):
654 raise ValueError('unsupported grid definition template number %s'%self.grid_definition_template_number)
655 else:
656 if self.scanmodeflags[2]:
657 storageorder='F'
658 else:
659 storageorder='C'
660 bitmapflag = self.bitmap_indicator_flag
661 drtnum = self.data_representation_template_number
662 drtmpl = self.data_representation_template
663
664 if order is None:
665 if ((drtnum == 3 or drtnum == 2) and drtmpl[6] != 0) or bitmapflag == 0:
666 order = 0
667 else:
668 order = 1
669 try:
670 f = open(self._grib_filename,'rb')
671 except (TypeError,IOError):
672 f = StringIO(self._grib_filename)
673 f.seek(self._grib_message_byteoffset)
674 gribmsg = f.read(self._grib_message_length)
675 f.close()
676 gdtnum = self.grid_definition_template_number
677 gdtmpl = self.grid_definition_template
678 ndpts = self.number_of_data_points_to_unpack
679 gdsinfo = self.grid_definition_info
680 ngrdpts = gdsinfo[1]
681 ipos = self._section7_byte_offset
682 fld1=g2clib.unpack7(gribmsg,gdtnum,gdtmpl,drtnum,drtmpl,ndpts,ipos,np.empty,storageorder=storageorder)
683
684 if bitmapflag == 0:
685 bitmap=self._bitmap
686 fld = fill_value*np.ones(ngrdpts,'f')
687 np.put(fld,np.nonzero(bitmap),fld1)
688 if masked_array:
689 fld = ma.masked_values(fld,fill_value)
690
691 elif masked_array and hasattr(self,'missing_value'):
692 if hasattr(self, 'missing_value2'):
693 mask = np.logical_or(fld1 == self.missing_value, fld1 == self.missing_value2)
694 else:
695 mask = fld1 == self.missing_value
696 fld = ma.array(fld1,mask=mask)
697 else:
698 fld = fld1
699 nx = None; ny = None
700 if hasattr(self,'points_in_x_direction'):
701 nx = self.points_in_x_direction
702 if hasattr(self,'points_in_y_direction'):
703 ny = self.points_in_y_direction
704 if nx is not None and ny is not None:
705 if ma.isMA(fld):
706 fld = ma.reshape(fld,(ny,nx))
707 else:
708 fld = np.reshape(fld,(ny,nx))
709 else:
710 if gdsinfo[2] and gdtnum == 40:
711 if expand:
712 nx = 2*ny
713 lonsperlat = self.grid_definition_list
714 if ma.isMA(fld):
715 fld = ma.filled(fld)
716 print fld.shape, lonsperlat.dtype, lonsperlat.shape, lonsperlat.sum()
717 fld = _redtoreg(nx, lonsperlat.astype(np.long),\
718 fld.astype(np.double), fill_value)
719 fld = ma.masked_values(fld,fill_value)
720 else:
721 fld = _redtoreg(nx, lonsperlat, fld.astype('f8'),\
722 fill_value)
723
724 if nx is not None and ny is not None:
725
726 if self.scanmodeflags[0]:
727 fldsave = fld.astype('f')
728 fld[:,:] = fldsave[:,::-1]
729
730 if not self.scanmodeflags[1]:
731 fldsave = fld.astype('f')
732 fld[:,:] = fldsave[::-1,:]
733
734
735 if self.scanmodeflags[3]:
736 fldsave = fld.astype('f')
737 fld[1::2,:] = fldsave[1::2,::-1]
738 return fld
739
740 values = property(data)
741
743 """alias for L{grid}"""
744 return self.grid()
745
747 """
748 return lats,lons (in degrees) of grid.
749 currently can handle reg. lat/lon, global gaussian, mercator, stereographic,
750 lambert conformal, albers equal-area, space-view and azimuthal
751 equidistant grids. L{latlons} method does the same thing.
752
753 @return: C{B{lats},B{lons}}, float32 numpy arrays
754 containing latitudes and longitudes of grid (in degrees).
755 """
756 from pygrib import gaulats
757 gdsinfo = self.grid_definition_info
758 gdtnum = self.grid_definition_template_number
759 gdtmpl = self.grid_definition_template
760 reggrid = gdsinfo[2] == 0
761 projparams = {}
762 projparams['a']=self.earthRmajor
763 projparams['b']=self.earthRminor
764 if gdtnum == 0:
765 lon1, lat1 = self.longitude_first_gridpoint, self.latitude_first_gridpoint
766 lon2, lat2 = self.longitude_last_gridpoint, self.latitude_last_gridpoint
767 delon = self.gridlength_in_x_direction
768 delat = self.gridlength_in_y_direction
769 lats = np.arange(lat1,lat2+delat,delat)
770 lons = np.arange(lon1,lon2+delon,delon)
771
772 if self.scanmodeflags[0]:
773 lons = lons[::-1]
774 if not self.scanmodeflags[1]:
775 lats = lats[::-1]
776 projparams['proj'] = 'cyl'
777 lons,lats = np.meshgrid(lons,lats)
778 elif gdtnum == 40:
779 lon1, lat1 = self.longitude_first_gridpoint, self.latitude_first_gridpoint
780 lon2, lat2 = self.longitude_last_gridpoint, self.latitude_last_gridpoint
781 nlats = self.points_in_y_direction
782 if not reggrid:
783 nlons = 2*nlats
784 delon = 360./nlons
785 else:
786 nlons = self.points_in_x_direction
787 delon = self.gridlength_in_x_direction
788 lons = np.arange(lon1,lon2+delon,delon)
789
790 lats = gaulats(nlats)
791 if lat1 < lat2:
792 lats = lats[::-1]
793
794 if self.scanmodeflags[0]:
795 lons = lons[::-1]
796 if not self.scanmodeflags[1]:
797 lats = lats[::-1]
798 projparams['proj'] = 'cyl'
799 lons,lats = np.meshgrid(lons,lats)
800
801 elif gdtnum in [10,20,30,31,110]:
802 nx = self.points_in_x_direction
803 ny = self.points_in_y_direction
804 dx = self.gridlength_in_x_direction
805 dy = self.gridlength_in_y_direction
806 lon1, lat1 = self.longitude_first_gridpoint, self.latitude_first_gridpoint
807 if gdtnum == 10:
808 projparams['lat_ts']=self.proj4_lat_ts
809 projparams['proj']=self.proj4_proj
810 projparams['lon_0']=self.proj4_lon_0
811 pj = pyproj.Proj(projparams)
812 llcrnrx, llcrnry = pj(lon1,lat1)
813 x = llcrnrx+dx*np.arange(nx)
814 y = llcrnry+dy*np.arange(ny)
815 x, y = np.meshgrid(x, y)
816 lons, lats = pj(x, y, inverse=True)
817 elif gdtnum == 20:
818 projparams['lat_ts']=self.proj4_lat_ts
819 projparams['proj']=self.proj4_proj
820 projparams['lat_0']=self.proj4_lat_0
821 projparams['lon_0']=self.proj4_lon_0
822 pj = pyproj.Proj(projparams)
823 llcrnrx, llcrnry = pj(lon1,lat1)
824 x = llcrnrx+dx*np.arange(nx)
825 y = llcrnry+dy*np.arange(ny)
826 x, y = np.meshgrid(x, y)
827 lons, lats = pj(x, y, inverse=True)
828 elif gdtnum in [30,31]:
829 projparams['lat_1']=self.proj4_lat_1
830 projparams['lat_2']=self.proj4_lat_2
831 projparams['proj']=self.proj4_proj
832 projparams['lon_0']=self.proj4_lon_0
833 pj = pyproj.Proj(projparams)
834 llcrnrx, llcrnry = pj(lon1,lat1)
835 x = llcrnrx+dx*np.arange(nx)
836 y = llcrnry+dy*np.arange(ny)
837 x, y = np.meshgrid(x, y)
838 lons, lats = pj(x, y, inverse=True)
839 elif gdtnum == 110:
840 projparams['proj']=self.proj4_proj
841 projparams['lat_0']=self.proj4_lat_0
842 projparams['lon_0']=self.proj4_lon_0
843 pj = pyproj.Proj(projparams)
844 llcrnrx, llcrnry = pj(lon1,lat1)
845 x = llcrnrx+dx*np.arange(nx)
846 y = llcrnry+dy*np.arange(ny)
847 x, y = np.meshgrid(x, y)
848 lons, lats = pj(x, y, inverse=True)
849 elif gdtnum == 90:
850 nx = self.points_in_x_direction
851 ny = self.points_in_y_direction
852 dx = self.gridlength_in_x_direction
853 dy = self.gridlength_in_y_direction
854 projparams['proj']=self.proj4_proj
855 projparams['lon_0']=self.proj4_lon_0
856 projparams['lat_0']=self.proj4_lat_0
857 projparams['h']=self.proj4_h
858 pj = pyproj.Proj(projparams)
859 x = dx*np.indices((ny,nx),'f')[1,:,:]
860 x = x - 0.5*x.max()
861 y = dy*np.indices((ny,nx),'f')[0,:,:]
862 y = y - 0.5*y.max()
863 lons, lats = pj(x,y,inverse=True)
864
865 abslons = np.fabs(lons); abslats = np.fabs(lats)
866 lons = np.where(abslons < 1.e20, lons, 1.e30)
867 lats = np.where(abslats < 1.e20, lats, 1.e30)
868 else:
869 raise ValueError('unsupported grid')
870 self.projparams = projparams
871 return lats.astype('f'), lons.astype('f')
872
874 """
875 Read the contents of a GRIB2 file.
876
877 @param filename: name of GRIB2 file (default, gribmsg=False) or binary string
878 representing a grib message (if gribmsg=True).
879
880 @return: a list of L{Grib2Message} instances representing all of the
881 grib messages in the file. Messages with multiple fields are split
882 into separate messages (so that each L{Grib2Message} instance contains
883 just one data field). The metadata in each GRIB2 message can be
884 accessed via L{Grib2Message} instance variables, the actual data
885 can be read using L{Grib2Message.data}, and the lat/lon values of the grid
886 can be accesses using L{Grib2Message.grid}. If there is only one grib
887 message, just the L{Grib2Message} instance is returned, instead of a list
888 with one element.
889 """
890 if gribmsg:
891 f = StringIO(filename)
892 else:
893 f = open(filename,'rb')
894 nmsg = 0
895
896 disciplines = []
897 startingpos = []
898 msglen = []
899 while 1:
900
901 nbyte = f.tell()
902 while 1:
903 f.seek(nbyte)
904 start = f.read(4).decode('ascii','ignore')
905 if start == '' or start == 'GRIB': break
906 nbyte = nbyte + 1
907 if start == '': break
908
909 startpos = f.tell()-4
910 f.seek(2,1)
911
912 disciplines.append(struct.unpack('>B',f.read(1))[0])
913
914 vers = struct.unpack('>B',f.read(1))[0]
915 if vers != 2:
916 raise IOError('not a GRIB2 file (version number %d)' % vers)
917 lengrib = struct.unpack('>q',f.read(8))[0]
918 msglen.append(lengrib)
919 startingpos.append(startpos)
920
921 f.seek(startpos)
922 gribmsg = f.read(lengrib)
923
924 end = gribmsg[-4:lengrib].decode('ascii','ignore')
925 if end != '7777':
926 raise IOError('partial GRIB message (no "7777" at end)')
927
928 nmsg=nmsg+1
929
930 if nmsg==0:
931 raise IOError('not a GRIB file')
932
933 numfields = []
934 f.seek(0)
935 for n in range(nmsg):
936 f.seek(startingpos[n])
937 gribmsg = f.read(msglen[n])
938 pos = 0
939 numflds = 0
940 while 1:
941 if gribmsg[pos:pos+4].decode('ascii','ignore') == 'GRIB':
942 sectnum = 0
943 lensect = 16
944 elif gribmsg[pos:pos+4].decode('ascii','ignore') == '7777':
945 break
946 else:
947 lensect = struct.unpack('>i',gribmsg[pos:pos+4])[0]
948 sectnum = struct.unpack('>B',gribmsg[pos+4:pos+5])[0]
949 if sectnum == 4: numflds=numflds+1
950
951 pos = pos + lensect
952
953
954 numfields.append(numflds)
955
956 gdtnum = []
957 gdtmpl = []
958 gdeflist = []
959 gdsinfo = []
960 pdtmpl = []
961 pdtnum = []
962 coordlist = []
963 drtmpl = []
964 drtnum = []
965 ndpts = []
966 bitmapflag = []
967 bitmap = []
968 pos7 = []
969 localsxn = []
970 msgstart = []
971 msglength = []
972 message = []
973 identsect = []
974 discipline = []
975 for n in range(nmsg):
976 spos = startingpos[n]
977 lengrib = msglen[n]
978
979 f.seek(spos)
980 gribmsg = f.read(lengrib)
981 discipl = disciplines[n]
982 lensect0 = 16
983
984
985
986
987
988
989
990 idsect,pos = g2clib.unpack1(gribmsg,lensect0,np.empty)
991
992 gdtnums = []
993 gdtmpls = []
994 gdeflists = []
995 gdsinfos = []
996 pdtmpls = []
997 coordlists = []
998 pdtnums = []
999 drtmpls = []
1000 drtnums = []
1001 ndptslist = []
1002 bitmapflags = []
1003 bitmaps = []
1004 sxn7pos = []
1005 localsxns = []
1006 while 1:
1007
1008 if gribmsg[pos:pos+4].decode('ascii','ignore') == '7777': break
1009 lensect = struct.unpack('>i',gribmsg[pos:pos+4])[0]
1010 sectnum = struct.unpack('>B',gribmsg[pos+4:pos+5])[0]
1011
1012 if sectnum == 2:
1013
1014
1015
1016
1017 localsxns.append(gribmsg[pos+5:pos+lensect])
1018 pos = pos + lensect
1019
1020 elif sectnum == 3:
1021 gds,gdtempl,deflist,pos = g2clib.unpack3(gribmsg,pos,np.empty)
1022 gdtnums.append(gds[4])
1023 gdtmpls.append(gdtempl)
1024 gdeflists.append(deflist)
1025 gdsinfos.append(gds)
1026
1027 elif sectnum == 4:
1028 pdtempl,pdtn,coordlst,pos = g2clib.unpack4(gribmsg,pos,np.empty)
1029 pdtmpls.append(pdtempl)
1030 coordlists.append(coordlst)
1031 pdtnums.append(pdtn)
1032
1033 elif sectnum == 5:
1034 drtempl,drtn,npts,pos = g2clib.unpack5(gribmsg,pos,np.empty)
1035 drtmpls.append(drtempl)
1036 drtnums.append(drtn)
1037 ndptslist.append(npts)
1038
1039 elif sectnum == 6:
1040 bmap,bmapflag = g2clib.unpack6(gribmsg,gds[1],pos,np.empty)
1041
1042 if bmapflag == 0:
1043 bitmaps.append(bmap.astype('b'))
1044
1045 elif bmapflag == 254:
1046 bmapflag = 0
1047 for bmp in bitmaps[::-1]:
1048 if bmp is not None: bitmaps.append(bmp)
1049 else:
1050 bitmaps.append(None)
1051 bitmapflags.append(bmapflag)
1052 pos = pos + lensect
1053
1054
1055 else:
1056 if sectnum != 7:
1057 msg = 'unknown section = %i' % sectnum
1058 raise ValueError(msg)
1059 sxn7pos.append(pos)
1060 pos = pos + lensect
1061
1062 gdtnum.append(_repeatlast(numfields[n],gdtnums))
1063 gdtmpl.append(_repeatlast(numfields[n],gdtmpls))
1064 gdeflist.append(_repeatlast(numfields[n],gdeflists))
1065 gdsinfo.append(_repeatlast(numfields[n],gdsinfos))
1066 pdtmpl.append(_repeatlast(numfields[n],pdtmpls))
1067 pdtnum.append(_repeatlast(numfields[n],pdtnums))
1068 coordlist.append(_repeatlast(numfields[n],coordlists))
1069 drtmpl.append(_repeatlast(numfields[n],drtmpls))
1070 drtnum.append(_repeatlast(numfields[n],drtnums))
1071 ndpts.append(_repeatlast(numfields[n],ndptslist))
1072 bitmapflag.append(_repeatlast(numfields[n],bitmapflags))
1073 bitmap.append(_repeatlast(numfields[n],bitmaps))
1074 pos7.append(_repeatlast(numfields[n],sxn7pos))
1075 if len(localsxns) == 0:
1076 localsxns = [None]
1077 localsxn.append(_repeatlast(numfields[n],localsxns))
1078 msgstart.append(_repeatlast(numfields[n],[spos]))
1079 msglength.append(_repeatlast(numfields[n],[lengrib]))
1080 identsect.append(_repeatlast(numfields[n],[idsect]))
1081 discipline.append(_repeatlast(numfields[n],[discipl]))
1082
1083 gdtnum = _flatten(gdtnum)
1084 gdtmpl = _flatten(gdtmpl)
1085 gdeflist = _flatten(gdeflist)
1086 gdsinfo = _flatten(gdsinfo)
1087 pdtmpl = _flatten(pdtmpl)
1088 pdtnum = _flatten(pdtnum)
1089 coordlist = _flatten(coordlist)
1090 drtmpl = _flatten(drtmpl)
1091 drtnum = _flatten(drtnum)
1092 ndpts = _flatten(ndpts)
1093 bitmapflag = _flatten(bitmapflag)
1094 bitmap = _flatten(bitmap)
1095 pos7 = _flatten(pos7)
1096 localsxn = _flatten(localsxn)
1097 msgstart = _flatten(msgstart)
1098 msglength = _flatten(msglength)
1099 identsect = _flatten(identsect)
1100 discipline = _flatten(discipline)
1101
1102 gribs = []
1103 for n in range(len(msgstart)):
1104 kwargs = {}
1105 kwargs['originating_center']=_table0[identsect[n][0]][0]
1106 wmo_code = _table0[identsect[n][0]][1]
1107 if wmo_code is not None:
1108 kwargs['center_wmo_code']=wmo_code
1109 kwargs['grid_definition_template_number']=gdtnum[n]
1110 kwargs['grid_definition_template']=gdtmpl[n]
1111 if gdeflist[n] != []:
1112 kwargs['grid_definition_list']=gdeflist[n]
1113 kwargs['grid_definition_info']=gdsinfo[n]
1114 kwargs['discipline_code']=discipline[n]
1115 kwargs['product_definition_template_number']=pdtnum[n]
1116 kwargs['product_definition_template']=pdtmpl[n]
1117 kwargs['data_representation_template_number']=drtnum[n]
1118 kwargs['data_representation_template']=drtmpl[n]
1119 if coordlist[n] != []:
1120 kwargs['extra_vertical_coordinate_info']=coordlist[n]
1121 kwargs['number_of_data_points_to_unpack']=ndpts[n]
1122 kwargs['bitmap_indicator_flag']=bitmapflag[n]
1123 if bitmap[n] is not []:
1124 kwargs['_bitmap']=bitmap[n]
1125 kwargs['_section7_byte_offset']=pos7[n]
1126 kwargs['_grib_message_byteoffset']=msgstart[n]
1127 kwargs['_grib_message_length']=msglength[n]
1128 kwargs['_grib_filename']=filename
1129 kwargs['identification_section']=identsect[n]
1130 kwargs['_grib_message_number']=n+1
1131 if localsxn[n] is not None:
1132 kwargs['has_local_use_section'] = True
1133 kwargs['_local_use_section']=localsxn[n]
1134 else:
1135 kwargs['has_local_use_section'] = False
1136 gribs.append(Grib2Message(**kwargs))
1137 f.close()
1138 if len(gribs) == 1:
1139 return gribs[0]
1140 else:
1141 return gribs
1142
1143 -def dump(filename, grbs):
1144 """
1145 write the given L{Grib2Message} instances to a grib file.
1146
1147 @param filename: file to write grib data to.
1148 @param grbs: a list of L{Grib2Message} instances.
1149 """
1150 gribfile = open(filename,'wb')
1151 for grb in grbs:
1152 try:
1153 f = open(grb._grib_filename,'rb')
1154 except TypeError:
1155 f = StringIO(grb._grib_filename)
1156 f.seek(grb._grib_message_byteoffset)
1157 gribmsg = f.read(grb._grib_message_length)
1158 f.close()
1159 gribfile.write(gribmsg)
1160 gribfile.close()
1161
1162
1163
1165 """return yyyy,mm,dd,min,ss from section 1"""
1166 yyyy=idsect[5]
1167 mm=idsect[6]
1168 dd=idsect[7]
1169 hh=idsect[8]
1170 min=idsect[9]
1171 ss=idsect[10]
1172 return yyyy,mm,dd,hh,min,ss
1173
1175 """unpack section 1 given starting point in bytes
1176 used to test pyrex interface to g2_unpack1"""
1177 idsect = []
1178 pos = pos + 5
1179 idsect.append(struct.unpack('>h',gribmsg[pos:pos+2])[0])
1180 pos = pos + 2
1181 idsect.append(struct.unpack('>h',gribmsg[pos:pos+2])[0])
1182 pos = pos + 2
1183 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1184 pos = pos + 1
1185 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1186 pos = pos + 1
1187 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1188 pos = pos + 1
1189 idsect.append(struct.unpack('>h',gribmsg[pos:pos+2])[0])
1190 pos = pos + 2
1191 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1192 pos = pos + 1
1193 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1194 pos = pos + 1
1195 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1196 pos = pos + 1
1197 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1198 pos = pos + 1
1199 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1200 pos = pos + 1
1201 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1202 pos = pos + 1
1203 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1204 pos = pos + 1
1205 return np.array(idsect,'i'),pos
1206
1208 """repeat last item in listin, until len(listin) = numfields"""
1209 if len(listin) < numfields:
1210 last = listin[-1]
1211 for n in range(len(listin),numfields):
1212 listin.append(last)
1213 return listin
1214
1216 try:
1217 flist = functools.reduce(operator.add,lst)
1218 except NameError:
1219 import functools
1220 flist = functools.reduce(operator.add,lst)
1221 return flist
1222
1223
1225 """
1226 Class for encoding data into a GRIB2 message.
1227 - Creating a class instance (L{__init__}) initializes the message and adds
1228 sections 0 and 1 (the indicator and identification sections),
1229 - method L{addgrid} adds a grid definition (section 3) to the messsage.
1230 - method L{addfield} adds sections 4-7 to the message (the product
1231 definition, data representation, bitmap and data sections).
1232 - method L{end} adds the end section (section 8) and terminates the message.
1233
1234
1235 A GRIB Edition 2 message is a machine independent format for storing
1236 one or more gridded data fields. Each GRIB2 message consists of the
1237 following sections:
1238 - SECTION 0: Indicator Section - only one per message
1239 - SECTION 1: Identification Section - only one per message
1240 - SECTION 2: (Local Use Section) - optional
1241 - SECTION 3: Grid Definition Section
1242 - SECTION 4: Product Definition Section
1243 - SECTION 5: Data Representation Section
1244 - SECTION 6: Bit-map Section
1245 - SECTION 7: Data Section
1246 - SECTION 8: End Section
1247
1248 Sequences of GRIB sections 2 to 7, 3 to 7, or sections 4 to 7 may be repeated
1249 within a single GRIB message. All sections within such repeated sequences
1250 must be present and shall appear in the numerical order noted above.
1251 Unrepeated sections remain in effect until redefined.
1252
1253 Note: Writing section 2 (the 'local use section') is
1254 not yet supported.
1255
1256 @ivar msg: A binary string containing the GRIB2 message.
1257 After the message has been terminated by calling
1258 the L{end} method, this string can be written to a file.
1259 """
1260
1261 - def __init__(self, discipline, idsect):
1262 """
1263 create a Grib2Enecode class instance given the GRIB2 discipline
1264 parameter and the identification section (sections 0 and 1).
1265
1266 The GRIB2 message is stored as a binary string in instance variable L{msg}.
1267
1268 L{addgrid}, L{addfield} and L{end} class methods must be called to complete
1269 the GRIB2 message.
1270
1271 @param discipline: Discipline or GRIB Master Table Number (Code Table 0.0).
1272 (0 for meteorlogical, 1 for hydrological, 2 for land surface, 3 for space,
1273 10 for oceanographic products).
1274
1275 @param idsect: Sequence containing identification section (section 1).
1276 - idsect[0]=Id of orginating centre (Common Code
1277 U{Table C-1<http://www.nws.noaa.gov/tg/GRIB_C1.htm>})
1278 - idsect[1]=Id of orginating sub-centre (local table)
1279 - idsect[2]=GRIB Master Tables Version Number (Code
1280 U{Table 1.0
1281 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-0.shtml>})
1282 - idsect[3]=GRIB Local Tables Version Number (Code
1283 U{Table 1.1
1284 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-1.shtml>})
1285 - idsect[4]=Significance of Reference Time (Code
1286 U{Table 1.2
1287 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-2.shtml>})
1288 - idsect[5]=Reference Time - Year (4 digits)
1289 - idsect[6]=Reference Time - Month
1290 - idsect[7]=Reference Time - Day
1291 - idsect[8]=Reference Time - Hour
1292 - idsect[9]=Reference Time - Minute
1293 - idsect[10]=Reference Time - Second
1294 - idsect[11]=Production status of data (Code
1295 U{Table 1.3
1296 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-3.shtml>})
1297 - idsect[12]=Type of processed data (Code
1298 U{Table
1299 1.4<http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-4.shtml>})
1300 """
1301 self.msg,msglen=g2clib.grib2_create(np.array([discipline,2],np.int32),np.array(idsect,np.int32))
1302
1303 - def addgrid(self,gdsinfo,gdtmpl,deflist=None):
1304 """
1305 Add a grid definition section (section 3) to the GRIB2 message.
1306
1307 @param gdsinfo: Sequence containing information needed for the grid definition section.
1308 - gdsinfo[0] = Source of grid definition (see Code
1309 U{Table 3.0
1310 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-0.shtml>})
1311 - gdsinfo[1] = Number of grid points in the defined grid.
1312 - gdsinfo[2] = Number of octets needed for each additional grid points defn.
1313 Used to define number of points in each row for non-reg grids (=0 for
1314 regular grid).
1315 - gdsinfo[3] = Interp. of list for optional points defn (Code
1316 U{Table 3.11
1317 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-11.shtml>})
1318 - gdsinfo[4] = Grid Definition Template Number (Code
1319 U{Table 3.1
1320 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-1.shtml>})
1321
1322 @param gdtmpl: Contains the data values for the specified Grid Definition
1323 Template ( NN=gdsinfo[4] ). Each element of this integer
1324 array contains an entry (in the order specified) of Grid
1325 Definition Template 3.NN
1326
1327 @param deflist: (Used if gdsinfo[2] != 0) Sequence containing the
1328 number of grid points contained in each row (or column)
1329 of a non-regular grid.
1330 """
1331 if deflist is not None:
1332 dflist = np.array(deflist,'i')
1333 else:
1334 dflist = None
1335 self.scanmodeflags = None
1336 gdtnum = gdsinfo[4]
1337 if gdtnum in [0,1,2,3,40,41,42,43,44,203,205,32768,32769]:
1338 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4]
1339 elif gdtnum == 10:
1340 self.scanmodeflags = _dec2bin(gdtmpl[15])[0:4]
1341 elif gdtnum == 20:
1342 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
1343 elif gdtnum == 30:
1344 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
1345 elif gdtnum == 31:
1346 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
1347 elif gdtnum == 90:
1348 self.scanmodeflags = _dec2bin(gdtmpl[16])[0:4]
1349 elif gdtnum == 110:
1350 self.scanmodeflags = _dec2bin(gdtmpl[15])[0:4]
1351 elif gdtnum == 120:
1352 self.scanmodeflags = _dec2bin(gdtmpl[6])[0:4]
1353 elif gdtnum == 204:
1354 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4]
1355 elif gdtnum in [1000,1100]:
1356 self.scanmodeflags = _dec2bin(gdtmpl[12])[0:4]
1357 self.msg,msglen=g2clib.grib2_addgrid(self.msg,np.array(gdsinfo,'i'),np.array(gdtmpl,'i'),dflist)
1358
1359 - def addfield(self,pdtnum,pdtmpl,drtnum,drtmpl,field,coordlist=None):
1360 """
1361 Add a product definition section, data representation section,
1362 bitmap section and data section to the GRIB2 message (sections 4-7).
1363 Must be called after grid definition section is created with L{addgrid}.
1364
1365 @param pdtnum: Product Definition Template Number (see Code U{Table
1366 4.0<http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-0.shtml>})
1367
1368 @param pdtmpl: Sequence with the data values for the specified Product Definition
1369 Template (N=pdtnum). Each element of this integer
1370 array contains an entry (in the order specified) of Product
1371 Definition Template 4.N
1372
1373 @param drtnum: Data Representation Template Number (see Code
1374 U{Table 5.0
1375 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table5-0.shtml>})
1376
1377 @param drtmpl: Sequence with the data values for the specified Data Representation
1378 Template (N=drtnum). Each element of this integer
1379 array contains an entry (in the order specified) of Data
1380 Representation Template 5.N
1381 Note that some values in this template (eg. reference
1382 values, number of bits, etc...) may be changed by the
1383 data packing algorithms.
1384 Use this to specify scaling factors and order of
1385 spatial differencing, if desired.
1386
1387 @param field: numpy array of data points to pack.
1388 If field is a masked array, then a bitmap is created from
1389 the mask.
1390
1391 @param coordlist: Sequence containing floating point values intended to document
1392 the vertical discretization with model data
1393 on hybrid coordinate vertical levels. Default None.
1394 """
1395 if not hasattr(self,'scanmodeflags'):
1396 raise ValueError('addgrid must be called before addfield')
1397
1398
1399 if self.scanmodeflags is not None:
1400 if self.scanmodeflags[0]:
1401
1402 fieldsave = field.astype('f')
1403 field[:,:] = fieldsave[:,::-1]
1404
1405 if not self.scanmodeflags[1]:
1406 fieldsave = field.astype('f')
1407 field[:,:] = fieldsave[::-1,:]
1408
1409
1410 if self.scanmodeflags[3]:
1411 fieldsave = field.astype('f')
1412 field[1::2,:] = fieldsave[1::2,::-1]
1413 fld = field.astype('f')
1414 if ma.isMA(field):
1415 bmap = 1 - np.ravel(field.mask.astype('i'))
1416 bitmapflag = 0
1417 else:
1418 bitmapflag = 255
1419 bmap = None
1420 if coordlist is not None:
1421 crdlist = np.array(coordlist,'f')
1422 else:
1423 crdlist = None
1424 self.msg,msglen=g2clib.grib2_addfield(self.msg,pdtnum,np.array(pdtmpl,'i'),crdlist,drtnum,np.array(drtmpl,'i'),np.ravel(fld),bitmapflag,bmap)
1425
1427 """
1428 Add an end section (section 8) to the GRIB2 message.
1429 A GRIB2 message is not complete without an end section.
1430 Once an end section is added, the GRIB2 message can be
1431 output to a file.
1432 """
1433 self.msg,msglen=g2clib.grib2_end(self.msg)
1434