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