def __init__(self, center=None, radius=None, forecast_date=None, gfs_data_file=None, hrrr_data_file=None, keep_files=False):
        """
        Construct a wind model centered at a given Latitude/Longitude with a radius defined in miles
        :param center:
        :param radius:
        :param forecast_date:
        """
        if gfs_data_file is None and hrrr_data_file is None and (center is None or radius is None or forecast_date is None):
            raise Exception("Invalid center or radius for wind model")

        self.forecast_date = forecast_date
        self.gfs_data_file = gfs_data_file
        self.gfs_height_map = None
        self.gfs_latlons = None
        self.hrrr_data_file = hrrr_data_file
        self.hrrr_height_map = None
        self.hrrr_latlons = None

        if gfs_data_file is not None:
            self.gfs_height_map, self.gfs_latlons = self._parse_grbs(pygrib.open(gfs_data_file))
        if hrrr_data_file is not None:
            self.hrrr_height_map, self.hrrr_latlons = self._parse_grbs(pygrib.open(hrrr_data_file))

        if gfs_data_file is None and hrrr_data_file is None:
            lat_radius = change_in_latitude_miles(radius)
            lon_radius = max(change_in_longitude_miles(center[0] - lat_radius, radius), change_in_longitude_miles(center[0] + lat_radius, radius))

            self.NW_bound = [center[0] + lat_radius, center[1] - lon_radius]
            self.SE_bound = [center[0] - lat_radius, center[1] + lon_radius]

            self._load_gfs(keep_files)

	def __init__(self,source,target):
		
		def lon_lat_to_cartesian(lon,lat,R=1):
			"""
			calculates lon, lat coordinates of a point on a sphere with
			radius R
			"""
			
			lon_r=np.radians(lon)
			lat_r=np.radians(lat)

			x=R*np.cos(lat_r)*np.cos(lon_r)
			y=R*np.cos(lat_r)*np.sin(lon_r)
			z=R*np.sin(lat_r)
			return x,y,z
	
		"""
		read data from files
		"""
		self.source=source
		self.target=target

		self.D1=pygrib.open(self.source)
		self.latSOURCE,self.lonSOURCE=self.D1[8].latlons()
				
		self.D2=pygrib.open(self.target)
		self.latTARGET,self.lonTARGET=self.D2[1].latlons()	

		"""
		flatten grid coordinates into cartesian x,y,z
		"""
		self.xs,self.ys,self.zs=lon_lat_to_cartesian(self.lonSOURCE.flatten(),self.latSOURCE.flatten())
		self.xt,self.yt,self.zt=lon_lat_to_cartesian(self.lonTARGET.flatten(),self.latTARGET.flatten())

	def plot_algorithms(file1,file2):	
		file1=pygrib.open(file1)
		temp_var_file1=file1[8].values

		file2=pygrib.open(file2)	
		temp_var_file2=file2[3].values
		
		file1.close()
		file2.close()
		
		plt.figure(figsize=(10,5))
		plt.subplot(121)
		plt.pcolormesh(temp_var_file1)
		plt.xlim([0, temp_var_file1.shape[0]])
		plt.ylim([0, temp_var_file1.shape[1]])
		plt.colorbar()
		plt.title("from file1")

		plt.subplot(122)
		plt.pcolormesh(temp_var_file2)
		plt.colorbar()
		plt.xlim([0, temp_var_file2.shape[0]])
		plt.ylim([0, temp_var_file2.shape[1]])
		plt.title("from file2");
	
		plt.show()

