Python basemap.Basemap类代码示例

    def _draw_basemap(self, proj_prop=None, drawparallels=True, vmin_polygons=None, vmax_polygons=None, **kwargs):
        """
        """
        if proj_prop is None:
            raise ValueError(
                'No projection properties are given! Please modify or choose a different backend!')

        the_map = Basemap(ax=self.pax, **proj_prop)
        xm = self.x.timmean()

        Z = xm
        lon = self.x.lon
        lat = self.x.lat

        X, Y = the_map(lon, lat)
        self.im = the_map.pcolormesh(X, Y, Z, **kwargs)

        self.__basemap_ancillary(the_map, drawparallels=drawparallels)

        # add polygons to map
        if self.polygons is not None:
            if False:  # individual polygons
                for p in self.polygons:
                    self._add_single_polygon_basemap(the_map, p)
            else:  # plot all polygons at once
                self._add_polygons_as_collection_basemap(
                    the_map, vmin=vmin_polygons, vmax=vmax_polygons)

    def mapMake(self,renderdpi,table,msize):
        cn = 0
        start = time()

        #initialize connection to database
        cn = psycopg2.connect(secret.DB_CONNECT)
        cr = cn.cursor()

        #get map ready to go
        fig = plt.figure(figsize=(8,4),dpi = renderdpi)
        fig.add_subplot(1,1,1)
        m = Basemap(projection='merc',llcrnrlat=-60,urcrnrlat=75,\
            llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='i')
        m.drawcoastlines(linewidth=.05)
        m.drawcountries(linewidth=.05)

        photoCnt = 0
        points = []

        cr.execute('SELECT latitude,longitude FROM %s;' %  table)
        for row in cr.fetchall():
            x,y = m(row[1],row[0])#convert to merc projection coords
            points.append((x,y))
            photoCnt += 1
        xs,ys = zip(*points)

        plt.title("%d %s" % (photoCnt,table))
        plt.scatter(xs,ys,s=msize,marker='.',c='green',edgecolors='none')
        plt.savefig(table,dpi = renderdpi)

        print "{:d} {:s}  mapped in {:f} seconds".format(photoCnt,\
                                                         table,time()-start)

def draw_basic_map_of_toronto(axis):
    """Draw a basic map of Toronto.

    :param axis: where to draw the map
    :returns: Basemap with Toronto
    """
    # These are the latitudes of the City of Toronto.
    # (Fiona can be better to find them, e.g., from the 'icitw_wgs84' Shapefile
    #  below.)

    low_left_corner_longitude = -79.75
    low_left_corner_latitude = 43.40
    up_right_corner_longitude = -79.10
    up_right_corner_latitude = 43.95

    to_map = Basemap(
        llcrnrlon=low_left_corner_longitude,
        llcrnrlat=low_left_corner_latitude,
        urcrnrlon=up_right_corner_longitude,
        urcrnrlat=up_right_corner_latitude,
        ellps="WGS84",
        resolution="h",
        area_thresh=0.1,
        ax=axis,
    )

    to_map.drawmapboundary(fill_color="white")
    return to_map

def basemap_raster_mercator(lon, lat, grid, cmap = None):
  """
  Render a raster in mercator projection.  Locations with no values are
  rendered transparent.
  """
  # longitude/latitude extent
  lons = (np.amin(lon), np.amax(lon))
  lats = (np.amin(lat), np.amax(lat))

  if cmap is None:
    cmap = mpl.cm.jet
    cmap.set_bad('w', 1.0)

  # construct spherical mercator projection for region of interest
  m = Basemap(projection='merc',llcrnrlat=lats[0], urcrnrlat=lats[1],
              llcrnrlon=lons[0],urcrnrlon=lons[1])

  vmin,vmax = np.nanmin(grid),np.nanmax(grid)
  masked_grid = np.ma.array(grid,mask=np.isnan(grid))
  fig = plt.figure(frameon=False)
  plt.axis('off')
  m.pcolormesh(lon,lat,masked_grid,latlon=True,cmap=cmap,vmin=vmin,vmax=vmax)

  str_io = StringIO.StringIO()
  plt.savefig(str_io,bbox_inches='tight',format='png',pad_inches=0,transparent=True)
  bounds = [ (lons[0],lats[0]),(lons[1],lats[0]),(lons[1],lats[1]),(lons[0],lats[1]) ]

  return str_io.getvalue(), bounds

