Python Basemap.transform_scalar方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
 def mk_regmap_contour(a2in, bnd, scm, stitle, soname, cbarname, miss_out):
   #------------------------
   # Basemap
   #------------------------
   print "Basemap"
   figmap   = plt.figure()
   axmap    = figmap.add_axes([0.1, 0.0, 0.8, 1.0])
   M        = Basemap( resolution="l", llcrnrlat=lllat, llcrnrlon=lllon, urcrnrlat=urlat, urcrnrlon=urlon, ax=axmap)
   
   #-- transform -----------
   print "transform"
   a2value_trans    = M.transform_scalar( a2in,   a1lon, a1lat, nnx, nny)
   a2value_trans    = a2value_trans.filled(miss_out)
 
   #-- to contour ----------
   a2value_trans    = dtanl_fsub.mk_a2contour_regional(a2value_trans.T, 0.0, 0.0, -9999.0).T
 
   #-- boundaries ----------
   bnd_cbar   = [-1.0e+40] + bnd + [1.0e+40]
   
   #-- color ---------------
   cminst   = matplotlib.cm.get_cmap(scm, len(bnd))
   acm      = cminst( arange( len(bnd) ) )
   lcm      = [[1,1,1,1]]+ acm.tolist()
   mycm     = matplotlib.colors.ListedColormap( lcm )
   
   #-- imshow    -----------
   im       = M.imshow(a2value_trans, origin="lower", norm=BoundaryNormSymm(bnd), cmap=mycm, interpolation="nearest")
 
   #-- shade     -----------
   a2shade_trans = ma.masked_not_equal(a2value_trans, miss_out)
   cmshade       = matplotlib.colors.ListedColormap([(0.8,0.8,0.8), (0.8,0.8,0.8)])
   im            = M.imshow(a2shade_trans, origin="lower", cmap=cmshade) 
   
   #-- coastline ---------------
   print "coastlines"
   M.drawcoastlines()
   
   #-- meridians and parallels
   M.drawmeridians(arange(0.0,360.0, meridians), labels=[0, 0, 0, 1])
   M.drawparallels(arange(-90.0,90.0, parallels), labels=[1, 0, 0, 0])
   #-- title -------------------
   axmap.set_title("%s"%(stitle))
   #-- save --------------------
   plt.savefig(soname)
   print soname 
   # for colorbar ---
   if cbarflag == "True":
     figmap   = plt.figure()
     axmap    = figmap.add_axes([0.1, 0.0, 0.8, 1.0])
     M        = Basemap( resolution="l", llcrnrlat=lllat, llcrnrlon=lllon, urcrnrlat=urlat, urcrnrlon=urlon, ax=axmap)
     a2v_trans    = M.transform_scalar( a2in,   a1lon, a1lat, nnx, nny)
     im       = M.imshow(a2value_trans, origin="lower", norm=BoundaryNormSymm(bnd), cmap=mycm)
   
     figcbar   = plt.figure(figsize=(5, 0.6))
     axcbar    = figcbar.add_axes([0, 0.4, 1.0, 0.6])
     bnd_cbar  = [-1.0e+40] + bnd + [1.0e+40]
     plt.colorbar(im, boundaries= bnd_cbar, extend="both", cax=axcbar, orientation="horizontal")
     figcbar.savefig(cbarname)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
 def background_map(self, ax):
     llcrnrlat, urcrnrlat, llcrnrlon, urcrnrlon, lat_ts = (31, 44, -126, -113, 37.5)
     m = Basemap(projection='merc', llcrnrlat=llcrnrlat,
                 urcrnrlat=urcrnrlat, llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon,
                 lat_ts=lat_ts, resolution='i', ax=ax)
     m.drawmapboundary(fill_color='lightblue', zorder=0)
     m.fillcontinents(zorder=0)
     etopofn = '/home/behry/uni/data/etopo1_central_europe_gmt.grd'
     etopodata = Dataset(etopofn, 'r')
     z = etopodata.variables['z'][:]
     x_range = etopodata.variables['x_range'][:]
     y_range = etopodata.variables['y_range'][:]
     spc = etopodata.variables['spacing'][:]
     lats = np.arange(y_range[0], y_range[1], spc[1])
     lons = np.arange(x_range[0], x_range[1], spc[0])
     topoin = z.reshape(lats.size, lons.size, order='C')
     # transform to nx x ny regularly spaced 5km native projection grid
     nx = int((m.xmax - m.xmin) / 5000.) + 1; ny = int((m.ymax - m.ymin) / 5000.) + 1
     topodat, x, y = m.transform_scalar(np.flipud(topoin), lons, lats, nx, ny, returnxy=True)
     ls = LightSource(azdeg=300, altdeg=15, hsv_min_sat=0.2, hsv_max_sat=0.3,
                      hsv_min_val=0.2, hsv_max_val=0.3)
     # shade data, creating an rgb array.
     rgb = ls.shade(np.ma.masked_less(topodat / 1000.0, 0.0), cm.gist_gray_r)
     m.imshow(rgb)
     m.drawmeridians(np.arange(6, 12, 2), labels=[0, 0, 0, 1], color='white',
                     linewidth=0.5, zorder=0)
     m.drawparallels(np.arange(44, 50, 2), labels=[1, 0, 0, 0], color='white',
                     linewidth=0.5, zorder=0)
     m.drawcoastlines(zorder=1)
     m.drawcountries(linewidth=1.5, zorder=1)
     m.drawstates()
     m.drawrivers(color='lightblue', zorder=1)
     return m

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
def plot_world_sst():
    # Read some NetCDF data
    import netCDF4 as nc
    
    ostia = nc.Dataset('ostia.nc')
    tmp = ostia.variables['analysed_sst'][0]
    ice = ostia.variables['sea_ice_fraction'][0]
    lon = ostia.variables['lon'][:]
    lat = ostia.variables['lat'][:]

    from mpl_toolkits.basemap import Basemap

    # Set up a map
    map = Basemap(projection='cyl')
    map.drawcoastlines()
    map.drawcountries()
    map.fillcontinents(color='lightgreen', lake_color='lightblue');
    map.drawmapboundary(fill_color='lightblue')

    # Re-project the data onto the map
    image = map.transform_scalar(tmp,lon,lat,200,200)

    # Plot the data
    map.imshow(image);

    plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
