1 __version__ = '2.0.0'
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 fld = _redtoreg(nx, lonsperlat.astype(np.long),\
717 fld.astype(np.double), fill_value)
718 fld = ma.masked_values(fld,fill_value)
719 else:
720 fld = _redtoreg(nx, lonsperlat.astype(np.long),\
721 fld.astype(np.double), fill_value)
722
723 if nx is not None and ny is not None:
724
725
726
727
728
729
730
731
732
733
734 if self.scanmodeflags[3]:
735 fldsave = fld.astype('f')
736 fld[1::2,:] = fldsave[1::2,::-1]
737 return fld
738
739 values = property(data)
740
742 """alias for L{grid}"""
743 return self.grid()
744
746 """
747 return lats,lons (in degrees) of grid.
748 currently can handle reg. lat/lon, global gaussian, mercator, stereographic,
749 lambert conformal, albers equal-area, space-view and azimuthal
750 equidistant grids. L{latlons} method does the same thing.
751
752 @return: C{B{lats},B{lons}}, float32 numpy arrays
753 containing latitudes and longitudes of grid (in degrees).
754 """
755 from pygrib import gaulats
756 gdsinfo = self.grid_definition_info
757 gdtnum = self.grid_definition_template_number
758 gdtmpl = self.grid_definition_template
759 reggrid = gdsinfo[2] == 0
760 projparams = {}
761 projparams['a']=self.earthRmajor
762 projparams['b']=self.earthRminor
763 if gdtnum == 0:
764 lon1, lat1 = self.longitude_first_gridpoint, self.latitude_first_gridpoint
765 lon2, lat2 = self.longitude_last_gridpoint, self.latitude_last_gridpoint
766 delon = self.gridlength_in_x_direction
767 delat = self.gridlength_in_y_direction
768 lats = np.arange(lat1,lat2+delat,delat)
769 lons = np.arange(lon1,lon2+delon,delon)
770
771
772
773
774
775 projparams['proj'] = 'cyl'
776 lons,lats = np.meshgrid(lons,lats)
777 elif gdtnum == 40:
778 lon1, lat1 = self.longitude_first_gridpoint, self.latitude_first_gridpoint
779 lon2, lat2 = self.longitude_last_gridpoint, self.latitude_last_gridpoint
780 nlats = self.points_in_y_direction
781 if not reggrid:
782 nlons = 2*nlats
783 delon = 360./nlons
784 else:
785 nlons = self.points_in_x_direction
786 delon = self.gridlength_in_x_direction
787 lons = np.arange(lon1,lon2+delon,delon)
788
789 lats = gaulats(nlats)
790 if lat1 < lat2:
791 lats = lats[::-1]
792
793
794
795
796
797 projparams['proj'] = 'cyl'
798 lons,lats = np.meshgrid(lons,lats)
799
800 elif gdtnum in [10,20,30,31,110]:
801 nx = self.points_in_x_direction
802 ny = self.points_in_y_direction
803 dx = self.gridlength_in_x_direction
804 dy = self.gridlength_in_y_direction
805 lon1, lat1 = self.longitude_first_gridpoint, self.latitude_first_gridpoint
806 if gdtnum == 10:
807 projparams['lat_ts']=self.proj4_lat_ts
808 projparams['proj']=self.proj4_proj
809 projparams['lon_0']=self.proj4_lon_0
810 pj = pyproj.Proj(projparams)
811 llcrnrx, llcrnry = pj(lon1,lat1)
812 x = llcrnrx+dx*np.arange(nx)
813 y = llcrnry+dy*np.arange(ny)
814 x, y = np.meshgrid(x, y)
815 lons, lats = pj(x, y, inverse=True)
816 elif gdtnum == 20:
817 projparams['lat_ts']=self.proj4_lat_ts
818 projparams['proj']=self.proj4_proj
819 projparams['lat_0']=self.proj4_lat_0
820 projparams['lon_0']=self.proj4_lon_0
821 pj = pyproj.Proj(projparams)
822 llcrnrx, llcrnry = pj(lon1,lat1)
823 x = llcrnrx+dx*np.arange(nx)
824 y = llcrnry+dy*np.arange(ny)
825 x, y = np.meshgrid(x, y)
826 lons, lats = pj(x, y, inverse=True)
827 elif gdtnum in [30,31]:
828 projparams['lat_1']=self.proj4_lat_1
829 projparams['lat_2']=self.proj4_lat_2
830 projparams['proj']=self.proj4_proj
831 projparams['lon_0']=self.proj4_lon_0
832 pj = pyproj.Proj(projparams)
833 llcrnrx, llcrnry = pj(lon1,lat1)
834 x = llcrnrx+dx*np.arange(nx)
835 y = llcrnry+dy*np.arange(ny)
836 x, y = np.meshgrid(x, y)
837 lons, lats = pj(x, y, inverse=True)
838 elif gdtnum == 110:
839 projparams['proj']=self.proj4_proj
840 projparams['lat_0']=self.proj4_lat_0
841 projparams['lon_0']=self.proj4_lon_0
842 pj = pyproj.Proj(projparams)
843 llcrnrx, llcrnry = pj(lon1,lat1)
844 x = llcrnrx+dx*np.arange(nx)
845 y = llcrnry+dy*np.arange(ny)
846 x, y = np.meshgrid(x, y)
847 lons, lats = pj(x, y, inverse=True)
848 elif gdtnum == 90:
849 nx = self.points_in_x_direction
850 ny = self.points_in_y_direction
851 dx = self.gridlength_in_x_direction
852 dy = self.gridlength_in_y_direction
853 projparams['proj']=self.proj4_proj
854 projparams['lon_0']=self.proj4_lon_0
855 projparams['lat_0']=self.proj4_lat_0
856 projparams['h']=self.proj4_h
857 pj = pyproj.Proj(projparams)
858 x = dx*np.indices((ny,nx),'f')[1,:,:]
859 x = x - 0.5*x.max()
860 y = dy*np.indices((ny,nx),'f')[0,:,:]
861 y = y - 0.5*y.max()
862 lons, lats = pj(x,y,inverse=True)
863
864 abslons = np.fabs(lons); abslats = np.fabs(lats)
865 lons = np.where(abslons < 1.e20, lons, 1.e30)
866 lats = np.where(abslats < 1.e20, lats, 1.e30)
867 else:
868 raise ValueError('unsupported grid')
869 self.projparams = projparams
870 return lats.astype('f'), lons.astype('f')
871
873 """
874 Read the contents of a GRIB2 file.
875
876 @param filename: name of GRIB2 file (default, gribmsg=False) or binary string
877 representing a grib message (if gribmsg=True).
878
879 @return: a list of L{Grib2Message} instances representing all of the
880 grib messages in the file. Messages with multiple fields are split
881 into separate messages (so that each L{Grib2Message} instance contains
882 just one data field). The metadata in each GRIB2 message can be
883 accessed via L{Grib2Message} instance variables, the actual data
884 can be read using L{Grib2Message.data}, and the lat/lon values of the grid
885 can be accesses using L{Grib2Message.grid}. If there is only one grib
886 message, just the L{Grib2Message} instance is returned, instead of a list
887 with one element.
888 """
889 if gribmsg:
890 f = StringIO(filename)
891 else:
892 f = open(filename,'rb')
893 nmsg = 0
894
895 disciplines = []
896 startingpos = []
897 msglen = []
898 while 1:
899
900 nbyte = f.tell()
901 while 1:
902 f.seek(nbyte)
903 start = f.read(4).decode('ascii','ignore')
904 if start == '' or start == 'GRIB': break
905 nbyte = nbyte + 1
906 if start == '': break
907
908 startpos = f.tell()-4
909 f.seek(2,1)
910
911 disciplines.append(struct.unpack('>B',f.read(1))[0])
912
913 vers = struct.unpack('>B',f.read(1))[0]
914 if vers != 2:
915 raise IOError('not a GRIB2 file (version number %d)' % vers)
916 lengrib = struct.unpack('>q',f.read(8))[0]
917 msglen.append(lengrib)
918 startingpos.append(startpos)
919
920 f.seek(startpos)
921 gribmsg = f.read(lengrib)
922
923 end = gribmsg[-4:lengrib].decode('ascii','ignore')
924 if end != '7777':
925 raise IOError('partial GRIB message (no "7777" at end)')
926
927 nmsg=nmsg+1
928
929 if nmsg==0:
930 raise IOError('not a GRIB file')
931
932 numfields = []
933 f.seek(0)
934 for n in range(nmsg):
935 f.seek(startingpos[n])
936 gribmsg = f.read(msglen[n])
937 pos = 0
938 numflds = 0
939 while 1:
940 if gribmsg[pos:pos+4].decode('ascii','ignore') == 'GRIB':
941 sectnum = 0
942 lensect = 16
943 elif gribmsg[pos:pos+4].decode('ascii','ignore') == '7777':
944 break
945 else:
946 lensect = struct.unpack('>i',gribmsg[pos:pos+4])[0]
947 sectnum = struct.unpack('>B',gribmsg[pos+4:pos+5])[0]
948 if sectnum == 4: numflds=numflds+1
949
950 pos = pos + lensect
951
952
953 numfields.append(numflds)
954
955 gdtnum = []
956 gdtmpl = []
957 gdeflist = []
958 gdsinfo = []
959 pdtmpl = []
960 pdtnum = []
961 coordlist = []
962 drtmpl = []
963 drtnum = []
964 ndpts = []
965 bitmapflag = []
966 bitmap = []
967 pos7 = []
968 localsxn = []
969 msgstart = []
970 msglength = []
971 message = []
972 identsect = []
973 discipline = []
974 for n in range(nmsg):
975 spos = startingpos[n]
976 lengrib = msglen[n]
977
978 f.seek(spos)
979 gribmsg = f.read(lengrib)
980 discipl = disciplines[n]
981 lensect0 = 16
982
983
984
985
986
987
988
989 idsect,pos = g2clib.unpack1(gribmsg,lensect0,np.empty)
990
991 gdtnums = []
992 gdtmpls = []
993 gdeflists = []
994 gdsinfos = []
995 pdtmpls = []
996 coordlists = []
997 pdtnums = []
998 drtmpls = []
999 drtnums = []
1000 ndptslist = []
1001 bitmapflags = []
1002 bitmaps = []
1003 sxn7pos = []
1004 localsxns = []
1005 while 1:
1006
1007 if gribmsg[pos:pos+4].decode('ascii','ignore') == '7777': break
1008 lensect = struct.unpack('>i',gribmsg[pos:pos+4])[0]
1009 sectnum = struct.unpack('>B',gribmsg[pos+4:pos+5])[0]
1010
1011 if sectnum == 2:
1012
1013
1014
1015
1016 localsxns.append(gribmsg[pos+5:pos+lensect])
1017 pos = pos + lensect
1018
1019 elif sectnum == 3:
1020 gds,gdtempl,deflist,pos = g2clib.unpack3(gribmsg,pos,np.empty)
1021 gdtnums.append(gds[4])
1022 gdtmpls.append(gdtempl)
1023 gdeflists.append(deflist)
1024 gdsinfos.append(gds)
1025
1026 elif sectnum == 4:
1027 pdtempl,pdtn,coordlst,pos = g2clib.unpack4(gribmsg,pos,np.empty)
1028 pdtmpls.append(pdtempl)
1029 coordlists.append(coordlst)
1030 pdtnums.append(pdtn)
1031
1032 elif sectnum == 5:
1033 drtempl,drtn,npts,pos = g2clib.unpack5(gribmsg,pos,np.empty)
1034 drtmpls.append(drtempl)
1035 drtnums.append(drtn)
1036 ndptslist.append(npts)
1037
1038 elif sectnum == 6:
1039 bmap,bmapflag = g2clib.unpack6(gribmsg,gds[1],pos,np.empty)
1040
1041 if bmapflag == 0:
1042 bitmaps.append(bmap.astype('b'))
1043
1044 elif bmapflag == 254:
1045 bmapflag = 0
1046 for bmp in bitmaps[::-1]:
1047 if bmp is not None: bitmaps.append(bmp)
1048 else:
1049 bitmaps.append(None)
1050 bitmapflags.append(bmapflag)
1051 pos = pos + lensect
1052
1053
1054 else:
1055 if sectnum != 7:
1056 msg = 'unknown section = %i' % sectnum
1057 raise ValueError(msg)
1058 sxn7pos.append(pos)
1059 pos = pos + lensect
1060
1061 gdtnum.append(_repeatlast(numfields[n],gdtnums))
1062 gdtmpl.append(_repeatlast(numfields[n],gdtmpls))
1063 gdeflist.append(_repeatlast(numfields[n],gdeflists))
1064 gdsinfo.append(_repeatlast(numfields[n],gdsinfos))
1065 pdtmpl.append(_repeatlast(numfields[n],pdtmpls))
1066 pdtnum.append(_repeatlast(numfields[n],pdtnums))
1067 coordlist.append(_repeatlast(numfields[n],coordlists))
1068 drtmpl.append(_repeatlast(numfields[n],drtmpls))
1069 drtnum.append(_repeatlast(numfields[n],drtnums))
1070 ndpts.append(_repeatlast(numfields[n],ndptslist))
1071 bitmapflag.append(_repeatlast(numfields[n],bitmapflags))
1072 bitmap.append(_repeatlast(numfields[n],bitmaps))
1073 pos7.append(_repeatlast(numfields[n],sxn7pos))
1074 if len(localsxns) == 0:
1075 localsxns = [None]
1076 localsxn.append(_repeatlast(numfields[n],localsxns))
1077 msgstart.append(_repeatlast(numfields[n],[spos]))
1078 msglength.append(_repeatlast(numfields[n],[lengrib]))
1079 identsect.append(_repeatlast(numfields[n],[idsect]))
1080 discipline.append(_repeatlast(numfields[n],[discipl]))
1081
1082 gdtnum = _flatten(gdtnum)
1083 gdtmpl = _flatten(gdtmpl)
1084 gdeflist = _flatten(gdeflist)
1085 gdsinfo = _flatten(gdsinfo)
1086 pdtmpl = _flatten(pdtmpl)
1087 pdtnum = _flatten(pdtnum)
1088 coordlist = _flatten(coordlist)
1089 drtmpl = _flatten(drtmpl)
1090 drtnum = _flatten(drtnum)
1091 ndpts = _flatten(ndpts)
1092 bitmapflag = _flatten(bitmapflag)
1093 bitmap = _flatten(bitmap)
1094 pos7 = _flatten(pos7)
1095 localsxn = _flatten(localsxn)
1096 msgstart = _flatten(msgstart)
1097 msglength = _flatten(msglength)
1098 identsect = _flatten(identsect)
1099 discipline = _flatten(discipline)
1100
1101 gribs = []
1102 for n in range(len(msgstart)):
1103 kwargs = {}
1104 kwargs['originating_center']=_table0[identsect[n][0]][0]
1105 wmo_code = _table0[identsect[n][0]][1]
1106 if wmo_code is not None:
1107 kwargs['center_wmo_code']=wmo_code
1108 kwargs['grid_definition_template_number']=gdtnum[n]
1109 kwargs['grid_definition_template']=gdtmpl[n]
1110 if gdeflist[n] != []:
1111 kwargs['grid_definition_list']=gdeflist[n]
1112 kwargs['grid_definition_info']=gdsinfo[n]
1113 kwargs['discipline_code']=discipline[n]
1114 kwargs['product_definition_template_number']=pdtnum[n]
1115 kwargs['product_definition_template']=pdtmpl[n]
1116 kwargs['data_representation_template_number']=drtnum[n]
1117 kwargs['data_representation_template']=drtmpl[n]
1118 if coordlist[n] != []:
1119 kwargs['extra_vertical_coordinate_info']=coordlist[n]
1120 kwargs['number_of_data_points_to_unpack']=ndpts[n]
1121 kwargs['bitmap_indicator_flag']=bitmapflag[n]
1122 if bitmap[n] is not []:
1123 kwargs['_bitmap']=bitmap[n]
1124 kwargs['_section7_byte_offset']=pos7[n]
1125 kwargs['_grib_message_byteoffset']=msgstart[n]
1126 kwargs['_grib_message_length']=msglength[n]
1127 kwargs['_grib_filename']=filename
1128 kwargs['identification_section']=identsect[n]
1129 kwargs['_grib_message_number']=n+1
1130 if localsxn[n] is not None:
1131 kwargs['has_local_use_section'] = True
1132 kwargs['_local_use_section']=localsxn[n]
1133 else:
1134 kwargs['has_local_use_section'] = False
1135 gribs.append(Grib2Message(**kwargs))
1136 f.close()
1137 if len(gribs) == 1:
1138 return gribs[0]
1139 else:
1140 return gribs
1141
1142 -def dump(filename, grbs):
1143 """
1144 write the given L{Grib2Message} instances to a grib file.
1145
1146 @param filename: file to write grib data to.
1147 @param grbs: a list of L{Grib2Message} instances.
1148 """
1149 gribfile = open(filename,'wb')
1150 for grb in grbs:
1151 try:
1152 f = open(grb._grib_filename,'rb')
1153 except TypeError:
1154 f = StringIO(grb._grib_filename)
1155 f.seek(grb._grib_message_byteoffset)
1156 gribmsg = f.read(grb._grib_message_length)
1157 f.close()
1158 gribfile.write(gribmsg)
1159 gribfile.close()
1160
1161
1162
1164 """return yyyy,mm,dd,min,ss from section 1"""
1165 yyyy=idsect[5]
1166 mm=idsect[6]
1167 dd=idsect[7]
1168 hh=idsect[8]
1169 min=idsect[9]
1170 ss=idsect[10]
1171 return yyyy,mm,dd,hh,min,ss
1172
1174 """unpack section 1 given starting point in bytes
1175 used to test pyrex interface to g2_unpack1"""
1176 idsect = []
1177 pos = pos + 5
1178 idsect.append(struct.unpack('>h',gribmsg[pos:pos+2])[0])
1179 pos = pos + 2
1180 idsect.append(struct.unpack('>h',gribmsg[pos:pos+2])[0])
1181 pos = pos + 2
1182 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1183 pos = pos + 1
1184 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1185 pos = pos + 1
1186 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1187 pos = pos + 1
1188 idsect.append(struct.unpack('>h',gribmsg[pos:pos+2])[0])
1189 pos = pos + 2
1190 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1191 pos = pos + 1
1192 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1193 pos = pos + 1
1194 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1195 pos = pos + 1
1196 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1197 pos = pos + 1
1198 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1199 pos = pos + 1
1200 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1201 pos = pos + 1
1202 idsect.append(struct.unpack('>B',gribmsg[pos:pos+1])[0])
1203 pos = pos + 1
1204 return np.array(idsect,'i'),pos
1205
1207 """repeat last item in listin, until len(listin) = numfields"""
1208 if len(listin) < numfields:
1209 last = listin[-1]
1210 for n in range(len(listin),numfields):
1211 listin.append(last)
1212 return listin
1213
1215 try:
1216 flist = functools.reduce(operator.add,lst)
1217 except NameError:
1218 import functools
1219 flist = functools.reduce(operator.add,lst)
1220 return flist
1221
1222
1224 """
1225 Class for encoding data into a GRIB2 message.
1226 - Creating a class instance (L{__init__}) initializes the message and adds
1227 sections 0 and 1 (the indicator and identification sections),
1228 - method L{addgrid} adds a grid definition (section 3) to the messsage.
1229 - method L{addfield} adds sections 4-7 to the message (the product
1230 definition, data representation, bitmap and data sections).
1231 - method L{end} adds the end section (section 8) and terminates the message.
1232
1233
1234 A GRIB Edition 2 message is a machine independent format for storing
1235 one or more gridded data fields. Each GRIB2 message consists of the
1236 following sections:
1237 - SECTION 0: Indicator Section - only one per message
1238 - SECTION 1: Identification Section - only one per message
1239 - SECTION 2: (Local Use Section) - optional
1240 - SECTION 3: Grid Definition Section
1241 - SECTION 4: Product Definition Section
1242 - SECTION 5: Data Representation Section
1243 - SECTION 6: Bit-map Section
1244 - SECTION 7: Data Section
1245 - SECTION 8: End Section
1246
1247 Sequences of GRIB sections 2 to 7, 3 to 7, or sections 4 to 7 may be repeated
1248 within a single GRIB message. All sections within such repeated sequences
1249 must be present and shall appear in the numerical order noted above.
1250 Unrepeated sections remain in effect until redefined.
1251
1252 Note: Writing section 2 (the 'local use section') is
1253 not yet supported.
1254
1255 @ivar msg: A binary string containing the GRIB2 message.
1256 After the message has been terminated by calling
1257 the L{end} method, this string can be written to a file.
1258 """
1259
1260 - def __init__(self, discipline, idsect):
1261 """
1262 create a Grib2Enecode class instance given the GRIB2 discipline
1263 parameter and the identification section (sections 0 and 1).
1264
1265 The GRIB2 message is stored as a binary string in instance variable L{msg}.
1266
1267 L{addgrid}, L{addfield} and L{end} class methods must be called to complete
1268 the GRIB2 message.
1269
1270 @param discipline: Discipline or GRIB Master Table Number (Code Table 0.0).
1271 (0 for meteorlogical, 1 for hydrological, 2 for land surface, 3 for space,
1272 10 for oceanographic products).
1273
1274 @param idsect: Sequence containing identification section (section 1).
1275 - idsect[0]=Id of orginating centre (Common Code
1276 U{Table C-1<http://www.nws.noaa.gov/tg/GRIB_C1.htm>})
1277 - idsect[1]=Id of orginating sub-centre (local table)
1278 - idsect[2]=GRIB Master Tables Version Number (Code
1279 U{Table 1.0
1280 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-0.shtml>})
1281 - idsect[3]=GRIB Local Tables Version Number (Code
1282 U{Table 1.1
1283 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-1.shtml>})
1284 - idsect[4]=Significance of Reference Time (Code
1285 U{Table 1.2
1286 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-2.shtml>})
1287 - idsect[5]=Reference Time - Year (4 digits)
1288 - idsect[6]=Reference Time - Month
1289 - idsect[7]=Reference Time - Day
1290 - idsect[8]=Reference Time - Hour
1291 - idsect[9]=Reference Time - Minute
1292 - idsect[10]=Reference Time - Second
1293 - idsect[11]=Production status of data (Code
1294 U{Table 1.3
1295 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-3.shtml>})
1296 - idsect[12]=Type of processed data (Code
1297 U{Table
1298 1.4<http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table1-4.shtml>})
1299 """
1300 self.msg,msglen=g2clib.grib2_create(np.array([discipline,2],np.int32),np.array(idsect,np.int32))
1301
1302 - def addgrid(self,gdsinfo,gdtmpl,deflist=None):
1303 """
1304 Add a grid definition section (section 3) to the GRIB2 message.
1305
1306 @param gdsinfo: Sequence containing information needed for the grid definition section.
1307 - gdsinfo[0] = Source of grid definition (see Code
1308 U{Table 3.0
1309 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-0.shtml>})
1310 - gdsinfo[1] = Number of grid points in the defined grid.
1311 - gdsinfo[2] = Number of octets needed for each additional grid points defn.
1312 Used to define number of points in each row for non-reg grids (=0 for
1313 regular grid).
1314 - gdsinfo[3] = Interp. of list for optional points defn (Code
1315 U{Table 3.11
1316 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-11.shtml>})
1317 - gdsinfo[4] = Grid Definition Template Number (Code
1318 U{Table 3.1
1319 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table3-1.shtml>})
1320
1321 @param gdtmpl: Contains the data values for the specified Grid Definition
1322 Template ( NN=gdsinfo[4] ). Each element of this integer
1323 array contains an entry (in the order specified) of Grid
1324 Definition Template 3.NN
1325
1326 @param deflist: (Used if gdsinfo[2] != 0) Sequence containing the
1327 number of grid points contained in each row (or column)
1328 of a non-regular grid.
1329 """
1330 if deflist is not None:
1331 dflist = np.array(deflist,'i')
1332 else:
1333 dflist = None
1334 self.scanmodeflags = None
1335 gdtnum = gdsinfo[4]
1336 if gdtnum in [0,1,2,3,40,41,42,43,44,203,205,32768,32769]:
1337 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4]
1338 elif gdtnum == 10:
1339 self.scanmodeflags = _dec2bin(gdtmpl[15])[0:4]
1340 elif gdtnum == 20:
1341 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
1342 elif gdtnum == 30:
1343 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
1344 elif gdtnum == 31:
1345 self.scanmodeflags = _dec2bin(gdtmpl[17])[0:4]
1346 elif gdtnum == 90:
1347 self.scanmodeflags = _dec2bin(gdtmpl[16])[0:4]
1348 elif gdtnum == 110:
1349 self.scanmodeflags = _dec2bin(gdtmpl[15])[0:4]
1350 elif gdtnum == 120:
1351 self.scanmodeflags = _dec2bin(gdtmpl[6])[0:4]
1352 elif gdtnum == 204:
1353 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4]
1354 elif gdtnum in [1000,1100]:
1355 self.scanmodeflags = _dec2bin(gdtmpl[12])[0:4]
1356 self.msg,msglen=g2clib.grib2_addgrid(self.msg,np.array(gdsinfo,'i'),np.array(gdtmpl,'i'),dflist)
1357
1358 - def addfield(self,pdtnum,pdtmpl,drtnum,drtmpl,field,coordlist=None):
1359 """
1360 Add a product definition section, data representation section,
1361 bitmap section and data section to the GRIB2 message (sections 4-7).
1362 Must be called after grid definition section is created with L{addgrid}.
1363
1364 @param pdtnum: Product Definition Template Number (see Code U{Table
1365 4.0<http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table4-0.shtml>})
1366
1367 @param pdtmpl: Sequence with the data values for the specified Product Definition
1368 Template (N=pdtnum). Each element of this integer
1369 array contains an entry (in the order specified) of Product
1370 Definition Template 4.N
1371
1372 @param drtnum: Data Representation Template Number (see Code
1373 U{Table 5.0
1374 <http://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_table5-0.shtml>})
1375
1376 @param drtmpl: Sequence with the data values for the specified Data Representation
1377 Template (N=drtnum). Each element of this integer
1378 array contains an entry (in the order specified) of Data
1379 Representation Template 5.N
1380 Note that some values in this template (eg. reference
1381 values, number of bits, etc...) may be changed by the
1382 data packing algorithms.
1383 Use this to specify scaling factors and order of
1384 spatial differencing, if desired.
1385
1386 @param field: numpy array of data points to pack.
1387 If field is a masked array, then a bitmap is created from
1388 the mask.
1389
1390 @param coordlist: Sequence containing floating point values intended to document
1391 the vertical discretization with model data
1392 on hybrid coordinate vertical levels. Default None.
1393 """
1394 if not hasattr(self,'scanmodeflags'):
1395 raise ValueError('addgrid must be called before addfield')
1396
1397
1398 if self.scanmodeflags is not None:
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409 if self.scanmodeflags[3]:
1410 fieldsave = field.astype('f')
1411 field[1::2,:] = fieldsave[1::2,::-1]
1412 fld = field.astype('f')
1413 if ma.isMA(field):
1414 bmap = 1 - np.ravel(field.mask.astype('i'))
1415 bitmapflag = 0
1416 else:
1417 bitmapflag = 255
1418 bmap = None
1419 if coordlist is not None:
1420 crdlist = np.array(coordlist,'f')
1421 else:
1422 crdlist = None
1423 self.msg,msglen=g2clib.grib2_addfield(self.msg,pdtnum,np.array(pdtmpl,'i'),crdlist,drtnum,np.array(drtmpl,'i'),np.ravel(fld),bitmapflag,bmap)
1424
1426 """
1427 Add an end section (section 8) to the GRIB2 message.
1428 A GRIB2 message is not complete without an end section.
1429 Once an end section is added, the GRIB2 message can be
1430 output to a file.
1431 """
1432 self.msg,msglen=g2clib.grib2_end(self.msg)
1433