Module ncepgrib2

# Source Code for Module ncepgrib2

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