Module ncepgrib2

Source Code for Module ncepgrib2

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