def pcpn(grids, valid, iarchive):
    """Attempt to use MRMS or stage IV pcpn here"""
    floor = datetime.datetime(2014, 11, 1)
    floor = floor.replace(tzinfo=pytz.timezone("UTC"))
    if valid < floor:
        # Use stageIV
        ts = (valid + datetime.timedelta(minutes=60)).replace(minute=0)
        gribfn = ts.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/stage4/ST4."
                              "%Y%m%d%H.01h.grib"))
        if not os.path.isfile(gribfn):
            return
        grbs = pygrib.open(gribfn)
        grib = grbs[1]
        lats, lons = grib.latlons()
        vals = grib.values
        nn = NearestNDInterpolator((lons.flatten(), lats.flatten()),
                                   vals.flatten())
        grids['pcpn'] = nn(XI, YI)
        return
    fn = None
    i = 0
    while i < 10:
        ts = valid - datetime.timedelta(minutes=i)
        if ts.minute % 2 == 0:
            testfn = mrms_util.fetch('PrecipRate', ts)
            if testfn is not None:
                fn = testfn
                break
        i += 1
    if fn is None:
        print("Warning, no PrecipRate data found!")
        return
    fp = gzip.GzipFile(fn, 'rb')
    (_, tmpfn) = tempfile.mkstemp()
    tmpfp = open(tmpfn, 'wb')
    tmpfp.write(fp.read())
    tmpfp.close()
    grbs = pygrib.open(tmpfn)
    values = grbs[1]['values']
    # just set -3 (no coverage) to 0 for now
    values = np.where(values < 0, 0, values)
    map(os.unlink, [fn, tmpfn])

    # 3500, 7000, starts in upper left
    top = int((55. - reference.IA_NORTH) * 100.)
    bottom = int((55. - reference.IA_SOUTH) * 100.)

    right = int((reference.IA_EAST - -130.) * 100.) - 1
    left = int((reference.IA_WEST - -130.) * 100.)

    # two minute accumulation is in mm/hr / 60 * 5
    # stage IV is mm/hr
    grids['pcpn'] = np.flipud(values[top:bottom, left:right]) / 12.0

def compute(valid):
    ''' Get me files '''
    prob = None
    for hr in range(-15,0):
        ts = valid + datetime.timedelta(hours=hr)
        fn = ts.strftime("hrrr.ref.%Y%m%d%H00.grib2")
        if not os.path.isfile(fn):
            continue

        grbs = pygrib.open(fn)
        gs = grbs.select(level=1000,forecastTime=(-1 * hr * 60))
        ref = generic_filter(gs[0]['values'], np.max, size=10)
        if prob is None:
            lats, lons = gs[0].latlons()
            prob = np.zeros( np.shape(ref) )
        
        prob = np.where(ref > 29, prob+1, prob)

    prob = np.ma.array(prob / 15. * 100.)
    prob.mask = np.ma.where(prob < 1, True, False)    
    
    m = MapPlot(sector='iowa',
                title='HRRR Composite Forecast 4 PM 20 May 2014 30+ dbZ Reflectivity',
                subtitle='frequency of previous 15 model runs all valid at %s, ~15km smoothed' % (valid.astimezone(pytz.timezone("America/Chicago")).strftime("%-d %b %Y %I:%M %p %Z"),))

    m.pcolormesh(lons, lats, prob, np.arange(0,101,10), units='%',
                     clip_on=False)
    m.map.drawcounties()
    m.postprocess(filename='test.ps')
    m.close()

def read_grib(filename, dataset, dt):
    '''
    read dataset in grib file filename
    with date and hour specified by datetime object dt
    '''
    localFilename = copyToLocal(filename)

    hour = dt.hour
    date = int(dt.strftime('%Y%m%d'))

    try:
        g = pygrib.open(localFilename)
    except IOError:
        raise Exception('Could not open grib file %s' % (localFilename))

    data = None
    grblist = []
    for grb in g:
        grblist.append(grb)
        if (grb.name == dataset) and (grb.hour == hour) and (grb.dataDate == date):
            data = grb.values

    g.close()

    if data == None:
        print 'List of datasets in grib file:'
        for grb in grblist:
            print grb
            # for j in sorted(grb.keys()): print '->', j
            # print '===', '"%s"' % (grb.name), grb.dataDate, type(grb.dataDate), grb.hour

        raise Exception('Could not read "%s" at %s in file %s' % (dataset, dt, localFilename))

    return data

def gettigge(center,idate,var,lev,ftyp='cntl'):
    YYYYMM=idate.strftime('%Y%m')
    YYYYMMDDHH=idate.strftime('%Y%m%d%H')
    varname = var+str(lev)
    fname=FNAMEBASE.format(center=center.lower(),
                           ftyp=ftyp,
                           varname=varname,
                           YYYYMM=YYYYMM,
                           YYYYMMDDHH=YYYYMMDDHH)
    print fname
    #
    shortname=var
    if var=='z': shortname='gh'
    grbs = pygrib.open(fname)
    grb = grbs.select(shortName=shortname,level=lev)
    xn = grb[0]['Ni']
    yn = grb[0]['Nj']
    tn = len(grb)
    lat,lon = grb[0].latlons()
    out = np.empty((tn,yn,xn))
    tyme = np.empty(tn,dtype=np.object)
    for t, gr in enumerate(grb):
        out[t,...] = gr.values
        tyme[t] = gr.validDate
    grbs.close()
    dims = ['tyme','lat','lon']
    grid = McGrid(out, lon=lon[0,:], lat=lat[::-1,0], lev=lev, tyme=tyme, dims=dims)
    out = McField(out,name=var,grid=grid)
    return out

