Module ncepgrib2

Source Code for Module ncepgrib2

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