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