def stere(data,latitude,longitude,filename):
    """
    :param data: CO浓度数据矩阵
    :param latitude: 纬度
    :param longitude: 经度
    :param filename:
    :return:
    """
    ax = plt.gca()
    fig = plt.figure()
    m = Basemap(width=10000000,height=7500000,
            resolution='l',projection='stere',\
            lat_ts=35,lat_0=35,lon_0=107.)
    nx = int((m.xmax-m.xmin)/5000.)+1
    ny = int((m.ymax-m.ymin)/5000.)+1
    topodat = m.transform_scalar(data,longitude,latitude,nx,ny)
    im = m.imshow(topodat,cm.GMT_haxby,vmin=0,vmax=4e18)
    m.drawcoastlines()
    m.drawcountries()
    m.drawparallels(np.arange(-80.,81.,20.),labels=[1,0,0,0])
    # 画平行的纬度,前一个参数是表示起点，终点，间距的序列，后一个参数是指在哪一个方向显示纬度、经度值
    m.drawmeridians(np.arange(-180.,181.,20.),labels=[0,0,0,1])
    # 画平行的经度
    m.colorbar(im)
    plt.title("CO"+filename[12:14])
    outname = filename+'.png'
    fig.savefig(outname, dpi=fig.dpi)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
def plot_data(data,lon_data, lat_data, periodname, AODcatname,maptype,cmapname,minv=0,maxv=0,folder=""):
    fig = plt.figure()
    #ax = fig.add_axes([0.1,0.1,0.8,0.8])

    m = Basemap(llcrnrlon=19,llcrnrlat=34,urcrnrlon=29,urcrnrlat=42,
                resolution='h',projection='cass',lon_0=24,lat_0=38)

    nx = int((m.xmax-m.xmin)/1000.)+1
    ny = int((m.ymax-m.ymin)/1000.)+1
    topodat = m.transform_scalar(data,lon_data,lat_data,nx,ny)
    
    if minv<>0 or maxv<>0 :
        im = m.imshow(topodat,cmap=plt.get_cmap(cmapname),vmin=minv,vmax=maxv)
    else:
        im = m.imshow(topodat,cmap=plt.get_cmap(cmapname))
 

    m.drawcoastlines()
    m.drawmapboundary()
    m.drawcountries()
    m.drawparallels(np.arange(35,42.,1.), labels=[1,0,0,1])
    m.drawmeridians(np.arange(-20.,29.,1.), labels=[1,0,0,1])
    cb = m.colorbar(im,"right", size="5%", pad='2%')
    title=maptype+" AOD "+AODcatname+" "+periodname+" 2007-2014"
    plt.title(title)
    pylab.savefig(folder+maptype+"AOD"+AODcatname+"_"+periodname + ".png")

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
    def map_slab1(self):
        '''
        Creates a map of the slab1.0 contours
        '''

        if self.slab1_details:
            bbox = self.slab1_details['Bounding_box']
            print(bbox)

            minlon = bbox[0]
            maxlon = bbox[1]
            minlat = bbox[2]
            maxlat = bbox[3]

            print('------------------------\nPlotting map...\n------------------------')

            figure = plt.figure(facecolor='white',figsize=(10,8))
            a = figure.add_subplot(111)
            m = Basemap(ax=a,lat_0=(minlat+((maxlat-minlat)/2.0)),lon_0=(minlon+((maxlon-minlon)/2.0)),resolution ='l',llcrnrlon=minlon,llcrnrlat=minlat,urcrnrlon=maxlon,urcrnrlat=maxlat)

            #Create grid uniformly spaced grid
            nx = int((m.xmax-m.xmin)/0.1)+1; ny = int((m.ymax-m.ymin)/0.1)+1

            slabdat = m.transform_scalar(self.slab1_details['Depth_array'],self.slab1_details['Lon_array'],self.slab1_details['Lat_array'],nx,ny)
            im = m.imshow(slabdat)

            m.drawparallels(np.arange(minlat,maxlat,((maxlat-minlat)/5)),labels=[1,1,0,0],linewidth=0.5,fontsize=10)
            m.drawmeridians(np.arange(minlon,maxlon,((maxlon-minlon)/5)),labels=[0,0,0,1],linewidth=0.5,fontsize=10)

            m.drawcoastlines()

            #Add tectonic plate boundaries - this doesn't work very well
            # self.faults = pts.read_faults('data/plates.xy')
            # for i in self.faults:
            #     faults_lons = self.faults[i][0]
            #     faults_lats = self.faults[i][1]
            #     x,y = m(faults_lons, faults_lats)
            #     m.plot(x,y,'b--',linewidth=1.0)

            #add a colorbar
            cb = m.colorbar(im,"bottom", size="5%", pad='15%')
            cb.set_label('Depth to slab surface [km]')

            #add a figure title
            a.set_title('Slab1.0 upper surface depths: %s' %self.name)

            return a,m

        else:
            print('Nothing to do: no slab1_details added!')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
