Python iemre.daily_offset函数代码示例

def do_month(year, month, routes):
    """ Generate a MRMS plot for the month!"""

    sts = datetime.datetime(year,month,1)
    ets = sts + datetime.timedelta(days=35)
    ets = ets.replace(day=1)

    today = datetime.datetime.now()
    if ets > today:
        ets = today

    idx0 = iemre.daily_offset(sts)
    idx1 = iemre.daily_offset(ets)

    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_mrms_daily.nc" % (year,),
                          'r')

    lats = nc.variables['lat'][:]
    lons = nc.variables['lon'][:]
    p01d = np.sum(nc.variables['p01d'][idx0:idx1,:,:],0) / 24.5
    nc.close()

    m = MapPlot(sector='iowa', title='MRMS %s - %s Total Precipitation' % (
            sts.strftime("%-d %b"), 
            (ets - datetime.timedelta(days=1)).strftime("%-d %b %Y")),
            subtitle='Data from NOAA MRMS Project')
    x,y = np.meshgrid(lons, lats)
    bins = [0.01, 0.1, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 7, 8, 12, 16, 20]
    m.pcolormesh(x, y, p01d, bins, units='inches')
    m.drawcounties()
    currentfn = "summary/iowa_mrms_q3_month.png"
    archivefn = sts.strftime("%Y/%m/summary/iowa_mrms_q3_month.png")
    pqstr = "plot %s %s00 %s %s png" % (
                routes, sts.strftime("%Y%m%d%H"), currentfn, archivefn)
    m.postprocess(pqstr=pqstr)

def main():
    """Go Main"""
    total = None
    years = 0.
    for yr in range(1981, 2018):
        print(yr)
        ncfn = "/mesonet/data/prism/%s_daily.nc" % (yr, )
        nc = netCDF4.Dataset(ncfn)
        if total is None:
            lons = nc.variables['lon'][:]
            lats = nc.variables['lat'][:]
            total = np.zeros(nc.variables['tmax'].shape[1:], np.float)
        days = np.zeros(nc.variables['tmax'].shape[1:], np.float)
        sidx = daily_offset(datetime.date(yr, 1, 1))
        eidx = daily_offset(datetime.date(yr, 7, 4))
        for idx in range(sidx, eidx):
            days += np.where(nc.variables['tmax'][idx, :, :] > THRESHOLD,
                             1, 0)
        nc.close()
        years += 1.
        total += days

    val = days - (total / years)
    print(np.max(val))
    print(np.min(val))
    mp = MapPlot(sector='conus',
                 title=("OSU PRISM 2017 Days with High >= 90$^\circ$F "
                        "Departure"),
                 subtitle=("2017 thru 4 July against 1981-2016 "
                           "Year to Date Average"))
    mp.contourf(lons, lats, val, np.arange(-25, 26, 5),
                units='days', cmap=plt.get_cmap('seismic'))
    mp.postprocess(filename='test.png')

def main():
    """Do Something Fun!"""
    form = cgi.FormContent()
    ts = datetime.datetime.strptime(form["date"][0], "%Y-%m-%d")
    lat = float(form["lat"][0])
    lon = float(form["lon"][0])
    fmt = form["format"][0]
    if fmt != 'json':
        sys.stdout.write("Content-type: text/plain\n\n")
        sys.stdout.write("ERROR: Service only emits json at this time")
        return

    i, j = iemre.find_ij(lon, lat)
    offset = iemre.daily_offset(ts)

    res = {'data': [], }

    fn = "/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year,)

    sys.stdout.write('Content-type: application/json\n\n')
    if not os.path.isfile(fn):
        sys.stdout.write(json.dumps(res))
        sys.exit()

    if i is None or j is None:
        sys.stdout.write(json.dumps({'error': 'Coordinates outside of domain'}
                                    ))
        return

    nc = netCDF4.Dataset(fn, 'r')

    c2000 = ts.replace(year=2000)
    coffset = iemre.daily_offset(c2000)

    cnc = netCDF4.Dataset("/mesonet/data/iemre/mw_dailyc.nc", 'r')

    res['data'].append({
        'daily_high_f': myrounder(
           datatypes.temperature(
                nc.variables['high_tmpk'][offset, j, i], 'K').value('F'), 1),
        'climate_daily_high_f': myrounder(
           datatypes.temperature(
                cnc.variables['high_tmpk'][coffset, j, i], 'K').value("F"), 1),
        'daily_low_f': myrounder(
           datatypes.temperature(
                nc.variables['low_tmpk'][offset, j, i], 'K').value("F"), 1),
        'climate_daily_low_f': myrounder(
           datatypes.temperature(
                cnc.variables['low_tmpk'][coffset, j, i], 'K').value("F"), 1),
        'daily_precip_in': myrounder(
           nc.variables['p01d'][offset, j, i] / 25.4, 2),
        'climate_daily_precip_in': myrounder(
           cnc.variables['p01d'][coffset, j, i] / 25.4, 2),
      })
    nc.close()
    cnc.close()

    sys.stdout.write(json.dumps(res))

