Python RPN.get_all_time_records_for_name方法代码示例

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]
def main():
    #path = "/b2_fs2/huziy/OMSC26_Can_long_new_v01/pm1958010100_02275344p"

    path = "/b2_fs2/huziy/OMSC26_Can_long_new_v01/pm1958010100_00008640p"
    r = RPN(path)

    res_date = c_int()
    res_time = c_int()
    mode = c_int(-3)

    #test1
    dateo = 10158030
    r._dll.newdate_wrapper(byref(c_int(dateo)), byref(res_date), byref(res_time), byref(mode))
    s_date = "{0:08d}{1:08d}".format(res_date.value, res_time.value)
    print("stamp: {0:09d}, result: {1}".format(dateo, s_date))

    dateo = 488069900
    r._dll.newdate_wrapper(byref(c_int(dateo)), byref(res_date), byref(res_time), byref(mode))
    s_date = "{0:08d}{1:08d}".format(res_date.value, res_time.value)
    print("stamp: {0:09d}, result: {1}".format(dateo, s_date))

    dateo = 1069261100
    r._dll.newdate_wrapper(byref(c_int(dateo)), byref(res_date), byref(res_time), byref(mode))
    s_date = "{0:08d}{1:08d}".format(res_date.value, res_time.value)
    print("stamp: {0:09d}, result: {1}".format(dateo, s_date))

    dateo = 632053700
    r._dll.newdate_wrapper(byref(c_int(dateo)), byref(res_date), byref(res_time), byref(mode))
    s_date = "{0:08d}{1:08d}".format(res_date.value, res_time.value)
    print("stamp: {0:09d}, result: {1}".format(dateo, s_date))

    ts = r.get_all_time_records_for_name("TRAF")

    times = list(sorted(ts.keys()))
    print(times[:20])
    print(times[0], times[-1])
    r.close()

    folderPath = "/b2_fs2/huziy/OMSC26_ERA40I_long_new_v02/"
    for fName in os.listdir(folderPath):
        if not fName.startswith("pm"): continue

        fPath = os.path.join(folderPath, fName)

        r = RPN(fPath)
        r.suppress_log_messages()

        data = r.get_all_time_records_for_name(varname="TDRA")
        r.close()

        print(fName)
        print(sorted(data.keys())[:5])
        print(25 * "*")
        input("press any key")

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]
def main():
    #path = "/home/huziy/skynet3_rech1/test/snw_LImon_NA_CRCM5_CanESM2_historical_r1i1p1_185001-200512.rpn"
    path = "/home/sheena/skynet3_exec2/RPN/src/permafrost/snw_NA_CRCM5_CanESM2_rcp45_r1i1p1_200601-210012.rpn"
    months = [1,2,12]

    varname = "I5"

    rObj = RPN( path )
    records = rObj.get_all_time_records_for_name(varname=varname)
    lons2d, lats2d = rObj.get_longitudes_and_latitudes()

    rObj.close()



    times = sorted(records.keys())
    vals =  np.array( [records[t] for t in times])

    year_range = list(range(2006, 2101))
    nc_file_name = "{0:s}_{1:d}_{2:d}.nc".format(varname, year_range[0], year_range[-1])
    nx, ny = vals[0].shape


    #create netcdf file
    ds = Dataset(nc_file_name, "w", format = 'NETCDF3_CLASSIC')
    ds.createDimension('lon', nx)
    ds.createDimension('lat', ny)
    ds.createDimension("year", len(year_range))
    the_var = ds.createVariable(varname, 'f', ("year",'lat','lon'))
    the_lon = ds.createVariable("xlon", 'f', ('lat','lon'))
    the_lat = ds.createVariable("xlat", 'f', ('lat','lon'))



    for i, the_year in enumerate(year_range):
        bool_vector = [t.year == the_year and t.month in months for t in times]
        bool_vector = np.array(bool_vector)
        the_var[i,:,:] = np.mean(vals[bool_vector], axis=0).transpose()

    the_lon[:] = lons2d[:,:].transpose()
    the_lat[:] = lats2d[:,:].transpose()

    ds.close()










    #TODO: implement
    pass

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]
    def get_seasonal_mean_for_year_of_2d_var(self, the_year, months=None, var_name=""):
        """
        Return mean over months of a given 2d field
        returns numpy array of dimensions (x, y)
        """
        monthly_means = []
        for the_month in months:

            key = (the_year, the_month)
            if key not in self.yearmonth_to_data_path:
                print(("Warning donot have data for {0}/{1}".format(the_year, the_month)))
                continue

            path = self.yearmonth_to_data_path[key]
            rpn_obj = RPN(path)
            records = rpn_obj.get_all_time_records_for_name(varname=var_name)
            monthly_means.append(np.mean(list(records.values()), axis=0))
            rpn_obj.close()

        print((the_year, np.min(np.mean(monthly_means, axis=0)), np.max(np.mean(monthly_means, axis=0))))
        return np.mean(monthly_means, axis=0)

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]
def plot_alt_from_recent_jpp_and_andrey():
    driver = "CanESM"
    path = "/home/huziy/skynet1_rech3/jpp_progs_recent/pf/Arctic_Andrey_ALT_canesm1"
    rObj = RPN(path)

    altt = rObj.get_all_time_records_for_name(varname="FALT")
    alt = np.mean( np.array(list(altt.values())), axis = 0)
    rObj.close()


    coord_file = path 
    basemap, lons2d, lats2d = draw_regions.get_basemap_and_coords(resolution="c",
        file_path = coord_file, lon1 = 60, lat1 = 89.999, lon2 = -30, lat2 = 0.00001, 
        projection = "omerc", round = True
    )
    assert isinstance(basemap, Basemap)

    lons2d[lons2d > 180] -= 360

    x, y = basemap(lons2d, lats2d)


    #permafrost_mask = draw_regions.get_permafrost_mask(lons2d, lats2d)
    #mask_cond = (permafrost_mask <= 0) | (permafrost_mask >= 3)
    alt = np.ma.masked_where(alt < 0, alt)

    fig = plt.figure()
    assert isinstance(fig, Figure)

    h_max = 10
    
    bounds = list(range(6))
    cmap = my_colormaps.get_lighter_jet_cmap(ncolors=len(bounds) - 1) #cm.get_cmap("jet",10)
    #bounds = [0,0.1,0.5,1,2,3,5,8,9,10,11]
    norm = BoundaryNorm(boundaries=bounds,ncolors=len(bounds)-1, clip = False)
    #norm = None

    cmap.set_over(cmap(1.0))
    clevels = np.arange(0,h_max+1,1)
    gs = gridspec.GridSpec(1,1)



    #ax = fig.add_axes([0,0,1,1], polar = True)
    ax = fig.add_subplot(gs[0,0])
    assert isinstance(ax, Axes)
    #hc = np.ma.masked_where(mask_cond | (np.min(altt.values(), axis = 0) < 0), alt)
    #hc = np.ma.masked_where( (hc < 0), hc)
    img = basemap.pcolormesh(x, y, alt, cmap = cmap, vmax = h_max, norm=norm)
    ax.set_title("ALT, JPP-Andrey, ({0} - CRCM5)  \n".format(driver))

    basemap.drawcoastlines(ax = ax, linewidth=0.5)
    #basemap.readshapefile("data/pf_4/permafrost8_wgs84/permaice", name="zone",
    #        ax=ax, linewidth=1.5)

    basemap.readshapefile("data/permafrost_lat-lon1/permafrost_latlon", name="zone",
            ax=ax, linewidth=1.5)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", "5%", pad="3%")
    cb = fig.colorbar(img,  cax = cax, extend = "max", ticks = bounds)
    cax.set_title("m \n")


    fig.tight_layout(w_pad=0.0)


    fig.savefig("alt_jpp_current_{0}.png".format(driver))

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]
def plot_mean_alt_from_jpp_results():
    start_year = 1981
    end_year = 2008
    path = "data/alts_by_jpp/MonTS_NA_ERA40_ALT_{0}_{1}".format(start_year, end_year)

    rObj = RPN(path)

    altt = rObj.get_all_time_records_for_name(varname="FALT")
    alt = np.mean( np.array(list(altt.values())), axis = 0)
    rObj.close()



    sim_data_folder = "/home/huziy/skynet1_rech3/cordex/CORDEX_DIAG/era40_driven_b1"
    coord_file = os.path.join(sim_data_folder, "pmNorthAmerica_0.44deg_ERA40-Int_B1_200812_moyenne")
    basemap, lons2d, lats2d = draw_regions.get_basemap_and_coords(resolution="c",
        file_path = coord_file, llcrnrlat=40.0, llcrnrlon=-145, urcrnrlon=-20, urcrnrlat=74
    )
    assert isinstance(basemap, Basemap)

    lons2d[lons2d > 180] -= 360

    x, y = basemap(lons2d, lats2d)


    permafrost_mask = draw_regions.get_permafrost_mask(lons2d, lats2d)
    mask_cond = (permafrost_mask <= 0) | (permafrost_mask >= 3)
    alt = np.ma.masked_where(mask_cond, alt)

    fig = plt.figure()
    assert isinstance(fig, Figure)


    h_max = 10
    cmap = my_colormaps.get_lighter_jet_cmap(ncolors=10) #cm.get_cmap("jet",10)
    bounds = [0,0.1,0.5,1,2,3,5,8,9,10,11]
    norm = BoundaryNorm(boundaries=bounds,ncolors=len(bounds), clip=True)
    #norm = None

    cmap.set_over(cmap(1.0))
    clevels = np.arange(0,h_max+1,1)
    gs = gridspec.GridSpec(3,1)

    ax = fig.add_subplot(gs[0,0])
    assert isinstance(ax, Axes)
    hc = np.ma.masked_where(mask_cond | (np.min(list(altt.values()), axis = 0) < 0), alt)
    #hc = np.ma.masked_where( (hc < 0), hc)
    img = basemap.pcolormesh(x, y, hc, cmap = cmap, vmax = h_max, norm=norm)
    ax.set_title("ALT, JPP ({0} - {1}) \n".format(start_year, end_year))

    basemap.drawcoastlines(ax = ax, linewidth=0.5)
    basemap.readshapefile("data/pf_4/permafrost8_wgs84/permaice", name="zone",
            ax=ax, linewidth=1.5)


    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", "5%", pad="3%")
    cb = fig.colorbar(img,  cax = cax, extend = "max", ticks = bounds)
    cax.set_title("m \n")


    fig.tight_layout(w_pad=0.0)


    fig.savefig("alt_jpp_current.png")

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]