def render_components(C, lats, lons, fname_tmpl = None, ndx = None):
    """
    Render the components in C [with dims comp_id x lats x lons] onto
    a world map of appropriate size.  The method assumes that the  lons x lats
    generate a rectangle.
    """
    m = Basemap(projection='mill',llcrnrlat=min(lats), urcrnrlat=max(lats), llcrnrlon=(min(lons)),
                urcrnrlon=max(lons),resolution='c')
    
    if ndx == None:
        ndx = np.arange(len(C)) + 1
        
    rmax = np.max(C)

    # lattitudes may not be sorted in ascending order, rectify this
    lat_ndx = np.argsort(lats)
    lats_s = lats[lat_ndx]

    # render each topo plot
    for ci in range(len(C)):
        
        plt.figure(figsize=(12, 8 * (max(lats) - min(lats)) / 180))
        plt.axes([0.05, 0.05, 0.9, 0.85])

        m.drawcoastlines()
        #m.fillcontinents(color='coral',lake_color='aqua', zorder = 0)
        m.etopo(scale = 0.2)

        # draw parallels and meridians.
        m.drawparallels(np.arange(-90.,91.,30.))
        m.drawmeridians(np.arange(0.,361.,60.))
        #m.drawmapboundary(fill_color='aqua')
        
        nx = int((m.xmax-m.xmin) / 20000) + 1
        ny = int((m.ymax-m.ymin) / 20000) + 1
        Ci = C[ci, lat_ndx, :]
        f = m.transform_scalar(Ci, lons, lats_s, nx, ny)
        
        # imlim seems to be equivalent to caxis() im MATLAB
        imgplt = m.imshow(f, alpha = 0.8)
        imgplt.set_clim(-rmax, rmax)
        plt.colorbar()
        
        plt.title('Component %d' % (ndx[ci]))
        
        if fname_tmpl:
            plt.savefig(fname_tmpl % (ndx[ci]))
            
    if not fname_tmpl:
        plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