def test_daily_offset():
    """ Compute the offsets """
    ts = utc(2013, 1, 1, 0, 0)
    offset = iemre.daily_offset(ts)
    assert offset == 0

    ts = datetime.date(2013, 2, 1)
    offset = iemre.daily_offset(ts)
    assert offset == 31

    ts = utc(2013, 1, 5, 12, 0)
    offset = iemre.daily_offset(ts)
    assert offset == 4

def do_precip12(nc, ts):
    """Compute the 24 Hour precip at 12 UTC, we do some more tricks though"""
    offset = iemre.daily_offset(ts)
    ets = datetime.datetime(ts.year, ts.month, ts.day, 12)
    ets = ets.replace(tzinfo=pytz.timezone("UTC"))
    sts = ets - datetime.timedelta(hours=24)
    offset1 = iemre.hourly_offset(sts)
    offset2 = iemre.hourly_offset(ets)
    if ts.month == 1 and ts.day == 1:
        print(("p01d_12z  for %s [idx:%s] %s(%s)->%s(%s) SPECIAL"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        ncfn = "/mesonet/data/iemre/%s_mw_hourly.nc" % (ets.year,)
        if not os.path.isfile(ncfn):
            print("Missing %s" % (ncfn,))
            return
        hnc = netCDF4.Dataset(ncfn)
        phour = np.sum(hnc.variables['p01m'][:offset2, :, :], 0)
        hnc.close()
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            sts.year,))
        phour += np.sum(hnc.variables['p01m'][offset1:, :, :], 0)
        hnc.close()
    else:
        print(("p01d_12z  for %s [idx:%s] %s(%s)->%s(%s)"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        ncfn = "/mesonet/data/iemre/%s_mw_hourly.nc" % (ts.year,)
        if not os.path.isfile(ncfn):
            print("Missing %s" % (ncfn,))
            return
        hnc = netCDF4.Dataset(ncfn)
        phour = np.sum(hnc.variables['p01m'][offset1:offset2, :, :], 0)
        hnc.close()
    nc.variables['p01d_12z'][offset] = phour

def do_hrrr(ts):
    """Convert the hourly HRRR data to IEMRE grid"""
    total = None
    xaxis = None
    yaxis = None
    for hr in range(5, 23):  # Only need 5 AM to 10 PM for solar
        utc = ts.replace(hour=hr).astimezone(pytz.timezone("UTC"))
        fn = utc.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/model/hrrr/%H/"
                           "hrrr.t%Hz.3kmf00.grib2"))
        if not os.path.isfile(fn):
            # print 'HRRR file %s missing' % (fn,)
            continue
        grbs = pygrib.open(fn)
        try:
            if utc >= SWITCH_DATE:
                grb = grbs.select(name='Downward short-wave radiation flux')
            else:
                grb = grbs.select(parameterNumber=192)
        except ValueError:
            print 'coop/hrrr_solarrad.py %s had no solar rad' % (fn,)
            continue
        if len(grb) == 0:
            print 'Could not find SWDOWN in HRR %s' % (fn,)
            continue
        g = grb[0]
        if total is None:
            total = g.values
            lat1 = g['latitudeOfFirstGridPointInDegrees']
            lon1 = g['longitudeOfFirstGridPointInDegrees']
            llcrnrx, llcrnry = LCC(lon1, lat1)
            nx = g['Nx']
            ny = g['Ny']
            dx = g['DxInMetres']
            dy = g['DyInMetres']
            xaxis = llcrnrx + dx * np.arange(nx)
            yaxis = llcrnry + dy * np.arange(ny)
        else:
            total += g.values

    if total is None:
        print 'coop/hrrr_solarrad.py found no HRRR data for %s' % (
                                                    ts.strftime("%d %b %Y"), )
        return

    # We wanna store as W m-2, so we just average out the data by hour
    total = total / 24.0

    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year,),
                         'a')
    offset = iemre.daily_offset(ts)
    data = nc.variables['rsds'][offset, :, :]
    for i, lon in enumerate(iemre.XAXIS):
        for j, lat in enumerate(iemre.YAXIS):
            (x, y) = LCC(lon, lat)
            i2 = np.digitize([x], xaxis)[0]
            j2 = np.digitize([y], yaxis)[0]
            data[j, i] = total[j2, i2]

    nc.variables['rsds'][offset] = data
    nc.close()

