Python Basemap.rotate_vector方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
def main():

    plot_utils.apply_plot_params(width_cm=20, height_cm=20, font_size=10)

    high_hles_years = [1993, 1995, 1998]
    low_hles_years = [1997, 2001, 2006]
    data_path = "/BIG1/skynet1_rech1/diro/sample_obsdata/eraint/eraint_uvslp_years_198111_201102_NDJmean_ts.nc"





    with xr.open_dataset(data_path) as ds:
        print(ds)


        u = get_composit_for_name(ds, "u10", high_years_list=high_hles_years, low_years_list=low_hles_years)
        v = get_composit_for_name(ds, "v10", high_years_list=high_hles_years, low_years_list=low_hles_years)
        msl = get_composit_for_name(ds, "msl", high_years_list=high_hles_years, low_years_list=low_hles_years)

        lons = ds["longitude"].values
        lats = ds["latitude"].values

        print(lats)
        print(msl.shape)
        print(lons.shape, lats.shape)


        lons2d, lats2d = np.meshgrid(lons, lats)


    fig = plt.figure()

    map = Basemap(llcrnrlon=-130, llcrnrlat=22, urcrnrlon=-28,
                  urcrnrlat=65, projection='lcc', lat_1=33, lat_2=45,
                  lon_0=-95, resolution='i', area_thresh=10000)


    clevs = np.arange(-11.5, 12, 1)
    cmap = cm.get_cmap("bwr", len(clevs) - 1)
    bn = BoundaryNorm(clevs, len(clevs) - 1)

    x, y = map(lons2d, lats2d)
    im = map.contourf(x, y, msl / 100, levels=clevs, norm=bn, cmap=cmap) # convert to mb (i.e hpa)
    map.colorbar(im)


    stride = 2
    ux, vy = map.rotate_vector(u, v, lons2d, lats2d)
    qk = map.quiver(x[::stride, ::stride], y[::stride, ::stride], ux[::stride, ::stride], vy[::stride, ::stride],
               scale=10, width=0.01, units="inches")
    plt.quiverkey(qk, 0.5, -0.1, 2, "2 m/s", coordinates="axes")


    map.drawcoastlines(linewidth=0.5)
    map.drawcountries()
    map.drawstates()
    #plt.show()

    fig.savefig("hles_wind_compoosits.png", bbox_inches="tight", dpi=300)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
 def test_npstere(self):
     # NP Stereographic case
     B=Basemap(projection='npstere', boundinglat=50., lon_0=0.)
     u,v,lat,lon=self.make_array()
     v = np.ones((len(lat), len(lon)))
     ru, rv = B.rotate_vector(u,v, lon, lat)
     assert_almost_equal(ru[2, :],[1,-1,-1,1], 6)
     assert_almost_equal(rv[2, :],[1,1,-1,-1], 6)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
 def test_cylindrical(self):
     # Cylindrical case
     B = Basemap()
     u,v,lat,lon=self.make_array()
     ru, rv = B.rotate_vector(u,v, lon, lat)
     # Check that the vectors are identical.
     assert_almost_equal(ru, u)
     assert_almost_equal(rv, v)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
 def test_nan(self):
     B = Basemap()
     u,v,lat,lon=self.make_array()
     # Set one element to 0, so that the vector magnitude is 0.
     u[1,1] = 0.
     ru, rv = B.rotate_vector(u,v, lon, lat)
     assert not np.isnan(ru).any()
     assert_almost_equal(u, ru)
     assert_almost_equal(v, rv)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
