Package ncepgrib2

Source Code for Package ncepgrib2

   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  # Code Table 3.2: Shape of the Earth. 
  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   
295 -def _dec2bin(val, maxbits = 8):
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
306 -def _putieeeint(r):
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
313 -def _getieeeint(i):
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
320 -def _isString(string):
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
326 -class Grib2Message:
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 """
419 - def __init__(self,**kwargs):
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 # grid information 428 gdsinfo = self.grid_definition_info 429 gdtnum = self.grid_definition_template_number 430 gdtmpl = self.grid_definition_template 431 reggrid = gdsinfo[2] == 0 # gdsinfo[2]=0 means regular 2-d grid 432 # shape of the earth. 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. # convert to m from km 454 scaledearthRminor = gdtmpl[6] 455 earthRminorscale = gdtmpl[5] 456 self.earthRminor = math.pow(10,-earthRminorscale)*scaledearthRminor 457 self.earthRminor = self.earthRminor*1000. # convert to m from km 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: # 'reduced' gaussian grid. 475 self.points_in_y_direction = gdtmpl[8] 476 if gdtnum in [0,1,203,205,32768]: # regular or rotated lat/lon grid 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: # mercator 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: # stereographic 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: # lambert conformal 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: # albers equal area. 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: # gaussian grid. 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: # spectral coefficients. 565 self.spectral_truncation_parameters = (gdtmpl[0],gdtmpl[1],gdtmpl[2]) 566 self.scanmodeflags = [None,None,None,None] # doesn't apply 567 elif gdtnum == 90: # near-sided vertical perspective satellite projection 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 # if lat_0 is equator, it's a geostationary view. 574 if self.proj4_lat_0 == 0.: # if lat_0 is equator, it's a 575 self.proj4_proj = 'geos' 576 # general case of 'near-side perspective projection' (untested) 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 # latitude of horizon on central meridian 584 lonmax = 90.-(180./N.pi)*N.arcsin(self.earthRmajor/self.proj4_h) 585 # longitude of horizon on equator 586 latmax = 90.-(180./N.pi)*N.arcsin(self.earthRminor/self.proj4_h) 587 # truncate to nearest thousandth of a degree (to make sure 588 # they aren't slightly over the horizon) 589 latmax = int(1000*latmax)/1000. 590 lonmax = int(1000*lonmax)/1000. 591 # h is measured from surface of earth at equator. 592 self.proj4_h = self.proj4_h - self.earthRmajor 593 # width and height of visible projection 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: # azimuthal equidistant. 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: # curvilinear orthogonal 610 self.scanmodeflags = _dec2bin(gdtmpl[18])[0:4] 611 # missing value. 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
619 - def __repr__(self):
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 # make sure scan mode is supported. 646 # if there is no 'scanmodeflags', then grid is not supported. 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 # default order=0 is missing values or bitmap exists. 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 # apply bitmap. 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 # missing values instead of bitmap 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: # rectangular grid. 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: # ECMWF 'reduced' global gaussian grid. 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 # check scan modes for rect grids. 716 if nx is not None and ny is not None: 717 # rows scan in the -x direction (so flip) 718 if self.scanmodeflags[0]: 719 fldsave = fld.astype('f') # casting makes a copy 720 fld[:,:] = fldsave[:,::-1] 721 # columns scan in the -y direction (so flip) 722 if not self.scanmodeflags[1]: 723 fldsave = fld.astype('f') # casting makes a copy 724 fld[:,:] = fldsave[::-1,:] 725 # adjacent rows scan in opposite direction. 726 # (flip every other row) 727 if self.scanmodeflags[3]: 728 fldsave = fld.astype('f') # casting makes a copy 729 fld[1::2,:] = fldsave[1::2,::-1] 730 return fld
731
732 - def grid(self):
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 # gdsinfo[2]=0 means regular 2-d grid 747 projparams = {} 748 projparams['a']=self.earthRmajor 749 projparams['b']=self.earthRminor 750 if gdtnum == 0: # regular lat/lon grid 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 # flip if scan mode says to. 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) # make 2-d arrays. 763 elif gdtnum == 40: # gaussian grid (only works for global!) 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: # ECMWF 'reduced' gaussian grid. 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 # compute gaussian lats (north to south) 775 lats = gaulats(nlats) 776 if lat1 < lat2: # reverse them if necessary 777 lats = lats[::-1] 778 # flip if scan mode says to. 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) # make 2-d arrays 784 # mercator, lambert conformal, stereographic, albers equal area, azimuthal equidistant 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: # mercator. 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: # stereographic 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]: # lambert, albers 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: # azimuthal equidistant 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: # satellite projection. 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 # set lons,lats to 1.e30 where undefined 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
857 -def Grib2Decode(filename,gribmsg=False):
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 # loop over grib messages, read section 0, get entire grib message. 880 disciplines = [] 881 startingpos = [] 882 msglen = [] 883 while 1: 884 # find next occurence of string 'GRIB' (or EOF). 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 # at EOF 892 # otherwise, start (='GRIB') contains indicator message (section 0) 893 startpos = f.tell()-4 894 f.seek(2,1) # next two octets are reserved 895 # get discipline info. 896 disciplines.append(struct.unpack('>B',f.read(1))[0]) 897 # check to see it's a grib edition 2 file. 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 # read in entire grib message. 905 f.seek(startpos) 906 gribmsg = f.read(lengrib) 907 # make sure the message ends with '7777' 908 end = gribmsg[-4:lengrib].decode('ascii','ignore') 909 if end != '7777': 910 raise IOError('partial GRIB message (no "7777" at end)') 911 # do next message. 912 nmsg=nmsg+1 913 # if no grib messages found, nmsg is still 0 and it's not GRIB. 914 if nmsg==0: 915 raise IOError('not a GRIB file') 916 # now for each grib message, find number of fields. 917 numfields = [] 918 f.seek(0) # rewind file. 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 #if sectnum == 2: numlocal=numlocal+1 935 pos = pos + lensect 936 #print sectnum,lensect,pos 937 #print n+1,len(gribmsg),numfields,numlocal 938 numfields.append(numflds) 939 # decode each section in grib message (sections 1 and above). 940 gdtnum = [] # grid defn template number from sxn 3 941 gdtmpl = [] # grid defn template from sxn 3 942 gdeflist = [] # optional grid definition list from sxn 3 943 gdsinfo = [] # grid definition section info from sxn3 944 pdtmpl = [] # product defn template from sxn 4 945 pdtnum = [] # product defn template number from sxn 4 946 coordlist = [] # vertical coordinate info from sxn 4 947 drtmpl = [] # data representation template from sxn 5 948 drtnum = [] # data representation template number from sxn 5 949 ndpts = [] # number of data points to be unpacked (from sxn 5) 950 bitmapflag = [] # bit-map indicator flag from sxn 6 951 bitmap = [] # bitmap from sxn 6. 952 pos7 = [] # byte offset for section 7. 953 localsxn = [] # local use sections. 954 msgstart = [] # byte offset in file for message start. 955 msglength = [] # length of the message in bytes. 956 message = [] # the actual grib message. 957 identsect = [] # identification section (section 1). 958 discipline = [] # discipline code. 959 for n in range(nmsg): 960 spos = startingpos[n] 961 lengrib = msglen[n] 962 #gribmsg = gribmsgs[n] 963 f.seek(spos) 964 gribmsg = f.read(lengrib) 965 discipl = disciplines[n] 966 lensect0 = 16 967 # get length of section 1 and section number. 968 #lensect1 = struct.unpack('>i',gribmsg[lensect0:lensect0+4])[0] 969 #sectnum1 = struct.unpack('>B',gribmsg[lensect0+4])[0] 970 #print 'sectnum1, lensect1 = ',sectnum1,lensect1 971 # unpack section 1, octets 1-21 (13 parameters). This section 972 # can occur only once per grib message. 973 #idsect,pos = _unpack1(gribmsg,lensect0) # python version 974 idsect,pos = g2clib.unpack1(gribmsg,lensect0,N.empty) # c version 975 # loop over rest of sections in message. 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 # check to see if this is the end of the message. 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 # section 2, local use section. 996 if sectnum == 2: 997 # "local use section", used by NDFD to store WX 998 # strings. This section is returned as a raw 999 # bytestring for further dataset-specific parsing, 1000 # not as a numpy array. 1001 localsxns.append(gribmsg[pos+5:pos+lensect]) 1002 pos = pos + lensect 1003 # section 3, grid definition section. 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 # section, product definition section. 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 # section 5, data representation section. 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 # section 6, bit-map section. 1023 elif sectnum == 6: 1024 bmap,bmapflag = g2clib.unpack6(gribmsg,gds[1],pos,N.empty) 1025 #bitmapflag = struct.unpack('>B',gribmsg[pos+5])[0] 1026 if bmapflag == 0: 1027 bitmaps.append(bmap.astype('b')) 1028 # use last defined bitmap. 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 # section 7, data section (nothing done here, 1038 # data unpacked when getfld method is invoked). 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 # extend by repeating last value for all remaining fields. 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 # private methods and functions below here. 1147
1148 -def _getdate(idsect):
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
1158 -def _unpack1(gribmsg,pos):
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
1191 -def _repeatlast(numfields,listin):
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
1199 -def _flatten(numfields, listin):
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
1206 -class Grib2Encode:
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
1374 - def end(self):
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