def plot_drainage_areas(path = "data/hydrosheds/test_using_splitting_amno.nc"):
    ds = Dataset(path)
    #basemap = polar_stereographic.basemap
    basemap =  Basemap()
    lons = ds.variables["lon"][:]
    lats = ds.variables["lat"][:]
    channel_slope = ds.variables["slope"][:]

    lons[lons < 0] += 360

    x, y = basemap(lons, lats)




    acc_area = ds.variables["accumulation_area"][:]
    acc_area = np.log(acc_area)
    acc_area = np.ma.masked_where(channel_slope < 0, acc_area)

    basemap.pcolormesh(x, y, acc_area)
    basemap.drawcoastlines()
    plt.colorbar()
    plt.xlim(x.min(), x.max())
    plt.ylim(y.min(), y.max())
    plt.show()

    def draw_samples(self, n_samples, burnin=100, thin=10):
        if not self.kde_model:
            self.fit_kde()

        # pick starting point
        start_x = np.array([[0., 31.]])

        # get region for restricting samples
        margin = 0.5
        m = Basemap(llcrnrlon=self.uganda_data.LONGITUDE.min() - margin,
            llcrnrlat=self.uganda_data.LATITUDE.min() - margin,
            urcrnrlon=self.uganda_data.LONGITUDE.max() + margin,
            urcrnrlat=self.uganda_data.LATITUDE.max() + margin,
            resolution='l',
            area_thresh=10000)
        m.readshapefile("data/regions/UGA_adm0", "regions", drawbounds=True)
        for xy, info in zip(m.regions, m.regions):
            p = path.Path(xy)

        slice_samples, _ = SliceSampler.SliceSampler.mvslice(self.kde_model.pdf, 
            start_x,
            sample_size=n_samples, 
            burnin=burnin, 
            thin=thin, 
            confirm_region=p)
        
        # we get lat, long; we want x, y
        return np.fliplr(slice_samples)

    def plot2D_parameter(self,grid,filename):
      
        fig = plt.figure(figsize=(figsize,figsize),dpi=self.cfg.map_dpi)
        
        Lon,Lat = np.meshgrid(self.lon, self.lat)
       
        map = Basemap(llcrnrlon=self.BB['lon'][0],llcrnrlat=self.BB['lat'][0],urcrnrlon=self.BB['lon'][1],urcrnrlat=self.BB['lat'][1],
                      rsphere=(6378137.00,6356752.3142), resolution='l',projection='merc')
    
        x, y = map(Lon, Lat)
        cs2 = map.contourf(x,y,grid,512)
        
        plt.subplots_adjust(left=0.0, right=figsize, top=figsize, bottom=0.0)
        
        self.saveplot(plt,filename)
        plt.close()

        fig1 = plt.figure(figsize=(11,11),dpi=100)
        cb = map.colorbar(cs2,"bottom", size="5%", pad="100%")
        
        name = filename + '_legend.jpg' 
        plt.savefig("maps/"+name, dpi=100,bbox_inches='tight',transparent=False)
        plt.close()
        list_of_maps_to_send.append("maps/"+name)

        img1 = Image.open("maps/"+name)
        w, h = img1.size
        box = (20, h-70, w-10, h-10)
        area = img1.crop(box)
        area.save("maps/"+name,'jpeg')
        plt.close()

def nice_plot(data,xmin,xmax,xint,centerlat,centerlon,stations,color,cmin,
              cmax,levels_t):
    """Make plots in map projection, requires input array, x-min max and 
       interval (also used for y), center of data in lat, lon, station 
       locations, color scale, min and max limits and levels array for color
       scale.
    """
    domain = xmax-xint/2.
    maps = Basemap(projection='laea',lat_0=centerlat,lon_0=centerlon,
                   width=domain*2,height=domain*2)
    s = plt.pcolormesh(np.arange(xmin-xint/2.,xmax+3*xint/2.,xint)+domain,
                       np.arange(xmin-xint/2.,xmax+3*xint/2.,xint)+domain,
                       data, cmap = color)
    s.set_clim(vmin=cmin,vmax=cmax)
    CS = plt.contour(np.arange(xmin,xmax+xint,xint)+domain,
                     np.arange(xmin,xmax+xint,xint)+domain,
                     data, colors='k',levels=levels_t)
    plt.clabel(CS, inline=1, fmt='%1.2f',fontsize=8)
    plt.scatter(stations[:,0]+domain, stations[:,1]+domain, color='k',s=2)
    maps.drawstates()
    fig = plt.gcf()
    circle=plt.Circle((domain,domain),100000,color='0.5',fill=False)
    fig.gca().add_artist(circle)
    circle=plt.Circle((domain,domain),200000,color='0.5',fill=False)
    fig.gca().add_artist(circle)