#.........这里部分代码省略.........
    r1 = RPN(path1)
    soiltemp1 = r1.get_4d_field(name=soiltemp_var_name)
    lons2d, lats2d = r1.get_longitudes_and_latitudes_for_the_last_read_rec()
    rll = RotatedLatLon(**r1.get_proj_parameters_for_the_last_read_rec())
    b = rll.get_basemap_object_for_lons_lats(lons2d=lons2d, lats2d=lats2d, resolution="c")
    r1.close()

    r2 = RPN(path2)
    soiltemp2 = r2.get_4d_field(name=soiltemp_var_name)
    r2.close()

    dates_sorted = list(sorted(soiltemp1.keys()))
    levels_sorted = list(sorted(list(soiltemp1.items())[0][1].keys()))

    # group dates for each year
    alt1_list = []
    alt2_list = []
    for year, dates_group in itertools.groupby(dates_sorted, lambda par: par.year):
        t1 = []
        t2 = []
        for d in dates_group:
            t1.append([soiltemp1[d][lev] for lev in levels_sorted])
            t2.append([soiltemp2[d][lev] for lev in levels_sorted])

        # Get maximum temperature profiles
        t1max = np.max(t1, axis=0).transpose(1, 2, 0)
        t2max = np.max(t2, axis=0).transpose(1, 2, 0)

        print(t1max.shape, t2max.shape)

        #calculate and save alts
        h1 = crcm_data_manager.get_alt(t1max)
        h2 = crcm_data_manager.get_alt(t2max)

        h1[h1 < 0] = np.nan
        h2[h2 < 0] = np.nan

        alt1_list.append(h1)
        alt2_list.append(h2)

    #take into account permafrost
    alt1_list_pf = []
    alt2_list_pf = []
    n_years_for_pf = 3
    for i in range(len(alt1_list) - n_years_for_pf):
        alt1_list_pf.append(np.max(alt1_list[i:i+n_years_for_pf], axis=0))
        alt2_list_pf.append(np.max(alt2_list[i:i+n_years_for_pf], axis=0))



    # calculate climatological mean
    alt1 = np.mean(alt1_list_pf, axis=0)
    alt2 = np.mean(alt2_list_pf, axis=0)

    print(np.isnan(alt1).any(), np.isnan(alt2).any())

    #mask nans
    alt1 = np.ma.masked_where(np.isnan(alt1), alt1)
    alt2 = np.ma.masked_where(np.isnan(alt2), alt2)

    #calculate change due to fixed depth to bedrock
    delta = alt2 - alt1

    #
    for i in range(len(alt1_list)):
        alt1_list[i] = np.ma.masked_where(np.isnan(alt1_list[i]), alt1_list[i])
        alt2_list[i] = np.ma.masked_where(np.isnan(alt2_list[i]), alt2_list[i])

    #tval, pval = ttest_ind(alt1_list, alt2_list)

    #mask oceans
    lons2d[lons2d > 180] -= 360
    dpth = maskoceans(lons2d, lats2d, dpth)
    delta = np.ma.masked_where(dpth.mask, delta)


    #calculate differences in SWE
    path_swe_1 = path1.replace("TBAR.rpn", "SNO.rpn")
    r1 = RPN(path_swe_1)
    swe1 = r1.get_all_time_records_for_name("SNO")
    r1.close()
    swe1_winter_clim = np.mean(
        [field for key, field in swe1.items() if key.month in [1, 2, 12]], axis=0)
    swe1_winter_clim = np.ma.masked_where((swe1_winter_clim >= 999) | dpth.mask, swe1_winter_clim)

    path_swe_2 = path2.replace("TBAR.rpn", "SNO.rpn")
    r2 = RPN(path_swe_2)
    swe2 = r2.get_all_time_records_for_name("SNO")
    r2.close()
    swe2_winter_clim = np.mean(
        [field for key, field in swe2.items() if key.month in [1, 2, 12]], axis=0)
    swe2_winter_clim = np.ma.masked_where((swe2_winter_clim >= 999) | dpth.mask, swe2_winter_clim)




    #plotting
    print("Start plotting ...")
    plot_differences(b, lons2d, lats2d, dpth, delta, label="ALT(DPTH=3.6m) - ALT(DPTH~)", pvalue=None,
                     swe_diff=swe2_winter_clim - swe1_winter_clim)

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]
def main():
    path = "/RESCUE/skynet3_rech1/huziy/CNRCWP/Calgary_flood/Global_NA_v1/Samples/Global_NA_v1_201306/pm2013010100_00017280p"
    varname = "PR"
    plot_units = "mm/day"
    mult_coeff = 1000 * 24 * 3600
    add_offset = 0

    img_folder = "/RESCUE/skynet3_rech1/huziy/CNRCWP/Calgary_flood/glob_sim/{}/monthly/{}".format(varname, os.path.basename(path))
    img_folder = Path(img_folder)
    if not img_folder.is_dir():
        img_folder.mkdir(parents=True)

    r = RPN(path=path)

    pr = r.get_all_time_records_for_name(varname=varname)

    lons2d, lats2d = r.get_longitudes_and_latitudes_for_the_last_read_rec()
    rll = RotatedLatLon(**r.get_proj_parameters_for_the_last_read_rec())



    bmp = rll.get_basemap_object_for_lons_lats(lons2d=lons2d,
                                               lats2d=lats2d,
                                               resolution="c", no_rot=True)