def test():
    plt.figure()
    b = Basemap()
    u = np.array([1, ])
    v = np.array([1, ])
    lon, lat = np.array([-90, ]), np.array([45, ])
    xx, yy = b(lon, lat)
    print(xx.shape)
    b.quiver(xx, yy, u, v, color="r")

    urot, vrot = b.rotate_vector(u, v, lon, lat)
    b.quiver(xx, yy, urot, vrot, color="g")

    b.drawcoastlines()



    # Plot the same in rotpole projection

    HL_LABEL = "CRCM5_HL"
    NEMO_LABEL = "CRCM5_NEMO"

    sim_label_to_path = OrderedDict(
        [(HL_LABEL, "/RESCUE/skynet3_rech1/huziy/CNRCWP/C5/2016/2-year-runs/coupled-GL+stfl_oneway/Samples"),
         (NEMO_LABEL, "/HOME/huziy/skynet3_rech1/CNRCWP/C5/2016/2-year-runs/coupled-GL+stfl/Samples")]
    )

    # get a coord file ...
    coord_file = ""
    found_coord_file = False
    for mdir in os.listdir(sim_label_to_path[HL_LABEL]):

        mdir_path = os.path.join(sim_label_to_path[HL_LABEL], mdir)
        if not os.path.isdir(mdir_path):
            continue

        for fn in os.listdir(mdir_path):
            print(fn)
            if fn[:2] not in ["pm", "dm", "pp", "dp"]:
                continue

            coord_file = os.path.join(mdir_path, fn)
            found_coord_file = True

        if found_coord_file:
            break

    bmp, lons, lats = nemo_hl_util.get_basemap_obj_and_coords_from_rpn_file(path=coord_file)


    plt.figure()
    urot, vrot = bmp.rotate_vector(u, v, lon, lat)
    xx, yy = bmp(lon, lat)

    bmp.quiver(xx, yy, urot, vrot, color="b")
    bmp.drawcoastlines()
    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