def plot_daily(date, interval, plotvar, mc, mckey):
    """ Generate the plot, please """
    opts = PLOT_OPS[plotvar]
    offset = iemre.daily_offset(date)
    
    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (date.year,))
    data = nc.variables[opts['ncvar_daily']][offset] / 25.4 # inches
    lons = nc.variables['lon'][:]
    lats = nc.variables['lat'][:]
    extra = lons[-1] + (lons[-1] - lons[-2])
    lons = np.concatenate([lons, [extra,]])

    extra = lats[-1] + (lats[-1] - lats[-2])
    lats = np.concatenate([lats, [extra,]])
    x,y = np.meshgrid(lons, lats)

    nc.close()

    p = plot.MapPlot(sector='midwest',
                 title='%s IEM Reanalysis %s [%s]' % (date.strftime("%-d %b %Y"),
                                                 opts['title'], opts['units'])
                 )
    p.pcolormesh(x, y, data, opts['clevs'],
                 units=opts['units'])
    p.postprocess(web=True, memcache=mc, memcachekey=mckey, memcacheexpire=0)

def estimate_hilo(ts):
    """Estimate the High and Low Temperature based on gridded data"""
    idx = iemre.daily_offset(ts)
    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year, ),
                         'r')
    highgrid12 = temperature(nc.variables['high_tmpk_12z'][idx, :, :],
                             'K').value('F')
    lowgrid12 = temperature(nc.variables['low_tmpk_12z'][idx, :, :],
                            'K').value('F')
    highgrid00 = temperature(nc.variables['high_tmpk'][idx, :, :],
                             'K').value('F')
    lowgrid00 = temperature(nc.variables['low_tmpk'][idx, :, :],
                            'K').value('F')
    nc.close()

    for sid in nt.sts.keys():
        if nt.sts[sid]['temp24_hour'] in [0, 22, 23]:
            val = highgrid00[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        else:
            val = highgrid12[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        if val > -80 and val < 140:
            nt.sts[sid]['high'] = "%.0f" % (val, )

        if nt.sts[sid]['temp24_hour'] in [0, 22, 23]:
            val = lowgrid00[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        else:
            val = lowgrid12[nt.sts[sid]['gridj'], nt.sts[sid]['gridi']]
        if val > -80 and val < 140:
            nt.sts[sid]['low'] = "%.0f" % (val, )

def run(ts):
    ''' Actually do the work, please '''
    nc = netCDF4.Dataset('/mesonet/data/iemre/%s_mw_mrms_daily.nc' % (
                                                            ts.year,),
                         'a')
    offset = iemre.daily_offset(ts)
    ncprecip = nc.variables['p01d']
    ts += datetime.timedelta(hours=24)
    gmtts = ts.astimezone(pytz.timezone("UTC"))

    fn = gmtts.strftime(("/mesonet/ARCHIVE/data/%Y/%m/%d/GIS/q2/"
                        "p24h_%Y%m%d%H00.png"))
    img = Image.open(fn)
    data = np.asarray(img)
    # data is 3500,7000 , starting at upper L
    data = np.flipud(data)
    # Anything over 254 is bad
    res = np.where(data > 254, 0, data)
    res = np.where(np.logical_and(data >= 0, data < 100), data * 0.25, res)
    res = np.where(np.logical_and(data >= 100, data < 180),
                   25. + ((data - 100) * 1.25), res)
    res = np.where(np.logical_and(data >= 180, data < 255),
                   125. + ((data - 180) * 5.), res)

    y1 = (iemre.NORTH - mrms.SOUTH) * 100.0
    y0 = (iemre.SOUTH - mrms.SOUTH) * 100.0
    x0 = (iemre.WEST - mrms.WEST) * 100.0
    x1 = (iemre.EAST - mrms.WEST) * 100.0
    ncprecip[offset, :, :] = res[y0:y1, x0:x1]
    nc.close()

def do_coop(ts):
    """Use COOP solar radiation data"""
    pgconn = psycopg2.connect(database='coop', host='iemdb', user='nobody')
    cursor = pgconn.cursor()

    cursor.execute("""SELECT ST_x(geom), ST_y(geom),
        coalesce(narr_srad, merra_srad) from alldata a JOIN stations t
        ON (a.station = t.id) WHERE
        day = %s and t.network ~* 'CLIMATE' and substr(id, 3, 1) != 'C'
        and substr(id, 3, 4) != '0000'
    """, (ts.strftime("%Y-%m-%d"), ))
    lons = []
    lats = []
    vals = []
    for row in cursor:
        if row[2] is None or row[2] < 0:
            continue
        lons.append(row[0])
        lats.append(row[1])
        vals.append(row[2])

    nn = NearestNDInterpolator((np.array(lons), np.array(lats)),
                               np.array(vals))
    xi, yi = np.meshgrid(iemre.XAXIS, iemre.YAXIS)

    nc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_daily.nc" % (ts.year,),
                         'a')
    offset = iemre.daily_offset(ts)
    # Data above is MJ / d / m-2, we want W / m-2
    nc.variables['rsds'][offset, :, :] = nn(xi, yi) * 1000000. / 86400.
    nc.close()

def do_precip(nc, ts):
    """Compute the 6 UTC to 6 UTC precip

    We need to be careful here as the timestamp sent to this app is today,
    we are actually creating the analysis for yesterday
    """
    sts = datetime.datetime(ts.year, ts.month, ts.day, 6)
    sts = sts.replace(tzinfo=pytz.timezone("UTC"))
    ets = sts + datetime.timedelta(hours=24)
    offset = iemre.daily_offset(ts)
    offset1 = iemre.hourly_offset(sts)
    offset2 = iemre.hourly_offset(ets)
    if ts.month == 12 and ts.day == 31:
        print(("p01d      for %s [idx:%s] %s(%s)->%s(%s) SPECIAL"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            ets.year,))
        phour = np.sum(hnc.variables['p01m'][:offset2, :, :], 0)
        hnc.close()
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            sts.year,))
        phour += np.sum(hnc.variables['p01m'][offset1:, :, :], 0)
        hnc.close()
    else:
        print(("p01d      for %s [idx:%s] %s(%s)->%s(%s)"
               ) % (ts, offset, sts.strftime("%Y%m%d%H"), offset1,
                    ets.strftime("%Y%m%d%H"), offset2))
        hnc = netCDF4.Dataset("/mesonet/data/iemre/%s_mw_hourly.nc" % (
            sts.year,))
        phour = np.sum(hnc.variables['p01m'][offset1:offset2, :, :], 0)
        hnc.close()
    nc.variables['p01d'][offset] = phour