#    bmp = Basemap(projection="robin", lon_0=180)

    xx, yy = bmp(lons2d, lats2d)


    dates = list(sorted(pr.keys()))
    data = np.array([pr[d] for d in dates])

    p = pd.Panel(data=data, items=dates, major_axis=range(data.shape[1]), minor_axis=range(data.shape[2]))

    # p_daily = p.groupby(lambda d: d.day, axis="items").mean()
    p_daily = p.apply(np.mean, axis="items")

    print(p_daily.head())

    lons2d[lons2d > 180] -= 360

    bmap_params = bmp.projparams
    bmap_params.update({
        'llcrnrlon': lons2d[0, 0], 'urcrnrlon': lons2d[-1, -1], 'llcrnrlat': lats2d[0, 0], 'urcrnrlat': lats2d[-1, -1]
    })
    rpole_crs = ccrs.RotatedPole(pole_longitude=bmap_params["lon_0"] + 180,
                                 pole_latitude=bmap_params["o_lat_p"])



    clevs = [0, 0.01, 0.1, 1, 1.5, 2, 5, 10, 20, 40, 60, 80]
    norm = BoundaryNorm(clevs, ncolors=len(clevs) - 1)
    cmap = cm.get_cmap(cm_basemap.s3pcpn, len(clevs) - 1)
    field = p_daily.values * mult_coeff + add_offset
    fig = plt.figure()
    plt.title("{}, {}/{}/{}, {}".format(varname, 1, dates[0].month, dates[0].year, plot_units))




    ax = plt.axes(projection=rpole_crs)
    ax.coastlines(resolution='110m')

    ax.gridlines()
    ax.gridlines()


    # cs = bmp.contourf(xx, yy, field, clevs, norm=norm, extend="max", cmap=cmap)
    cs = ax.pcolormesh(lons2d[:-1, :-1], lats2d[:-1, :-1], field[:-1, :-1], norm=norm, cmap=cmap, transform=rpole_crs)
    plt.colorbar(cs, ticks=clevs, extend="max", ax=ax)
    img_file = img_folder.joinpath("{:02d}-{:02d}-{}.png".format(1, dates[0].month, dates[0].year))
    # bmp.drawcoastlines()
    # bmp.drawstates()
    # bmp.drawcounties()
    # bmp.drawcountries()

    plt.savefig(img_file.open("wb"))
    plt.close(fig)
    print(pr[dates[0]].mean())