71964,60.7,-135.1,5310.0,250.0,36.5
'''
    # Extracts data from IDV output and masks out stations outside of North America
    lat,lon,height,wdir,speed = np.loadtxt(StringIO(data), delimiter=',',
        unpack=True, usecols=(1,2,3,4,5))

    # Create a map for plotting
    bm = Basemap(projection='tmerc', lat_0=90.0, lon_0=-100.0, lat_ts=40.0,
        llcrnrlon=-121, llcrnrlat=24, urcrnrlon=-65, urcrnrlat=46,
        resolution='l')

    # Get U,V components from wind speed and direction
    u,v = metpy.get_wind_components(speed, wdir)

    # Rotate the vectors to be properly aligned in the map projection
    u,v = bm.rotate_vector(u, v, lon, lat)

    # Generate grid of x,y positions
    lon_grid, lat_grid, x_grid, y_grid = bm.makegrid(130, 60, returnxy=True)

    # Transform the obs to basemap space for gridding
    obx,oby = bm(lon, lat)

    # Perform analysis of height obs using Cressman weights
    heights_oban = grid_data(height, x_grid, y_grid, obx, oby,
        obans[which_oban][0], obans[which_oban][1])

    heights_oban = maskoceans(lon_grid, lat_grid, heights_oban)

    # Map plotting
    contours = np.arange(5000., 5800., 60.0)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
lats1 = -90.+dellat*np.arange(nlats)
lons1 = -180.+dellon*np.arange(nlons)
lons, lats = np.meshgrid(lons1, lats1)

# plot vectors in geographical (lat/lon) coordinates.

# north polar projection.
m = Basemap(lon_0=-135,boundinglat=25,round=True,
            resolution='c',area_thresh=10000.,projection='npstere')
# create a figure, add an axes.
fig=plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1,0.1,0.7,0.7])
# rotate wind vectors to map projection coordinates.
# (also compute native map projections coordinates of lat/lon grid)
# only do Northern Hemisphere.
urot,vrot,x,y = m.rotate_vector(u[36:,:],v[36:,:],lons[36:,:],lats[36:,:],returnxy=True)
# plot filled contours over map.
cs = m.contourf(x,y,p[36:,:],15,cmap=plt.cm.jet)
# plot wind vectors over map.
Q = m.quiver(x,y,urot,vrot) #or specify, e.g., width=0.003, scale=400)
qk = plt.quiverkey(Q, 0.95, 1.05, 25, '25 m/s', labelpos='W')
m.colorbar(pad='12%') # draw colorbar
m.drawcoastlines()
m.drawcountries()
# draw parallels
delat = 20.
circles = np.arange(0.,90.+delat,delat).tolist()+\
          np.arange(-delat,-90.-delat,-delat).tolist()
m.drawparallels(circles,labels=[1,1,1,1])
# draw meridians
delon = 45.

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
def main():
    start_year = 1979
    end_year = 1981

    HL_LABEL = "CRCM5_HL"
    NEMO_LABEL = "CRCM5_NEMO"

    dx = 0.1
    dy = 0.1

    file_prefix = "pm"
    PR_level = -1
    PR_level_type = level_kinds.ARBITRARY

    tprecip_vname = "PR"
    sprecip_vname = "SN"

    TT_level = 1
    TT_level_type = level_kinds.HYBRID

    sim_label_to_path = OrderedDict(
        [(HL_LABEL, "/RESCUE/skynet3_rech1/huziy/CNRCWP/C5/2016/2-year-runs/coupled-GL+stfl_oneway/Samples"),
         (NEMO_LABEL, "/HOME/huziy/skynet3_rech1/CNRCWP/C5/2016/2-year-runs/coupled-GL+stfl/Samples")]
    )

    # get a coord file ... (use pm* files, since they contain NEM1 variable)
    # Should be NEMO_LABEL, since the hostetler case does not calculate NEM? vars
    coord_file = ""
    found_coord_file = False
    for mdir in os.listdir(sim_label_to_path[NEMO_LABEL]):

        mdir_path = os.path.join(sim_label_to_path[NEMO_LABEL], mdir)
        if not os.path.isdir(mdir_path):
            continue

        for fn in os.listdir(mdir_path):

            if fn[:2] not in ["pm", ]:
                continue

            if fn[-9:-1] == "0" * 8:
                continue

            coord_file = os.path.join(mdir_path, fn)
            found_coord_file = True

        if found_coord_file:
            break

    bmp, lons, lats = nemo_hl_util.get_basemap_obj_and_coords_from_rpn_file(path=coord_file)
    xx, yy = bmp(lons, lats)


    stride = 3

    #
    rot_path = "/HOME/huziy/skynet3_rech1/CNRCWP/C5/2016/2-year-runs/coupled-GL+stfl/rotated_wind_CRCM5_NEMO.nc"
    fig = plt.figure()

    with Dataset(rot_path) as ds:
        ncvars = ds.variables
        uu_rot, vv_rot = ncvars["UU"][10, ...], ncvars["VV"][10, ...]

    plt.title("rotated in the file")
    bmp.quiver(xx[::stride, ::stride], yy[::stride, ::stride], uu_rot[::stride, ::stride], vv_rot[::stride, ::stride], scale=1000, color="r")


    lons[lons > 180] -= 360

    not_rot_path = "/HOME/huziy/skynet3_rech1/CNRCWP/C5/2016/2-year-runs/coupled-GL+stfl/not_rotated_wind_CRCM5_NEMO.nc"
    with Dataset(not_rot_path) as ds:
        ncvars = ds.variables
        uu, vv = ncvars["UU"][10, ...], ncvars["VV"][10, ...]

        uu_rot1, vv_rot1 = rotate_vecs_from_geo_to_rotpole(uu, vv, lons, lats, bmp=bmp)


    plt.title("not rotated in the file")
    bmp.quiver(xx[::stride, ::stride], yy[::stride, ::stride], uu_rot1[::stride, ::stride], vv_rot1[::stride, ::stride],
               scale=1000)
    bmp.drawcoastlines()





    plt.figure()
    b = Basemap(lon_0=0)
    xx1, yy1 = b(lons, lats)

    uu_rot2, vv_rot2 = b.rotate_vector(uu, vv, lons, lats)
    b.quiver(xx1[::stride, ::stride], yy1[::stride, ::stride], uu_rot2[::stride, ::stride], vv_rot2[::stride, ::stride],
               scale=1000)

    b.drawcoastlines()

    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]

#.........这里部分代码省略.........
        alpha = kw.pop('alpha', 0.5)
        plot_obj = m.contour(X, Y, data, cmap=cmap, colors=colors,
            levels=levels, alpha=alpha, **kw)
        label_contour = kw.pop('label_contour', False)
        if label_contour:
            plt.clabel(plot_obj,plot_obj.levels[::2],fmt='%.2G')

    # contourf + contour
    elif plot_type in ('contourf+',):
        if proj in ('ortho','npstere', 'nplaea', 'npaeqd', 'spstere',
            'splaea', 'spaeqd')\
            and np.isclose(np.abs(lon_edge[-1]-lon_edge[0]), 360):
            data, lon = addcyclic(data, lon)
            Lon, Lat = np.meshgrid(lon,lat)
            X, Y = m(Lon, Lat)
        extend = kw.pop('extend', 'both')
        linewidths = kw.pop('linewidths', 1)
        plot_obj = m.contourf(X, Y, data, extend=extend, cmap=cmap,
            levels=levels, **kw)
        colors = kw.pop('colors', 'k')
        if colors is not None:
            cmap = None
        alpha = kw.pop('alpha', 0.5)
        m.contour(X, Y, data, cmap=cmap, colors=colors, alpha=alpha,
            levels=levels, linewidths=linewidths, **kw)

    # quiverplot
    elif plot_type in ('quiver',):
        stride = kw.pop('stride', 1)
        stride_lon = kw.pop('stride_lon', stride)
        stride_lat = kw.pop('stride_lat', stride)
        lon_ = lon[::stride_lon] # subset of lon
        lat_ = lat[::stride_lat]
        u_ = u[::stride_lat, ::stride_lon]
        v_ = v[::stride_lat, ::stride_lon]
        Lon_, Lat_ = np.meshgrid(lon_, lat_)
        u_rot, v_rot, X_, Y_ = m.rotate_vector(
            u_, v_, Lon_, Lat_, returnxy=True
        )
        quiver_color = kw.pop('quiver_color', 'g')
        quiver_scale = kw.pop('quiver_scale', None)
        hide_qkey = kw.pop('hide_qkey', False)
        qkey_kw = kw.pop('qkey_kw', {})
        qkey_X = kw.pop('qkey_X', 0.85)
        qkey_X = qkey_kw.pop('X', qkey_X)
        qkey_Y = kw.pop('qkey_Y', 1.02)
        qkey_Y = qkey_kw.pop('Y', qkey_Y)
        qkey_U = kw.pop('qkey_U', 2)
        qkey_U = qkey_kw.pop('U', qkey_U)
        qkey_label = kw.pop('qkey_label', '{:g} '.format(qkey_U) + units)
        qkey_label = qkey_kw.pop('label', qkey_label)
        qkey_labelpos = kw.pop('qkey_labelpos', 'W')
        qkey_labelpos = qkey_kw.pop('labelpos', qkey_labelpos)
        plot_obj = m.quiver(X_, Y_, u_rot, v_rot, color=quiver_color,
            scale=quiver_scale, **kw)
        if not hide_qkey:
            # quiverkey plot
            plt.quiverkey(plot_obj, qkey_X, qkey_Y, qkey_U,
                label=qkey_label, labelpos=qkey_labelpos, **qkey_kw)

    # hatch plot
    elif plot_type in ('hatch', 'hatches'):
        hatches = kw.pop('hatches', ['///'])
        plot_obj = m.contourf(X, Y, data, colors='none', hatches=hatches,
            extend='both', **kw)
    else:
        print('Please choose a right plot_type from ("pcolor", "contourf", "contour")!')

    # set clim
    if plot_type in ('pcolor', 'pcolormesh', 'imshow'):
        plt.clim(clim)

    # plot colorbar
    if plot_type in ('pcolor', 'pcolormesh', 'contourf', 'contourf+',
        'imshow'):
        ax_current = plt.gca()
        divider = make_axes_locatable(ax_current)
        cax = divider.append_axes(cbar_position, size=cbar_size, pad=cbar_pad)
        cbar = plt.colorbar(plot_obj, cax=cax, extend=cbar_extend,
            orientation=cbar_orientation, **cbar_kw)
        if cbar_type in ('v', 'vertical'):
            # put the units on the top of the vertical colorbar
            cbar.ax.xaxis.set_label_position('top')
            cbar.ax.set_xlabel(units)
            cbar.ax.set_ylabel(long_name)
        elif cbar_type in ('h', 'horizontal'):
            # cbar.ax.yaxis.set_label_position('right')
            # cbar.ax.set_ylabel(units, rotation=0, ha='left', va='center')
            if long_name == '' or units =='':
                cbar.ax.set_xlabel('{}{}'.format(long_name, units))
            else:
                cbar.ax.set_xlabel('{} [{}]'.format(long_name, units))
        # remove the colorbar to avoid repeated colorbars
        if hide_cbar:
            cbar.remove()
        # set back the main axes as the current axes
        plt.sca(ax_current)


    return plot_obj

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
fig=plt.figure(figsize=(9, 3))


map = Basemap(projection='sinu', 
              lat_0=0, lon_0=0)

lons = np.linspace(-180, 180, 10)
lats = np.linspace(-90, 90, 10)

lons, lats = np.meshgrid(lons, lats)

v10 = np.ones((lons.shape)) * 15
u10 = np.zeros((lons.shape))

u10_rot, v10_rot, x, y = map.rotate_vector(u10, v10, lons, lats, returnxy=True)

ax = fig.add_subplot(121)
ax.set_title('Without rotation')

map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
map.drawcoastlines(color = '0.15')


map.barbs(x, y, u10, v10, 
    pivot='middle', barbcolor='#333333')

ax = fig.add_subplot(122)
ax.set_title('Rotated vectors')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
def plot_flow_vectors_basemap(lons, lats, uu, vv, flow_speed, subregion:list=None, grid_shape=None, ax:Axes=None,
                              streamplot=False, draw_colorbar=True):
    from mpl_toolkits.basemap import Basemap

    if grid_shape is None:
        grid_shape = (300, 300)

    if subregion is None:
        subregion = [0, 1, 0, 1]

    if ax is None:
        fig = plt.figure()



    nx, ny = lons.shape





    b = Basemap(lon_0=180,
                llcrnrlon=lons[35, 35],
                llcrnrlat=lats[35, 35],
                urcrnrlon=lons[nx // 2, ny // 2],
                urcrnrlat=lats[nx // 2, ny // 2],
                resolution="i", area_thresh=2000)

    # im = b.pcolormesh(xx, yy, flow_speed)
    # b.colorbar(im)
    stride = 6



    uu1, vv1 = b.rotate_vector(uu, vv, lons, lats)

    lons_g, lats_g, xx_g, yy_g = b.makegrid(*grid_shape, returnxy=True)

    nx, ny = lons_g.shape
    i_start, i_end = int(nx * subregion[0]), int(nx * subregion[1])
    j_start, j_end = int(ny * subregion[2]), int(ny * subregion[3])



    xt, yt, zt = lat_lon.lon_lat_to_cartesian(lons_g.flatten(), lats_g.flatten())
    xs, ys, zs = lat_lon.lon_lat_to_cartesian(lons.flatten(), lats.flatten())
    ktree = KDTree(data=list(zip(xs, ys, zs)))

    dists, inds = ktree.query(list(zip(xt, yt, zt)))

    uu_to_plot = uu1.flatten()[inds].reshape(lons_g.shape)
    vv_to_plot = vv1.flatten()[inds].reshape(lons_g.shape)
    flow_speed_to_plot = flow_speed.flatten()[inds].reshape(lons_g.shape)



    clevs = [0, 0.01, 0.02, 0.04, 0.06, 0.08, 0.12, 0.16]
    ncolors = len(clevs) - 1
    norm = BoundaryNorm(clevs, ncolors)
    cmap = cm.get_cmap("gist_ncar_r", ncolors)

    im = b.pcolormesh(xx_g, yy_g, flow_speed_to_plot, alpha=0.5, cmap=cmap, norm=norm)
    if not streamplot:

        b.quiver(xx_g[i_start:i_end:stride, j_start:j_end:stride], yy_g[i_start:i_end:stride, j_start:j_end:stride],
                      uu_to_plot[i_start:i_end:stride, j_start:j_end:stride], vv_to_plot[i_start:i_end:stride, j_start:j_end:stride],
                      headlength=2, headaxislength=2, headwidth=4, units="inches", color="k")



    else:
        b.streamplot(xx_g, yy_g, uu_to_plot, vv_to_plot, linewidth=0.4, density=3, arrowstyle="fancy", arrowsize=0.4, ax=ax, color="k")


    # im = b.contourf(xx_g, yy_g, flow_speed_to_plot, levels=clevs, alpha=0.5, cmap=cmap, norm=norm)


    cb = b.colorbar(im, location="bottom")
    cb.ax.set_visible(draw_colorbar)


    b.drawcoastlines(linewidth=0.3, ax=ax)

    if ax is None:
        fig.savefig("nemo/circ_annual_mean_basemap.png", bbox_inches="tight", dpi=300)
        plt.close(fig)

    return im

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
  dates = [dates[tidx],]
  tidx_offset=tidx

idxs =arange(len(x))
np.random.shuffle(idxs)
idxs = idxs[:20000]

os.system('mkdir -p jpgs/vel')
if True:
  if True:
    fig=figure()
    fig.subplots_adjust(left=0.0,right=1.0,bottom=0.0,top=1.0)
    pc=tripcolor(x,y,nv,vel,cmap=cmap,rasterized=True)

    del hvel,zcor,uvel,vvel
    urot,vrot = proj.rotate_vector(u[idxs],v[idxs],lon[idxs],lat[idxs])
    del u,v
    #proj.quiver(x[idxs],y[idxs],u[idxs],v[idxs],pivot='mid',alpha=0.5)
    #urot=u[idxs]; vrot=v[idxs]
    if (uselim):
      upper_scale=vmax
    else:
      upper_scale=10.0 # deep transports
      upper_scale=0.1 # surface transports
    urot = ma.masked_where(vel[idxs]<0.05*upper_scale,urot)
    vrot = ma.masked_where(vel[idxs]<0.05*upper_scale,vrot)
    #quiverkws={'scale_units':'x','scale':10./10000.} # for deep transports
    quiverkws={'scale_units':'x','scale':upper_scale/10000.}
    proj.quiver(x[idxs],y[idxs],urot,vrot,pivot='mid',alpha=0.5,**quiverkws)

    # test vector

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
latg=nvargrd['lat_rho'][:]
long=nvargrd['lon_rho'][:]
nlat,nlon=long.shape


cang=nvargrd['tide_Cangle'][:]
cangr=pl.zeros_like(cang)


u0=pl.ones_like (long)
v0=pl.zeros_like(latg)

ur=pl.zeros_like (long)
vr=pl.zeros_like(latg)

ur,vr=proj.rotate_vector(u0,v0,long,latg)
alpha=pl.arccos(ur)*180./pi

for id in range(len(cang)):
    cangr[id,:]=cang[id,:]-alpha

cangr[cangr<0.]  =cangr[cangr<0.]  +360.
cangr[cangr>360.]=cangr[cangr>360.]-360.


nvargrd['tide_Cangle'][:]=cangr

#dimensions:
#    two = 2 ;
#    eta_rho = 300 ;
#    xi_rho = 240 ;

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
        num_days = size(files) #int(size(files)/int(num_days_lag)+1)
        print num_days
        if size(files)>0:
            drift_days = ma.masked_all((num_days,2,  nx, ny))
            drift_days_int_mask = ma.masked_all((num_days, 2, nx, ny))
            curl_days = ma.masked_all((num_days, nx, ny))
            for x in xrange(0, size(files)):
                print x
                day = '%02d' % x 
                f = Dataset(files[x], 'r')
                u = f.variables['zonal_motion'][0]/(60.*60.*24.*num_days_lag)
                v = f.variables['meridional_motion'][0]/(60.*60.*24.*num_days_lag)
                q = f.variables['quality_flag'][0]

                #ROTATE VECTORS TO X/Y GRID on NSIDC GRID FOR CURL CALC (SO IT IS FLAT ALONG THE BOTTOM)
                u_r,v_r = m.rotate_vector(u,v,lon,lat)
                #grid data onto 100km grid - also matches wind forcing fields
                drift_days[x, 0], drift_days[x, 1], q_int = BGF.interp_uvCSAT(u_r, v_r, q, xpts, ypts, xpts2m, ypts2m)
            
                curl_days[x] = BGF.calc_curl_sq_2d_xy_gradient(drift_days[x, 0], drift_days[x, 1], dx_res)
        drift_days_mask_count = ma.count_masked(drift_days[:, 0], axis=0)
        days_in_month=15
        curl_months[y, mon] = ma.masked_where(drift_days_mask_count>days_in_month, ma.mean(curl_days, axis=0))
        drift_month[y, mon, 0]  = ma.masked_where(drift_days_mask_count>days_in_month, ma.mean(drift_days[:, 0] , axis=0))
        drift_month[y, mon, 1]  = ma.masked_where(drift_days_mask_count>days_in_month, ma.mean(drift_days[:, 1] , axis=0))
  

curl_months.dump(outpath+str(start_year)+'-'+str(end_year)+'-curl_data_months'+grid_str+'.txt')


plot_curl=0

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import rotate_vector [as 别名]
vmin=np.min(speed); vmax=np.max(speed)

print(speed.max()) ; print(np.mean(speed)); print(np.std(speed))
#UN = uvel/speed ; VN = vvel/speed

#print fh.variables['hi'].units

#lon_0=-60.; lat_0=90.
m = Basemap(projection='npstere',boundinglat=55,lon_0=340,resolution='l', round=False)

xi, yi = m(lons, lats)

print(xi[0, :])
print(xi.min(), xi.max())

urot, vrot = m.rotate_vector(uvel, vvel, lons, lats)



#yy = np.arange(0, yi.shape[0], 10)
#xx = np.arange(0, xi.shape[1], 10)
#points = np.meshgrid(yy, xx)


plt.figure(figsize=(12,12))
ax = plt.gca()

#levs=np.arange(0.01,0.2,0.005)
cs = m.pcolor(xi,yi,urot,vmin=None, vmax=None)

stride = 10