def Mapshow(dims, data, type, para1, para2, tit, unit):

    # Prepare for drawing
    lons = np.arange(dims["minlon"], dims["maxlon"] + dims["res"] / 2, dims["res"])
    lats = np.arange(dims["minlat"], dims["maxlat"] + dims["res"] / 2, dims["res"])
    x, y = np.meshgrid(lons, lats)
    # draw Chile Basemap with lambert projection at normal x, y settings
    m = Basemap(
        llcrnrlon=dims["minlon"],
        llcrnrlat=dims["minlat"],
        urcrnrlon=dims["maxlon"],
        urcrnrlat=dims["maxlat"],
        projection="cyl",
        fix_aspect=True,
        lat_1=-10,
        lat_2=10,
        lon_0=20,
    )  # projection='lcc'
    # draw boundaries
    m.drawcoastlines()
    m.drawcountries(linewidth=2)
    m.drawstates()
    m.drawparallels(arange(-20, 30, 20), labels=[1, 0, 0, 0])  # only left ytick
    m.drawmeridians(arange(-10, 60, 20), labels=[0, 0, 0, 1])  # only bottom xtick
    # for  the classified figure
    X, Y = m(x, y)
    if type == "imshow":
        # im = m.contourf(X, Y, data, cmap=plt.cm.bwr, extend='both')
        plotdata = m.transform_scalar(data, lons, lats, dims["nlon"], dims["nlat"])
        im = m.imshow(plotdata, vmin=para2, vmax=para1, cmap=plt.cm.bwr)
    elif type == "contour":
        if para1 is not None:
            im = m.contourf(X, Y, data, para1, cmap=plt.cm.bwr, extend="both")
        else:
            im = m.contourf(X, Y, data, cmap=plt.cm.bwr, extend="both")
        cb = m.colorbar(im, pad="3%", ticks=para2)  # cb = m.colorbar(im, location='bottom', pad='16%')
        # The following method will return the zero as a small number close to zero
        # if para2 is not None:
        # 	cb.set_ticks(para2)
        # 	cb.set_ticklabels(para2)

        # map data with lon and lat position
    plt.title(tit, fontsize=20)
    plt.xlabel(unit, fontsize=18, labelpad=15)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