def arc_map(region=[-33.7, 149.4, -31.8, 151.8], res=1000, service='World_Topo_Map',epsg=4326,
            resolution='h',projection='mill'):
    '''
    Create and return the thingy
    Inputs:
        region:  [bot,left,top,right]
        service= one of 'World_Imagery', or 'World_Topo_Map'
            Imagery is more like blue marble, with no site labels
            Topo has inbuilt site labels, which render based on the resolution -> higher res for more zoomed in image.
        epsg: tell arcGIS where to pull image from, affects the image orientation
            Use 4326 for Sydney centric maps
            Use 3577 for Australia centred view
    '''

    # set up basemap projection
    # epsg is used to inform basemap that we will pull an image from ArcGIS servers
    # epsg=28355 does not work...
    # 3577 is centred at middle of Australia...
    # 4326 works for East Australia!
    m = Basemap(projection=projection, epsg=epsg,
                resolution=resolution,
                llcrnrlon=region[1], llcrnrlat=region[0],
                urcrnrlon=region[3], urcrnrlat=region[2])

    # Download backgroud to basemap's projection.
    #m.arcgisimage(service='World_Topo_Map', xpixels=res, verbose=True)
    m.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service=service,
                  xpixels=res, verbose=True)

    return m

def get_basemap():
# Lambert Conformal map of lower 48 states
    m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
    shp_info = m.readshapefile('maps/st99_d00','states',drawbounds=True)
#   print(shp_info)
    return m

    def drawMap(self, var, strMonth):                

        m = Basemap( projection = 'cyl',llcrnrlon = p.lon0, \
                     llcrnrlat = self.lat0, urcrnrlon = self.lon0 + self.dxy * self.nX,\
                 urcrnrlat = self.lat0 + self.dxy * self.nY, resolution='i')

        var = np.flipud( np.ma.masked_array( var, var >= -0.2) )

        var = np.ma.masked_array( var, var < -10)

        levels = [ -0.1, -0.3, -0.5, -0.7, -0.9]

        im = m.pcolormesh( self.longitude, self.latitude,var, vmin = -0.9, vmax = 0.1)

        m.drawmeridians(np.arange( -10.0, 30, 5), labels = [1,0,0,1], fontsize = 10, rotation = 20)
        m.drawparallels(np.arange( 30.0, 60.0, 5), labels = [1,0,0,1], fontsize = 10, rotation = 20)
        
        m.drawcoastlines()
        m.drawcountries()
        plt.colorbar(im, orientation = 'vertical' , pad = 0.05 )


        plt.title( "CorrCoeff for month " + strMonth)


        plt.savefig( "corrCoeff" + strMonth + ".png", orientation = "landscape")
        plt.clf()

def ray_density(lat1, lon1, lat2, lon2,
                dt=1, gr_x=360, gr_y=180, 
                npts=180, projection='robin', 
                ray_coverage=False):
    '''
    Create the DATA array which contains the
    info for ray density
    '''
    global long_0

    mymap = Basemap(projection=projection, lon_0=long_0, lat_0=0)
    #npts=max(gr_x, gr_y)
    # grd[2]: longitude
    # grd[3]: latitude
    grd = mymap.makegrid(gr_x, gr_y, returnxy=True)

    lons, lats = mymap.gcpoints(lon1, lat1, lon2, lat2, npts)
    dist = locations2degrees(lat1, lon1, lat2, lon2)
    bap = int((dist - 97.0)*npts/dist)/2
    midlon = len(lons)/2
    midlat = len(lats)/2
    lons = lons[midlon-bap:midlon+1+bap]
    lats = lats[midlat-bap:midlat+1+bap]
    data = np.zeros([len(grd[2]), len(grd[3])])
    for i in range(len(lons)):
        xi, yi = point_finder(lons[i], lats[i], grd)
        # first one is latitude and second longitude
        try:
            #data[yi][xi] = dt/float(dist-97.0)
            data[yi][xi] += dt/len(lons)
        except Exception, e:
            print e