def do():
    now = sts
    while now < ets:
        fn = now.strftime("/mesonet/ARCHIVE/data/%Y/%m/%d/model/hrrr/%H/hrrr.t%Hz.3kmf00.grib2")
        if not os.path.isfile(fn):
            print fn
            now += interval
            continue
        grbs = pygrib.open(fn)
        try:
            gs = grbs.select(name='2 metre temperature')
        except:
            print fn
            now += interval
            continue

        g = gs[0]['values']
        if now == sts:
            lats,lons = gs[0].latlons()
            maxinterval = np.zeros(np.shape(g))
            current = np.zeros(np.shape(g))
        print np.max(g), np.min(g)
        current = np.where(g < 273, current + 1, 0)
        maxinterval = np.where(current > maxinterval, current, maxinterval)
        now += interval

    np.save('maxinterval.npy', np.array(maxinterval))
    np.save('lats.npy', lats)
    np.save('lons.npy', lons)

def doday():
    """
    Create a plot of precipitation stage4 estimates for some day
    """
    sts = mx.DateTime.DateTime(2013,5,25,12)
    ets = mx.DateTime.DateTime(2013,5,31,12)
    interval = mx.DateTime.RelativeDateTime(days=1)
    now = sts
    total = None
    while now < ets:
        fp = "/mesonet/ARCHIVE/data/%s/stage4/ST4.%s.24h.grib" % (
            now.strftime("%Y/%m/%d"), 
            now.strftime("%Y%m%d%H") )
        if os.path.isfile(fp):
            lts = now
            grbs = pygrib.open(fp)

            if total is None:
                g = grbs[1]
                total = g["values"]
                lats, lons = g.latlons()
            else:
                total += grbs[1]["values"]
            grbs.close()
        now += interval
        
    m = MapPlot(sector='iowa', title='NOAA Stage IV & Iowa ASOS Precipitation',
                subtitle='25-30 May 2013')
    m.pcolormesh(lons, lats, total / 25.4, numpy.arange(0,14.1,1), latlon=True,
                 units='inch')
    m.drawcounties()
    m.plot_values(dlons, dlats, dvals, '%.02f')
    m.postprocess(filename='test.svg')
    import iemplot
    iemplot.makefeature('test')

def readfile(infile):
    '''Read the GRIB file and return the regional (SCS) sst data'''
    grbs=pygrib.open(infile)
    grb=grbs.select(name='Temperature')[0]
    sst, lats, lons = grb.data(lat1=14,lat2=30,lon1=106,lon2=126)
    grbs.close()
    return sst

def ncepgfs(lonStart=3, lonEnd=4, latStart=41, latEnd=42, \
            m=None, name='ncepgfs', contour=None):
    """
    Plot latest NCEP GFS field from \
    http://nomad1.ncep.noaa.gov/pub/gfs_master/
    """

    base = '/media/SOLabNFS/SERVERS/media/hyrax/data/auxdata/model/ncep/gfs/'
    outdir = '/home/mag/'
    
    fn = 'gfs20101218/gfs.t18z.master.grbf00'
#    fn = 'gfs20101218/gfs.t18z.master.grbf03'
    
    refDate = datetime.datetime(1978,1,1,0,0,0)
    wantedDate = datetime.datetime(2010,12,18,17,39,0)

    # Get the file date
    grbs = pygrib.open(base + fn)
    grb = grbs.message(1)
    fileDate = grb.analDate
    filetime = fileDate.strftime("%Y%m%d_%H%M")
    print "File time: ", filetime
    print "File date: ", fileDate
    
    lats, lons = grb.latlons()
    u10 = grbs.message(9)['values']
    v10 = grbs.message(9)['values']
#    u = grbs.message(3)['values']
#    v = grbs.message(4)['values']
    w = sqrt(u10**2+v10**2)