def millshow(zz,lon,lat,lonmin=-180.,lonmax=180.,latmin=-70,latmax=90, \
            theTitle=None,cmap=None,clim=None,map_type='coastline', \
            cbar='v',cbartitle=None,units=None):
    """ mill projection, imshow plot 
    """
    if cmap is None:
        cmap = plt.cm.jet
    elif cmap.startswith('GMT',0,3):
        cmap = eval('cm.' + cmap) # use cmap from Basemap module cm
    if clim is None:
        if isinstance(zz,np.ma.core.MaskedArray):
            zz1d = zz.compressed()
        else:
            zz1d = zz.ravel()
        clim = np.percentile(zz1d,1), np.percentile(zz1d,99)
    
    # base map
    lon_0 = (lon[0] + lon[-1])/2.
    m = Basemap(projection='mill', lon_0=lon_0, llcrnrlon=lonmin, \
                llcrnrlat=latmin, urcrnrlon=lonmax, urcrnrlat=latmax)
    if map_type=='coastline':
        m.drawcoastlines(linewidth=0.25)
    elif map_type=='filled':
        m.fillcontinents(color=[.33, .33, .33])
    
    # show data
    nx = len(lon); ny = len(lat)
    data = m.transform_scalar(zz,lon,lat,nx,ny)
    im = m.imshow(data,cmap)
    
    # labels        
    plt.clim(clim)
    m.drawmapboundary(fill_color='grey')
    m.drawparallels(np.arange(-90.,120.,30.),labels=[1,0,0,0])
    m.drawmeridians(np.arange(0.,420.,60.),labels=[0,0,0,1])
    if theTitle is not None:
        plt.title(theTitle)
    if cbar=='v':
        vcolorbar(im=im,cbartitle=cbartitle,units=units)
    elif cbar=='h':
        hcolorbar(im=im,cbartitle=cbartitle,units=units)
    return im    

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

#.........这里部分代码省略.........
#        dir_name = os.path.dirname(in_file)
#        if not os.path.exists(dir_name):
#            os.makedirs(dir_name)

        vgFile = open(in_file+'.vg'+str(vg_files[plot_ctr]))
        pattern = re.compile(r"[a-z]+")

        for line in vgFile:
            match = pattern.search( line )
            if not match:
                value = float(line[int(self.image_defs['first_char']):int(self.image_defs['last_char'])]) 
                # TODO Does this need to be normalised here if it's also being done in the plot?
                value = max(self.image_defs['min'], value)
                value = min(self.image_defs['max'], value)
                #if value < self.image_defs[2] : value = self.image_defs[2] 
                #if value > self.image_defs[3] : value = self.image_defs[3] 
                points[int(line[3:6])-1][int(line[0:3])-1] = value
        vgFile.close()

        map = Basemap(\
            llcrnrlon=area_rect.get_sw_lon(), llcrnrlat=area_rect.get_sw_lat(),\
            urcrnrlon=area_rect.get_ne_lon(), urcrnrlat=area_rect.get_ne_lat(),\
            projection=projection,\
            lat_0=plot_centre_location.get_latitude(),\
            lon_0=plot_centre_location.get_longitude(),\
            resolution=resolution,
            ax=ax)
            
        map.drawcoastlines(color='black')
        map.drawcountries(color='grey')
        map.drawmapboundary(color='black', linewidth=1.0)

        warped = ma.zeros((grid, grid),float)
        warped, warped_lon, warped_lat = map.transform_scalar(points,lons,lats,grid,grid, returnxy=True, checkbounds=False, masked=True)
        warped = warped.filled(self.image_defs['min']-1.0)
        
        colMap.set_under(color ='k', alpha=0.0)
    
        im = map.imshow(warped,
            cmap=colMap,
            extent = (-180, 180, -90, 90),
            origin = 'lower',
            norm = P.Normalize(clip = False,
            vmin=self.image_defs['min'],
            vmax=self.image_defs['max']))


        #######################
        # Plot greyline
        #######################
        if plot_terminator:
            the_sun = Sun()
            the_month = plot_parameters.get_month(vg_files[plot_ctr]-1)
            the_day = plot_parameters.get_day(vg_files[plot_ctr]-1)
            the_hour = plot_parameters.get_utc(vg_files[plot_ctr]-1)
            if (the_day == 0):
                the_day = 15
            the_year = datetime.date.today().year
            num_days_since_2k = the_sun.daysSince2000Jan0(the_year, the_month, the_day)

            res =  the_sun.sunRADec(num_days_since_2k)
            declination = res[1]
            if(declination==0.0):
                declination=-0.001

            tau = the_sun.computeGHA(the_day, the_month, the_year, the_hour);

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
lons = np.linspace(20,380,num=hsize);#etopodata.variables['ETOPO05_X'][:]
lats = np.linspace(-90,90,num=vsize);#etopodata.variables['ETOPO05_Y'][:]
topoin,lons = shiftgrid(180.,topoin,lons,start=False)
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])

