Module ncepgrib2

Source Code for Module ncepgrib2

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