def run(ts, routes):
    """ Run for a given UTC timestamp """
    fn = ts.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/model/rtma/%H/"
                      "rtma.t%Hz.awp2p5f000.grib2"))
    if not os.path.isfile(fn):
        print 'wind_power.py missing', fn
        return

    grb = pygrib.open(fn)
    try:
        u = grb.select(name='10 metre U wind component')[0]
        v = grb.select(name='10 metre V wind component')[0]
    except:
        print('Missing u/v wind for wind_power.py\nFN: %s' % (fn,))
        return
    mag = (u['values']**2 + v['values']**2)**.5

    mag = (mag * 1.35)**3 * 0.002641
    # 0.002641

    lats, lons = u.latlons()
    lts = ts.astimezone(pytz.timezone("America/Chicago"))
    pqstr = ("plot %s %s00 midwest/rtma_wind_power.png "
             "midwest/rtma_wind_power_%s00.png png"
             ) % (routes, ts.strftime("%Y%m%d%H"), ts.strftime("%H"))
    m = MapPlot(sector='midwest',
                title=(r'Wind Power Potential :: '
                       '(speed_mps_10m * 1.35)$^3$ * 0.002641'),
                subtitle=('valid: %s based on NOAA Realtime '
                          'Mesoscale Analysis'
                          ) % (lts.strftime("%d %b %Y %I %p")))
    m.pcolormesh(lons, lats, mag, numpy.array(levels), units='MW')

    m.postprocess(pqstr=pqstr)

def create_sub(full_grib, data_dir='.'):
  ''' Extracts a subset of GRIB records from full_grib and 
  saves them in a file named the same as full_grib, but in 
  the data_dir folder.
  '''
  #names=['Temperature', 'Soil temperature', 'Total precipitation', 'Wind speed (gust)', 'Water equivalent of accumulated snow depth']
  name, ext = os.path.splitext(full_grib)
  out_file = '{2}/{0}_sub{1}'.format(name, ext, data_dir)
  grbout = open(out_file, 'wb')

  grbs = pygrib.open(full_grib)
  surface_names=['t', 'acpcp', 'gust', 'sdwe']
  for name in surface_names:
    print 'finding {0}'.format(name)
    grbs.seek(0)
    for grb in grbs.select(shortName=name, typeOfLevel='surface'):
      grbout.write(grb.tostring())
  agl_names=['10u', '10v']
  grbs.seek(0)
  for name in agl_names:
    print 'finding {0}'.format(name)
    for grb in grbs.select(shortName=name, typeOfLevel='heightAboveGround'):
      grbout.write(grb.tostring())
  grbout.close()
  return out_file

def containsMessages(filename):
    grbs = pg.open(filename)
    count = grbs.messages
    #print "messages: %s" % (count)
    if count:
        return grbs
    return None

def run(ts, routes):
    """ Run for a given UTC timestamp """
    fn = ts.strftime("/mesonet/ARCHIVE/data/%Y/%m/%d/model/rtma/%H/rtma.t%Hz.awp2p5f000.grib2")
    if not os.path.isfile(fn):
        print "wind_power.py missing", fn
        return

    grb = pygrib.open(fn)
    u = grb.select(name="10 metre U wind component")[0]
    v = grb.select(name="10 metre V wind component")[0]
    mag = (u["values"] ** 2 + v["values"] ** 2) ** 0.5

    mag = (mag * 1.35) ** 3 * 0.002641
    # 0.002641

    lats, lons = u.latlons()
    lts = ts.astimezone(pytz.timezone("America/Chicago"))
    pqstr = "plot %s %s00 midwest/rtma_wind_power.png midwest/rtma_wind_power_%s00.png png" % (
        routes,
        ts.strftime("%Y%m%d%H"),
        ts.strftime("%H"),
    )
    m = MapPlot(
        sector="midwest",
        title=r"Wind Power Potential :: (speed_mps_10m * 1.35)$^3$ * 0.002641",
        subtitle="valid: %s based on NOAA Realtime Mesoscale Analysis" % (lts.strftime("%d %b %Y %I %p")),
    )
    m.pcolormesh(lons, lats, mag, numpy.array(levels), units="MW")

    m.postprocess(pqstr=pqstr)

 def load_data(self):
     for m, member in enumerate(self.members):
         for f, forecast_date in enumerate(self.forecast_dates.to_pydatetime()):
             dt = int((forecast_date - self.run_date).total_seconds() / 3600)
             filename_args = (self.ensemble_name, member, self.ml_model, self.variable,
                              forecast_date.strftime("%Y%m%d%H%M"))
             filename = self.path + self.run_date.strftime("%Y%m%d") + \
                        "/{0}_{1}_{2}_{3}_{4}.grib2".format(*filename_args)
             if not exists(filename):
                 filename_args = (self.ensemble_name, member, self.ml_model, self.variable,
                                  self.run_date.strftime("%Y%m%d%H") + "f{0:02d}".format(dt))
                 filename = self.path + self.run_date.strftime("%Y%m%d") + \
                            "/{0}_{1}_{2}_{3}_{4}.grib2".format(*filename_args)
                 if not exists(filename):
                     continue
             grbs = pygrib.open(filename)
             if self.lon is None:
                 self.lat, self.lon = grbs[self.message_number].latlons()
                 self.projparams = grbs[self.message_number].projparams
                 self.proj = Proj(grbs[self.message_number].projparams)
                 self.x, self.y = self.proj(self.lon, self.lat)
                 self.x /= 1000.0
                 self.y /= 1000.0
                 self.dx = grbs[self.message_number]['DxInMetres'] / 1000.0
                 self.i, self.j = np.indices(self.lon.shape)
             data = grbs[self.message_number].values
             data *= 1000.0
             if self.data is None:
                 self.data = np.empty((len(self.members), len(self.forecast_dates),
                                       data.shape[0], data.shape[1]), dtype=float)
             self.data[m, f] = data.filled(0)
             grbs.close()
     return