# ('lcc' = lambert conformal conic).
# use major and minor sphere radii from WGS84 ellipsoid.
m = Basemap(llcrnrlon=-145.5,llcrnrlat=1.,urcrnrlon=-2.566,urcrnrlat=46.352,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='l',area_thresh=1000.,projection='lcc',\
            lat_1=50.,lon_0=-107.,ax=ax)

# transform to nx x ny regularly spaced 5km native projection grid
nx = int((m.xmax-m.xmin)/5000.)+1; ny = int((m.ymax-m.ymin)/5000.)+1
topodat = m.transform_scalar(topoin,lons,lats,nx,ny)

# plot image over map with imshow.
im = m.imshow(topodat,cm.GMT_haxby)

# draw coastlines and boundaries.
m.drawcoastlines()
m.drawcountries()
m.drawstates()
# draw parallels and meridians.
# label on left and bottom of map.
parallels = np.arange(0.,80,20.)
m.drawparallels(parallels,labels=[1,0,0,1])
meridians = np.arange(10.,360.,30.)
m.drawmeridians(meridians,labels=[1,0,0,1])
# add colorbar

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
xres = 360.0 / nx
lats = linspace(-90.0, 90.0, ny)
lons = linspace(0.0, 360.0 - 360.0/nx, nx)
lllat  = 20.0
lllon  = 120.0
urlat  = 50.0
urlon  = 160.0

nnx    = int( (urlon - lllon)/ xres)
nny    = int( (urlat - lllat)/ yres)
#-----------------------
dir_his = "/media/disk2/out/CMIP5/day/NorESM1-M/historical/r1i1p1/tracks/dura24/wfpr/mp"
shis   = dir_his + "/mp.p00.00.c01.04.r1000.nw17_DJF_day_NorESM1-M_historical_r1i1p1.bn"
adatin = fromfile(shis, float32).reshape(17,96,144)[0]
#
figmap = plt.figure()
axmap  = figmap.add_axes([0.1, 0.1, 0.8, 0.8])
M      = Basemap(resolution="l", llcrnrlat=lllat, llcrnrlon=lllon, urcrnrlat=urlat, urcrnrlon=urlon, ax=axmap)
#--

adat   = M.transform_scalar(adatin, lons, lats, nnx, nny)



#--
im     = M.imshow(adat, origin="lower", cmap="RdBu", interpolation="nearest")
M.drawcoastlines()
#
plt.savefig("a2.png")
plt.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
if POPDIAGPY == 'TRUE':
  outfile = os.path.join(outdir,'map_nutlim_sp.png')
else:
  outfile = os.path.join(outdir,'map_%04d-%04d_nutlim_sp.png'%(yrstart,yrend))

print('plotting %s'%(outfile))
fig = pl.figure()
ax  = fig.add_subplot(111)
map = Basemap(llcrnrlon=30,urcrnrlon=360+30,llcrnrlat=-85,urcrnrlat=85,
        resolution='l',projection='cyl',ax=ax)