# 需要导入模块: from rpn.rpn import RPN [as 别名]
# 或者: from rpn.rpn.RPN import get_all_time_records_for_name [as 别名]
def plot_mean(data_path = "", file_prefix = "pm"):
    """
    """
    var_name = "CLDP"
    means = []
    means_dict = {}
    for file in os.listdir(data_path):
        if not file.startswith(file_prefix): continue
        file_path = os.path.join(data_path, file)

        r = RPN(file_path)
        levels = r.get_all_time_records_for_name(varname=var_name)
        means_dict.update(levels)
        means.append(np.mean( list(levels.values()) , axis = 0))
    mean_level = np.array(means).mean(axis = 0)
    b, lons2d, lats2d = draw_regions.get_basemap_and_coords(
        file_path="data/from_guillimin/vary_lake_level1/pm1985010100_00000000p",
        lon1=-68, lat1=52, lon2=16.65, lat2=0
    )

    std_level = np.std( np.array(means), axis=0)


    #plot mean
    plt.figure()
    _plot_depth(mean_level, lons2d, lats2d, basemap=b, clevels=range(0,310, 10))
    plt.savefig("mean_lake_levels.png")


    #plot std
    plt.figure()
    _plot_depth(std_level, lons2d, lats2d, basemap=b, clevels=np.arange(0,2.2, 0.2), lowest_value=1e-3)
    plt.savefig("std_lake_levels.png")



    #plot initial lake depth
    plt.figure()
    plot_initial_lake_depth(lons2d=lons2d, lats2d=lats2d, basemap=b)
    plt.savefig("initial_lake_depth.png")

    #plot lake fraction
    plt.figure()
    plot_lake_fraction(lons2d=lons2d, lats2d=lats2d, basemap=b)
    plt.savefig("lake_fraction.png")

    fig = plt.figure()
    assert isinstance(fig, Figure)
    gs = GridSpec(2,2)

    seasons = ["(a) Winter (DJF)", "(b) Spring (MAM)", "(c) Summer (JJA)", "(d) Fall (SON)"]
    months_list = [[12,1,2], [3,4,5], [6,7,8], [9,10,11]]

    i = 0
    for season, months in zip(seasons, months_list):
        ax = fig.add_subplot(gs[i // 2, i % 2])
        assert isinstance(ax, Axes)
        data = _get_seasonal_mean_anomaly(means_dict, mean_level, months=months)
        the_max = np.round(np.max(np.abs(data)) * 10 ) / 10.0

        data = np.ma.masked_where(mean_level < 0.1, data)

        _plot_depth(data, lons2d, lats2d, ax = ax, basemap=b,
            clevels=np.arange(-0.5, 0.6, 0.1),
            lowest_value=0.001
        )
        ax.set_title(season)
        i += 1
    fig.tight_layout()
    fig.savefig("seasonal_anomalies.png")