def main():
    """Go!"""
    title = 'NOAA MRMS Q3: RADAR + Guage Corrected Rainfall Estimates + NWS Storm Reports'
    mp = MapPlot(sector='custom',
                 north=42.3, east=-93.0, south=41.65, west=-94.1,
                 axisbg='white',
                 titlefontsize=14,
                 title=title,
                 subtitle='Valid: 14 June 2018')

    shp = shapefile.Reader('cities.shp')
    for record in shp.shapeRecords():
        geo = shape(record.shape)
        mp.ax.add_geometries([geo], ccrs.PlateCarree(), zorder=Z_OVERLAY2,
                             facecolor='None', edgecolor='k', lw=2)

    grbs = pygrib.open('MRMS_GaugeCorr_QPE_24H_00.00_20180614-200000.grib2')
    grb = grbs.message(1)
    pcpn = distance(grb['values'], 'MM').value('IN')
    lats, lons = grb.latlons()
    lons -= 360.
    clevs = [0.01, 0.1, 0.3, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10]
    cmap = nwsprecip()
    cmap.set_over('k')

    mp.pcolormesh(lons, lats, pcpn, clevs, cmap=cmap, latlon=True,
                  units='inch')
    lons, lats, vals, labels = get_data()
    mp.drawcounties()
    mp.plot_values(lons, lats, vals, "%s", labels=labels,
                   labelbuffer=1, labelcolor='white')

    mp.drawcities(labelbuffer=5, minarea=0.2)
    mp.postprocess(filename='test.png')

    def get(self, gmessage, key='values'):
        '''
        Returns the value for the 'key' for a given message number 'gmessage' or
        message field name 'gmessage'.
        '''

        grbs = pygrib.open(self._abspath)

        if type(gmessage) == int:
            mnbr = gmessage
        elif type(gmessage) == str:
            msg_found = False
            msgnum = 1
            while msgnum < self.nmsgs + 1:
                if grbs[msgnum]['parameterName'] == gmessage:
                    msg_found = True
                    break
                msgnum = msgnum + 1

            if msg_found:
                mnbr = msgnum
            else:
                print("No Grib message found with parameter name = %s" %
                      gmessage)
                return None

        if grbs[mnbr].valid_key(key):

            arr = grbs[mnbr][key]
            grbs.close()
            return arr
        else:
            grbs.close()
            return

def GribVals(infile):
    grbs = pygrib.open(infile)    
    for grb in grbs:
          print grb
          
    grb_vals  = pd.DataFrame(grb.values[:])
    return grb_vals

    def __init__(self, filename, filename_info, filetype_info):
        super(GRIBFileHandler, self).__init__(filename, filename_info, filetype_info)

        self._msg_datasets = {}
        self._start_time = None
        self._end_time = None
        try:
            with pygrib.open(filename) as grib_file:
                first_msg = grib_file.message(1)
                last_msg = grib_file.message(grib_file.messages)
                start_time = self._convert_datetime(
                    first_msg, 'validityDate', 'validityTime')
                end_time = self._convert_datetime(
                    last_msg, 'validityDate', 'validityTime')
                self._start_time = start_time
                self._end_time = end_time
                if 'keys' not in filetype_info:
                    self._analyze_messages(grib_file)
                    self._idx = None
                else:
                    self._create_dataset_ids(filetype_info['keys'])
                    self._idx = pygrib.index(filename,
                                             *filetype_info['keys'].keys())
        except (RuntimeError, KeyError):
            raise IOError("Unknown GRIB file format: {}".format(filename))