map.drawcoastlines(linewidth=0.5)
map.drawrivers(linewidth=0.75,color='white')
map.drawparallels(N.arange(-90,120,30),labels=[0,0,0,0],linewidth=0.25)
map.drawmeridians(N.arange(-720,750,30),labels=[0,0,0,0],linewidth=0.25)
map.fillcontinents('white')
data, x, y = map.transform_scalar(sp_lim,tlon,tlat,nlon,nlat,returnxy=True)
IM  = map.imshow(data,norm=colors.normalize(vmin=0,vmax=6.5),cmap=cmap)
CB  = fig.colorbar(IM,ax=ax,shrink=0.5,aspect=12,drawedges=True,ticks=cbticks,
        boundaries=cblevels,values=cbticks)
CB.ax.set_yticklabels(cbticklabels)
for line in CB.ax.get_yticklines(): # remove tick lines
    line.set_markersize(0)

bottom_labels(ax,case.replace('_',' '),'',year_string,fontsize=10)
ax.set(title='Small Phyto. growth limitation factor in upper %.0f m'%(zref))
fig.savefig(outfile,dpi=150,bbox_inches='tight',pad_inches=0.2)
pl.close(fig)

# Diatom nutrient limitation plot
if POPDIAGPY == 'TRUE':
  outfile = os.path.join(outdir,'map_nutlim_diat.png')

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
nnx        = int( (urlon - lllon)/dlon)
nny        = int( (urlat - lllat)/dlat)
a1lon_loc  = linspace(lllon, urlon, nnx)
a1lat_loc  = linspace(lllat, urlat, nny)
LONS, LATS = meshgrid( a1lon_loc, a1lat_loc)
#------------------------
# Basemap
#------------------------
print "Basemap"
figmap   = plt.figure()
axmap    = figmap.add_axes([0.1, 0.0, 0.8, 1.0])
M        = Basemap( resolution="l", llcrnrlat=lllat, llcrnrlon=lllon, urcrnrlat=urlat, urcrnrlon=urlon, ax=axmap)

#-- transform -----------
print "transform"
a2u_trans    = M.transform_scalar( a2u,   a1lon, a1lat, nnx, nny) 
a2v_trans    = M.transform_scalar( a2v,   a1lon, a1lat, nnx, nny) 
#--- 
a2psl_trans  = M.transform_scalar( a2psl, a1lon, a1lat, nnx, nny)
#

#-- boundaries ----------
#bnd        = [1,3,5,7,9,11,13,15,17]
bnd        = list(arange(0,21+1,3))
bnd_cbar   = [-1.0e+40] + bnd + [1.0e+40]

#-- color ---------------
#scm      = "rainbow"
scm      = "Spectral_r"
#scm      = "gist_rainbow"
cminst   = matplotlib.cm.get_cmap(scm, len(bnd))

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import transform_scalar [as 别名]
# longitudes go from 20 to 380.
topoin = np.loadtxt('etopo20data.gz')
lons = np.loadtxt('etopo20lons.gz')
lats = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topoin,lons,start=False)

# setup of basemap ('lcc' = lambert conformal conic).
# use major and minor sphere radii from WGS84 ellipsoid.
m = Basemap(llcrnrlon=-145.5,llcrnrlat=1.,urcrnrlon=-2.566,urcrnrlat=46.352,\
            rsphere=(6378137.00,6356752.3142),\
            resolution='l',area_thresh=1000.,projection='lcc',\
            lat_1=50.,lon_0=-107.)
# transform to nx x ny regularly spaced native projection grid
nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
topodat,x,y = m.transform_scalar(topoin,lons,lats,nx,ny,returnxy=True)
# create the figure.
fig=plt.figure(figsize=(8,8))
# add an axes, leaving room for colorbar on the right.
ax = fig.add_axes([0.1,0.1,0.7,0.7])
# associate this axes with the Basemap instance.
m.ax = ax
# plot image over map with imshow.
im = m.imshow(topodat,plt.cm.jet)
# setup colorbar axes instance.
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.075, b, 0.05, h])
plt.colorbar(im,cax=cax) # draw colorbar
# plot blue dot on boulder, colorado and label it as such.
xpt,ypt = m(-104.237,40.125)