def main():
    unusable_name = 'unusable_%s.csv' % AREA
    unusable_data = csv.writer(open(unusable_name,'wb'))

    for i in range(2004,2013):
        m = Basemap(projection='merc',resolution=None,llcrnrlat=small_lat, \
        urcrnrlon=big_long,urcrnrlat=big_lat,llcrnrlon=small_long)
        fig = plt.figure(figsize=(11,11))
        ax = plt.subplot(111)
        m.drawmapboundary(fill_color='#1B1B1B')
    
        crime_dictionary = get_crimes(i,m,unusable_data)
        plot_crimes(crime_dictionary['nonviolent'],m)
        plot_crimes_violent(crime_dictionary['violent'],m) 
        plot_crimes_murder(crime_dictionary['murder'],m)
    
        plot_buildings(m)
        plot_streets(m)
        label_streets()

        output_file_name = '%s/%s_%s.png' % (AREA,AREA,str(i))
        ax.text(0.01,0.01,'https://github.com/oneschirm',color='white',\
        fontsize=12,transform=ax.transAxes)
        ax.text(0.9, 0.01, str(i),color='white', fontsize=20,transform=ax.transAxes)
        plt.savefig(output_file_name, dpi=150,bbox_inches='tight')

def drawMap(report):
	entries = get_site_status(report)
#	entries = {'unl.edu': [276, 0, 246, 0], 'desy.de': [107, 0, 0, 0], 'fnal.gov': [16, 0, 294, 0], 'N/A': [0, 0, 0, 0]}
	posList = get_positions(entries)

	#bounds = [(-60, -120), (60, 120)]
	#bounds = [(30, -10), (60, 40)]
	bounds = get_bounds(posList, margin = 10)

	matplotlib.pyplot.subplots_adjust(left=0, right=1, top=1, bottom=0, wspace=0, hspace=0)
	fig = matplotlib.pyplot.figure(figsize=(12, 6))
	ax = matplotlib.pyplot.subplot(111)
	m = Basemap(projection='cyl', lat_0=0, lon_0=0,
		llcrnrlon=bounds[0][0], llcrnrlat=bounds[0][1],
		urcrnrlon=bounds[1][0], urcrnrlat=bounds[1][1])

	map_positions(m, posList)
	#posList = remove_all_overlap(posList)
	#print posList

	m.bluemarble()
	for pos in posList:
		draw_pie(ax, pos['info'], (pos['x'], pos['y']), pos['size'])
		ax.text(pos['x']+5, pos['y']+5, pos['site'], color='white', fontsize=8)
	fig.savefig(os.path.expanduser('~/map.png'), dpi=300)

def run(FILE_NAME):
    
    with h5py.File(FILE_NAME, mode='r') as f:

        name = '/Grid/IRprecipitation'
        data = f[name][:]
        units = f[name].attrs['units']
        _FillValue = f[name].attrs['_FillValue']
        data[data == _FillValue] = np.nan
        data = np.ma.masked_where(np.isnan(data), data)

        
        # Get the geolocation data
        latitude = f['/Grid/lat'][:]
        longitude = f['/Grid/lon'][:]

        
    m = Basemap(projection='cyl', resolution='l',
                llcrnrlat=-90, urcrnrlat=90,
                llcrnrlon=-180, urcrnrlon=180)
    m.drawcoastlines(linewidth=0.5)
    m.drawparallels(np.arange(-90, 91, 45))
    m.drawmeridians(np.arange(-180, 180, 45), labels=[True,False,False,True])
    m.pcolormesh(longitude, latitude, data.T, latlon=True)
    cb = m.colorbar()    
    cb.set_label(units)

    basename = os.path.basename(FILE_NAME)
    plt.title('{0}\n{1}'.format(basename, name))
    fig = plt.gcf()
    # plt.show()
    pngfile = "{0}.py.png".format(basename)
    fig.savefig(pngfile)