def grid_day(nc, ts):
    """
    I proctor the gridding of data on an hourly basis
    @param ts Timestamp of the analysis, we'll consider a 20 minute window
    """
    offset = iemre.daily_offset(ts)
    if ts.day == 29 and ts.month == 2:
        ts = datetime.datetime(2000, 3, 1)

    sql = """SELECT * from ncdc_climate71 WHERE valid = '%s' and
             substr(station,3,4) != '0000' and substr(station,3,1) != 'C'
             """ % (ts.strftime("%Y-%m-%d"), )
    cursor.execute(sql)
    if cursor.rowcount > 4:
        res = generic_gridder(cursor, 'high')
        if res is not None:
            nc.variables['high_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        cursor.scroll(0, mode='absolute')
        res = generic_gridder(cursor, 'low')
        if res is not None:
            nc.variables['low_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        cursor.scroll(0, mode='absolute')
        res = generic_gridder(cursor, 'precip')
        if res is not None:
            nc.variables['p01d'][offset] = res * 25.4
    else:
        print "%s has %02i entries, FAIL" % (ts.strftime("%Y-%m-%d"),
                                             cursor.rowcount)

def grid_day(nc, ts):
    """
    """
    offset = iemre.daily_offset(ts)
    print(('cal hi/lo for %s [idx:%s]') % (ts, offset))
    sql = """
       SELECT ST_x(s.geom) as lon, ST_y(s.geom) as lat, s.state,
       s.name, s.id as station,
       (CASE WHEN pday >= 0 then pday else null end) as precipdata,
       (CASE WHEN max_tmpf > -50 and max_tmpf < 130
           then max_tmpf else null end) as highdata,
       (CASE WHEN min_tmpf > -50 and min_tmpf < 95
           then min_tmpf else null end) as lowdata
       from summary_%s c, stations s WHERE day = '%s' and
       s.network in ('IA_ASOS', 'MN_ASOS', 'WI_ASOS', 'IL_ASOS', 'MO_ASOS',
        'KS_ASOS', 'NE_ASOS', 'SD_ASOS', 'ND_ASOS', 'KY_ASOS', 'MI_ASOS',
        'OH_ASOS', 'AWOS') and c.iemid = s.iemid
        """ % (ts.year, ts.strftime("%Y-%m-%d"))
    df = read_sql(sql, pgconn)

    if len(df.index) > 4:
        res = generic_gridder(df, 'highdata')
        nc.variables['high_tmpk'][offset] = datatypes.temperature(
                                                res, 'F').value('K')
        res = generic_gridder(df, 'lowdata')
        nc.variables['low_tmpk'][offset] = datatypes.temperature(
                                            res, 'F').value('K')
    else:
        print "%s has %02i entries, FAIL" % (ts.strftime("%Y-%m-%d"),
                                             cursor.rowcount)

def grid_day(nc, ts):
    """
    """
    offset = iemre.daily_offset(ts)
    icursor.execute("""
       SELECT ST_x(s.geom) as lon, ST_y(s.geom) as lat, 
       (CASE WHEN pday >= 0 then pday else null end) as precipdata,
       (CASE WHEN max_tmpf > -50 and max_tmpf < 130 then max_tmpf else null end) as highdata,
       (CASE WHEN min_tmpf > -50 and min_tmpf < 95 then min_tmpf else null end) as lowdata 
       from summary_%s c, stations s WHERE day = '%s' and 
       s.network in ('IA_ASOS', 'MN_ASOS', 'WI_ASOS', 'IL_ASOS', 'MO_ASOS',
        'KS_ASOS', 'NE_ASOS', 'SD_ASOS', 'ND_ASOS', 'KY_ASOS', 'MI_ASOS',
        'OH_ASOS', 'AWOS') and c.iemid = s.iemid
        """ % (ts.year, ts.strftime("%Y-%m-%d")))

    if icursor.rowcount > 4:
        res = generic_gridder(icursor, 'highdata')
        nc.variables['high_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        icursor.scroll(0, mode='absolute')
        res = generic_gridder(icursor, 'lowdata')
        nc.variables['low_tmpk'][offset] = datatypes.temperature(res, 'F').value('K')
        icursor.scroll(0, mode='absolute')
        #res = generic_gridder(icursor, 'precipdata')
        #nc.variables['p01d'][offset] = res * 25.4
    else:
        print "%s has %02i entries, FAIL" % (ts.strftime("%Y-%m-%d"), 
            icursor.rowcount)

def do_var(varname):
    """
    Run our estimator for a given variable
    """
    currentnc = None
    sql = """select day, station from alldata_%s WHERE %s is null
        and day >= '1893-01-01' ORDER by day ASC""" % (state.lower(), varname)
    ccursor.execute(sql)
    for row in ccursor:
        day = row[0]
        station = row[1]
        if station not in nt.sts:
            continue

        sql = """
            SELECT station, %s from alldata_%s WHERE %s is not NULL
            and station in %s and day = '%s'
            """ % (varname, state, varname, tuple(friends[station]), day)
        ccursor2.execute(sql)
        weight = []
        value = []
        for row2 in ccursor2:
            idx = friends[station].index(row2[0])
            weight.append(weights[station][idx])
            value.append(row2[1])

        if len(weight) < 3:
            # Nearest neighbors failed, so lets look at our grided analysis
            # and sample from it
            if currentnc is None or currentnc.title.find(str(day.year)) == -1:
                currentnc = netCDF4.Dataset(("/mesonet/data/iemre/"
                                             "%s_mw_daily.nc") % (day.year,))
            tidx = iemre.daily_offset(datetime.datetime(day.year, day.month,
                                                        day.day))
            iidx, jidx = iemre.find_ij(nt.sts[station]['lon'],
                                       nt.sts[station]['lat'])
            iemreval = currentnc.variables[vnameconv[varname]][tidx, jidx,
                                                               iidx]
            if varname in ('high', 'low'):
                interp = temperature(iemreval, 'K').value('F')
            else:
                interp = distance(iemreval, 'MM').value('IN')
            print '--> Neighbor failure, %s %s %s' % (station, day, varname)
        else:
            mass = sum(weight)
            interp = np.sum(np.array(weight) * np.array(value) / mass)

        dataformat = '%.2f'
        if varname in ['high', 'low']:
            dataformat = '%.0f'
        print(('Set station: %s day: %s varname: %s value: %s'
               ) % (station, day, varname, dataformat % (interp,)))
        sql = """
            UPDATE alldata_%s SET estimated = true, %s = %s WHERE
            station = '%s' and day = '%s'
            """ % (state.lower(), varname,
                   dataformat % (interp,), station, day)
        sql = sql.replace(' nan ', ' null ')
        ccursor2.execute(sql)