我们从Python开源项目中,提取了以下50个代码示例,用于说明如何使用mpl_toolkits.basemap.Basemap()。
def plot_grid2D(lons, lats, tec_grid2D, datetime, title_label = ''): LATS, LONS = np.meshgrid(lats, lons) m = Basemap(llcrnrlon=-180, llcrnrlat=-55, urcrnrlon=180, urcrnrlat=75, projection='merc', area_thresh=1000, resolution='i') m.drawstates() m.drawcountries() m.drawcoastlines() parallels = np.arange(-90,90,20) m.drawparallels(parallels,labels=[True,False,False,True]) meridians = np.arange(0,360,40) m.drawmeridians(meridians,labels=[True,False,False,True]) m.scatter(LONS, LATS, c=tec_grid2D, latlon = True, linewidths=0, s=5) m.colorbar() plt.title('%s\n%s' % (title_label, datetime.isoformat(' ')))
def prepare_basemap(min_lat, min_lon, max_lat, max_lon, delta_lat, delta_lon): """ :param min_lat: float :param min_lon: float :param max_lat: float :param max_lon: float :param delta_lat: float :param delta_lon: float :return: A Basemap instance """ m = Basemap(projection='cyl', llcrnrlat=min_lat, urcrnrlat=max_lat, urcrnrlon=max_lon, llcrnrlon=min_lon, resolution='l') m.drawcountries() m.drawcoastlines() m.drawparallels(np.arange(np.around(min_lat, 1), np.around(max_lat, 1), -round(delta_lat) / 5.0), labels=[1, 0, 0, 0]) m.drawmeridians(np.arange(np.around(min_lon, 1), np.around(max_lon, 1), -round(delta_lon) / 5.0), labels=[0, 0, 0, 1]) return m
def make_stations(network,latmin,latmax,lonmin,lonmax,dlat,dlon,plot=True,**kwargs): head = kwargs.get('head','RM') elev = kwargs.get('elev',0) dep = kwargs.get('dep',0) lats = np.arange(latmin,latmax+dlat,dlat) lons = np.arange(lonmin,lonmax+dlon,dlon) lats_pts,lons_pts = np.meshgrid(lats,lons) if plot: m = Basemap(projection='hammer',lon_0=0,resolution='l') m.drawcoastlines() m.drawmeridians(np.arange(0,351,10)) m.drawparallels(np.arange(-80,81,10)) lons_pts2,lats_pts2 = m(lons_pts,lats_pts) m.scatter(lons_pts2,lats_pts2) plt.show() f = open('STATIONS','w') print 'total number of stations : ',len(lats_pts) st = 1 for i in range(0,len(lats)): for j in range(0,len(lons)): f.write('{}{:04d} {} {:5.2f} {:5.2f} {:5.2f} {:5.2f}'.format(head,st,network,lats[i],lons[j],elev,dep)+'\n') st += 1
def plot_geo_config(stations,events): ''' plot geometrical configuration of tomography experiment ''' m = Basemap(projection='hammer',lon_0=0,resolution='l') m.drawcoastlines() m.drawmeridians(np.arange(0,361,30)) m.drawparallels(np.arange(-90,91,30)) station_lons = stations[0,:] station_lats = stations[1,:] x,y = m(station_lons,station_lats) m.scatter(x,y,s=50,marker='^',c='r') for event in events: lon = event[1] lat = event[2] x,y = m(lon,lat) m.scatter(x,y,s=100,marker='*',c='y') plt.savefig('geo_config.pdf',format='pdf')
def plot_pierce_points(self,depth=410.0,ax='None'): ''' Plots pierce points for a given depth. If an axes object is supplied as an argument it will use it for the plot. Otherwise, a new axes object is created. ''' if ax == 'None': m = Basemap(llcrnrlon=0.0,llcrnrlat=-35.0,urcrnrlon=120.0,urcrnrlat=35.0) else: m = ax m.drawmapboundary() m.drawcoastlines() found_points = False for i in range(0,len(self.pierce_dict)): if self.pierce_dict[i]['depth'] == depth: x,y = m(self.pierce_dict[i]['lon'],self.pierce_dict[i]['lat']) found_points = True if found_points == True: m.scatter(x,y,50,marker='+',color='r') else: print "no pierce points found for the given depth"
def plot_contiguous_US_tweets(lon, lat, file_path): ''' INPUT: List of longitudes (lon), list of latitudes (lat), file path to save the plot (file_path) OUTPUT: Plot of tweets in the contiguous US. ''' map = Basemap(projection='merc', resolution = 'h', area_thresh = 10000, llcrnrlon=-140.25, # lower left corner longitude of contiguous US llcrnrlat=5.0, # lower left corner latitude of contiguous US urcrnrlon=-56.25, # upper right corner longitude of contiguous US urcrnrlat=54.75) # upper right corner latitude of contiguous US x,y = map(lon, lat) map.plot(x, y, 'bo', markersize=2, alpha=.3) map.drawcoastlines() map.drawstates() map.drawcountries() map.fillcontinents(color = '#DAF7A6', lake_color='#a7cdf2') map.drawmapboundary(fill_color='#a7cdf2') plt.gcf().set_size_inches(15,15) plt.savefig(file_path, format='png', dpi=1000)
def main2(): from mpl_toolkits.basemap import Basemap import numpy as np import matplotlib.pyplot as plt # llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon # are the lat/lon values of the lower left and upper right corners # of the map. # resolution = 'i' means use intermediate resolution coastlines. # lon_0, lat_0 are the central longitude and latitude of the projection. m = Basemap(llcrnrlon=9.5, llcrnrlat=63.2, urcrnrlon=10.9, urcrnrlat=64., resolution='i', projection='tmerc', lon_0=10.7, lat_0=63.4) # can get the identical map this way (by specifying width and # height instead of lat/lon corners) # m = Basemap(width=894887,height=1116766,\ # resolution='i',projection='tmerc',lon_0=-4.36,lat_0=54.7) m.drawcoastlines() m.fillcontinents(color='coral', lake_color='aqua') # m.drawparallels(np.arange(-40, 61., 2.)) # m.drawmeridians(np.arange(-20., 21., 2.)) m.drawmapboundary(fill_color='aqua') plt.title("Transverse Mercator Projection") plt.show()
def _initialize_background_map(self): # calculate the correct aspect ratio by drawing a map of the substrate fig = plt.figure() ax = plt.subplot(111) nodes_to_map_coordinates = [MapCoordinate(lat=data['Latitude'], lon=data['Longitude']) for n, data in self._substrate.node.items()] lower_left, upper_right = ExtendedGraphVisualizer._get_corner_map_coordinates(nodes_to_map_coordinates) self._background_map = Basemap(resolution="c", projection='merc', anchor="W", ax=ax, llcrnrlat=lower_left.lat, urcrnrlat=upper_right.lat, llcrnrlon=lower_left.lon, urcrnrlon=upper_right.lon, fix_aspect=True) self._background_map.fillcontinents(ax=ax) (xmin, xmax), (ymin, ymax) = ax.get_xlim(), ax.get_ylim() dx, dy = float(abs(xmax - xmin)), float(abs(ymax - ymin)) self._map_width = (dx / dy) * MAP_HEIGHT plt.close(fig)
def plot_map(): lon=grd.lon; lat=grd.lat lat_ts=((lat[0,0]+lat[-1,0])*0.5) cbox=[lon[0, 0],lat[0, 0],lon[-1, -1],lat[-1, -1]] map = Basemap(projection='merc', resolution='l', lat_ts=lat_ts, llcrnrlon=cbox[0], llcrnrlat=cbox[1], urcrnrlon=cbox[2], urcrnrlat=cbox[3]) map.drawmapboundary(fill_color='white') map.drawcoastlines() map.fillcontinents(color='#cc9966', lake_color='#99ffff') lon_tcks = np.arange(-100.,30.,20.) lat_tcks = np.arange(0.,70.,10.) map.drawmeridians(lon_tcks, labels=[True, False, False, True]) map.drawparallels(lat_tcks, labels=[False, True, True, False]) # get projected coordinates x,y=map(lon,lat) return map, x, y
def basemap_raster_mercator(lon, lat, grid, cmin, cmax, cmap_name): # longitude/latitude extent lons = (np.amin(lon), np.amax(lon)) lats = (np.amin(lat), np.amax(lat)) # 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,figsize=(12,8),dpi=72) plt.axis('off') cmap = mpl.cm.get_cmap(cmap_name) m.pcolormesh(lon,lat,masked_grid,latlon=True,cmap=cmap,vmin=cmin,vmax=cmax) str_io = StringIO.StringIO() plt.savefig(str_io,bbox_inches='tight',format='png',pad_inches=0,transparent=True) plt.close() numpy_bounds = [ (lons[0],lats[0]),(lons[1],lats[0]),(lons[1],lats[1]),(lons[0],lats[1]) ] float_bounds = [ (float(x), float(y)) for x,y in numpy_bounds ] return str_io.getvalue(), float_bounds
def basemap_barbs_mercator(u,v,lat,lon): # lon/lat extents lons = (np.amin(lon), np.amax(lon)) lats = (np.amin(lat), np.amax(lat)) # 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) fig = plt.figure(frameon=False,figsize=(12,8),dpi=72*4) plt.axis('off') m.quiver(lon,lat,u,v,latlon=True) str_io = StringIO.StringIO() plt.savefig(str_io,bbox_inches='tight',format='png',pad_inches=0,transparent=True) plt.close() numpy_bounds = [ (lons[0],lats[0]),(lons[1],lats[0]),(lons[1],lats[1]),(lons[0],lats[1]) ] float_bounds = [ (float(x), float(y)) for x,y in numpy_bounds ] return str_io.getvalue(), float_bounds
def __init__(self,ASN): self.filename = str(ASN) + 'network graph' self.alias_filename = str(ASN) + 'alias_candidates.txt' self.ISP_Network = networkx.Graph() self.files_read = 0 self.border_points = set() if os.path.isfile(self.filename): load_file = shelve.open(self.filename, 'r') self.files_read = load_file['files_read'] self.ISP_Network = load_file['ISP_Network'] self.border_points = load_file['border_points'] self.ASN = ASN plt.ion() country = raw_input('enter country to focus on map [Israel/Usa/Australia/other]: ') if country == 'Israel' or country == 'israel' or country == 'ISRAEL': self.wmap = Basemap(projection='aeqd', lat_0 = 31.4, lon_0 = 35, width = 200000, height = 450000, resolution = 'i') elif country == 'USA' or country == 'usa': self.wmap = Basemap(projection='aeqd', lat_0 = 40, lon_0 = -98, width = 4500000, height = 2700000, resolution = 'i') elif country == 'Australia' or 'australia' or 'AUSTRALIA': self.wmap = Basemap(projection='aeqd', lat_0 = -23.07, lon_0 = 132.08, width = 4500000, height = 3500000, resolution = 'i') else: self.wmap = Basemap(projection='cyl', resolution = 'c') plt.hold(False)
def make_a_basemap(ax, ll_cnr, ru_cnr): bm = Basemap( projection=Choroplether.EPSG_4283_APPROXIMATION['projection'], ellps=Choroplether.EPSG_4283_APPROXIMATION['ellps'], lon_0=(ru_cnr[0] - ll_cnr[0]) / 2, lat_0=(ru_cnr[1] - ll_cnr[1]) / 2, llcrnrlon=ll_cnr[0], llcrnrlat=ll_cnr[1], urcrnrlon=ru_cnr[0], urcrnrlat=ru_cnr[1], lat_ts=0, resolution='i', suppress_ticks=True, ax=ax ) return bm
def get_basemap(shapefile): shp = fiona.open(shapefile+'.shp') bds = shp.bounds shp.close() extra = 0.01 ll = (bds[0], bds[1]) ur = (bds[2], bds[3]) coords = list(chain(ll, ur)) w, h = coords[2] - coords[0], coords[3] - coords[1] m = Basemap(projection='tmerc', lon_0= -118.2437, lat_0= 34.0522, ellps = 'WGS84', llcrnrlon=coords[0], # llcrnrlat=coords[1] - extra + 0.01 * h, llcrnrlat=33.6, urcrnrlon=coords[2], urcrnrlat=coords[3], lat_ts=0, resolution='i', suppress_ticks=True) m.readshapefile(shapefile, 'LA', color='black', zorder=2) return m, coords # Convenience functions for working with colour ramps and bars
def main(): bm = Basemap(projection = "rotpole", o_lat_p = 36.0, o_lon_p = 180.0, llcrnrlat = -10.590603, urcrnrlat = 46.591976, llcrnrlon = -139.08585, urcrnrlon = 22.661009, lon_0 = -106.0, rsphere = 6370000, resolution = 'l') fig = plt.figure(figsize=(8,8)) ax = fig.add_axes([0.1,0.1,0.8,0.8]) bm.drawcoastlines(linewidth=.5) print bm.proj4string plt.savefig("basemap_map.png") plt.close(fig)
def _basemap(self, geobounds, **kwargs): if not basemap_enabled(): return None local_kwargs = dict(projection = "lcc", lon_0 = self.stand_lon, lat_0 = self.moad_cen_lat, lat_1 = self.truelat1, lat_2 = self.truelat2, llcrnrlat = geobounds.bottom_left.lat, urcrnrlat = geobounds.top_right.lat, llcrnrlon = geobounds.bottom_left.lon, urcrnrlon = geobounds.top_right.lon, rsphere = Constants.WRF_EARTH_RADIUS, resolution = 'l') local_kwargs.update(kwargs) _basemap = Basemap(**local_kwargs) return _basemap
def _basemap(self, geobounds, **kwargs): if not basemap_enabled(): return None local_kwargs = dict(projection = "merc", lon_0 = self.stand_lon, lat_0 = self.moad_cen_lat, lat_ts = self._lat_ts, llcrnrlat = geobounds.bottom_left.lat, urcrnrlat = geobounds.top_right.lat, llcrnrlon = geobounds.bottom_left.lon, urcrnrlon = geobounds.top_right.lon, rsphere = Constants.WRF_EARTH_RADIUS, resolution = "l") local_kwargs.update(kwargs) _basemap = Basemap(**local_kwargs) return _basemap
def _basemap(self, **kwargs): if not basemap_enabled(): return None local_kwargs = dict(projection = "stere", lon_0 = self.stand_lon, lat_0 = self._hemi, lat_ts = self._lat_ts, llcrnrlat = geobounds.bottom_left.lat, urcrnrlat = geobounds.top_right.lat, llcrnrlon = geobounds.bottom_left.lon, urcrnrlon = geobounds.top_right.lon, rsphere = Constants.WRF_EARTH_RADIUS, resolution = "l") local_kwargs.update(kwargs) _basemap = Basemap(**local_kwargs) return _basemap
def _basemap(self, geobounds, **kwargs): if not basemap_enabled(): return None local_kwargs = dict(projection = "rotpole", o_lat_p = self._bm_cart_pole_lat, o_lon_p = self.pole_lon, llcrnrlat = geobounds.bottom_left.lat, urcrnrlat = geobounds.top_right.lat, llcrnrlon = geobounds.bottom_left.lon, urcrnrlon = geobounds.top_right.lon, lon_0 = self._bm_lon_0, rsphere = Constants.WRF_EARTH_RADIUS, resolution = "l") local_kwargs.update(kwargs) _basemap = Basemap(**local_kwargs) return _basemap
def create_basemap(self, projection='merc'): """Creates a very basic basemap: 'bbox': (x1, x2, y1, y2) If you want to use the power of basemap you can write your own basemap and pass it to you object """ bm = Basemap( llcrnrlat=self.bbox[2], urcrnrlat=self.bbox[3], llcrnrlon=self.bbox[0], urcrnrlon=self.bbox[1], projection=projection) bm.drawcoastlines(linewidth=0) return bm
def create_map(): """ create the base map for the choropleth. """ logging.debug('create_map() -- Please wait while I create the world.') degrees_width = 118.0 degrees_height = 55.0 center_lat = 44.5 center_lon = -110.0 my_map = Basemap( # ax=ax, projection='merc', # default is cyl ellps='WGS84', lat_0=center_lat, lon_0=center_lon, llcrnrlat=center_lat - degrees_height / 2.0, llcrnrlon=center_lon - degrees_width / 2.0, urcrnrlat=center_lat + degrees_height / 2.0, urcrnrlon=center_lon + degrees_width / 2.0, resolution='i', # 'c', 'l', 'i', 'h', 'f' ) logging.debug('created map') logging.debug('loading shapes...') for section_name in CONTEST_SECTIONS.keys(): # logging.debug('trying to load shape for %s', section_name) try: my_map.readshapefile('shapes/%s' % section_name, section_name, drawbounds=False) except IOError, err: logging.error('Could not load shape for %s' % section_name) logging.debug('loaded section shapes') return my_map
def Heatmap(self, grib_object, file_name): """Generate Heatmap with Data from grib object grib_object: an object containing raw data to be visualized file_name: a string representing the name of generated picture """ data = grib_object.values lat,lon = grib_object.latlons() unit = grib_object['units'] data_type = grib_object['name'] date = self.GetTime(grib_object) fig = plt.figure(figsize=(8,8)) ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) m = Basemap( resolution='c', # c, l, i, h, f or None projection='cyl', lat_0=39.72, lon_0=-86.29, llcrnrlon=-87.79, llcrnrlat= 38.22, urcrnrlon=-84.79, urcrnrlat=41.22) parallels = np.arange(38.22, 41.22, 0.5) m.drawparallels(parallels,labels=[False,True,True,False]) meridians = np.arange(-87.79, -84.79, 0.5) m.drawmeridians(meridians,labels=[True,False,False,True]) x,y = m(lon, lat) cs = m.pcolormesh(x,y,data, shading='flat', cmap=plt.cm.jet) cbar = plt.colorbar(cs,location='bottom', fraction=0.046, pad=0.06) # Adjust the position of Unit cbar_ax = cbar.ax cbar_ax.text(0.0, -1.3, unit, horizontalalignment='left') m.readshapefile(self.shapeFile,'areas') plt.title(data_type + ' ' + date, fontsize = 'x-large') plt.savefig(file_name) plt.close(fig)
def map_trace(tr,ax='None',showpath=True,showplot=True,Lat_0=0.0,Lon_0=60.0): from mpl_toolkits.basemap import Basemap if ax == 'None': m = Basemap(projection='ortho',lat_0=Lat_0,lon_0=Lon_0,resolution='l') m.drawmapboundary() m.drawcoastlines() m.fillcontinents(color='gray',lake_color='white') else: m = ax x1,y1 = m(tr.stats.sac['evlo'],tr.stats.sac['evla']) x2,y2 = m(tr.stats.sac['stlo'],tr.stats.sac['stla']) m.scatter(x1,y1,s=200.0,marker='*',facecolors='y',edgecolors='k',zorder=99) m.scatter(x2,y2,s=20.0,marker='^',color='b',zorder=99) if showpath == True: m.drawgreatcircle(tr.stats.sac['evlo'],tr.stats.sac['evla'], tr.stats.sac['stlo'],tr.stats.sac['stla'], linewidth=1,color='k',alpha=0.5) if showplot == True: plt.show() else: return m ######################################################################### # Plot map of earthquake and station for all traces in stream #########################################################################
def plot_eq(self,ax='None',showpath=True,showplot=True,lon_0=0.0,lat_0=0.0): from mpl_toolkits.basemap import Basemap if ax == 'None': #m = Basemap(projection='hammer',lon_0=self.ry) m = Basemap(projection='ortho',lat_0=lat_0,lon_0=lon_0,resolution='l') m.drawmapboundary() m.drawcoastlines() m.fillcontinents(color='gray',lake_color='white') else: m = ax x1,y1 = m(self.sy,90.0-self.sx) x2,y2 = m(self.ry,90.0-self.rx) m.scatter(x1,y1,s=200.0,marker='*',facecolors='y',edgecolors='k',zorder=99) m.scatter(x2,y2,s=20.0,marker='^',color='b',zorder=99) if showpath == True: lon_s = self.sy lat_s = 90.0-self.sx lon_r = self.ry lat_r = 90.0-self.rx print "lon_s,lat_s,lon_r,lat_r", lon_s, lat_s, lon_r, lat_r m.drawgreatcircle(lon_s,lat_s,lon_r,lat_r,linewidth=1,color='k',alpha=0.5) if showplot == True: plt.show() else: return m
def draw_basemap(resolution='l', ax=None,country_linewidth=0.5, coast_linewidth= 1.0, zorder=None, **kwds): if ax is None: ax = plt.gca() m = Basemap(*(ax.viewLim.min + ax.viewLim.max), resolution=resolution, ax=ax, **kwds) m.drawcoastlines(linewidth=coast_linewidth, zorder=zorder) m.drawcountries(linewidth=country_linewidth, zorder=zorder) return m
def _create_map(self, coordinates, labels, n_clusters): """ Graphs GPS coordinates on map (used ti verify clustering algo) """ map = Basemap(llcrnrlat=22, llcrnrlon=-119, urcrnrlat=49, urcrnrlon=-64, projection='lcc',lat_1=33,lat_2=45,lon_0=-95) base_color = 'white' border_color = 'lightgray' boundary_color = 'gray' map.fillcontinents(color=base_color, lake_color=border_color) map.drawstates(color=border_color) map.drawcoastlines(color=border_color) map.drawcountries(color=border_color) map.drawmapboundary(color=boundary_color) lat, lon = zip(*coordinates) map.scatter(lon, lat, latlon=True, s=1, alpha=0.5, zorder=2) plt.title('Estimated number of clusters: %d' % n_clusters) #Upload an image to Cloud Storage plt.savefig(constants.C_MAP_FPATH) if os.path.exists(constants.C_MAP_FPATH): print "HERERE" else: print "NOT HEEEEEEEEEEEEEEEEEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRREEEEEEEEEEEE"
def setProjection(self,gid,axs=None,west=None,north=None,east=None,south=None): i = gid-1 #If there is input on grid extent change the corners if (west != None and north != None and east != None and south != None): llcrnrlon = west llcrnrlat = south urcrnrlon = east urcrnrlat = north else: llcrnrlon = self.ll_lon[i] llcrnrlat = self.ll_lat[i] urcrnrlon = self.ur_lon[i] urcrnrlat = self.ur_lat[i] #GOES is not in a specific projection so using cyl # self.map[i] = Basemap(ax=axs,projection=self.projectionType, # llcrnrlon=self.ll_lon[i],llcrnrlat=self.ll_lat[i], # urcrnrlat = self.ur_lat[i],urcrnrlon = self.ur_lon[i], # resolution=self.resolution) self.map[i] = Basemap(ax=axs, projection=self.projectionType, lat_ts = self.lat0[i], llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat, urcrnrlon = urcrnrlon, urcrnrlat = llcrnrlon, resolution=self.resolution) self.maplon[i],self.maplat[i] = self.map[i].makegrid(self.nx[i],self.ny[i])
def setProjection(self,gid=None,axs=None,west=None,north=None,east=None,south=None): self.map = [None]*1 lon = self.variables['lon'] lat = self.variables['lat'] self.glons[0], self.glats[0] = np.meshgrid(lon,lat) self.nx = [self.glons[0].shape[1]] self.ny = [self.glons[0].shape[0]] self.dx = [50000.] self.dy = [62500.] if self.projectionType == None or self.projectionType == 'robin': self.projectionType = 'robin' self.lon0 = [0.0] self.lat0 = [0.0] self.ll_lon = [-179.375] self.ll_lat = [-90] self.ur_lon = [180] self.ur_lat = [90] self.map[0] = Basemap(ax=axs, projection=self.projectionType, lat_0 = self.lat0[0], lon_0 = self.lon0[0], llcrnrlon = self.ll_lon[0], llcrnrlat = self.ll_lat[0], urcrnrlon = self.ur_lon[0], urcrnrlat = self.ur_lat[0], resolution=self.resolution) else: self.lon0 = [270.0] self.lat0 = [0.0] self.map[0] = Basemap(ax=axs, projection=self.projectionType, lon_0 = self.lon0[0], boundinglat = self.lat0[0], resolution=self.resolution)
def setProjection(self,gid,axs=None,west=None,north=None,east=None,south=None): i = gid-1 #If there is input on grid extent change the corners if (west != None and north != None and east != None and south != None): llcrnrlon = west llcrnrlat = south urcrnrlon = east urcrnrlat = north else: llcrnrlon = self.ll_lon[i] llcrnrlat = self.ll_lat[i] urcrnrlon = self.ur_lon[i] urcrnrlat = self.ur_lat[i] #GOES is not in a specific projection so using cyl # self.map[i] = Basemap(ax=axs,projection=self.projectionType, # llcrnrlon=self.ll_lon[i],llcrnrlat=self.ll_lat[i], # urcrnrlat = self.ur_lat[i],urcrnrlon = self.ur_lon[i], # resolution=self.resolution) self.map[i] = Basemap(ax=axs, projection=self.projectionType, lat_ts = self.lat0[i], llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat, urcrnrlon = urcrnrlon, urcrnrlat = urcrnrlat, resolution=self.resolution) self.maplon[i],self.maplat[i] = self.map[i].makegrid(self.nx[i],self.ny[i]) #This function navigate the input GOES file # i.e. determines the lat/long for each pixel # Based on the GOES-R Product Definition and User's Guide (PUG)
def setProjection(self,gid=None,axs=None,west=None,north=None,east=None,south=None): #If there is input on grid extent change the corners if (west != None and north != None and east != None and south != None): llcrnrlon = west llcrnrlat = south urcrnrlon = east urcrnrlat = north else: #Hardcode for U.S. for now self.ll_lon = -122.0 self.ll_lat = 20.0 self.ur_lon = -63.0 self.ur_lat = 50.0 llcrnrlon = self.ll_lon llcrnrlat = self.ll_lat urcrnrlon = self.ur_lon urcrnrlat = self.ur_lat self.map = [None]*1 self.nx = [None] self.ny = [None] #Hardcode for U.S. self.lat0 = [38.0] self.lon0 = [-95.0] self.map[0] = Basemap(projection=self.projectionType, lat_0 = self.lat0[0], lon_0 = self.lon0[0], llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat, urcrnrlon = urcrnrlon, urcrnrlat = urcrnrlat, resolution=self.resolution, ax = axs) ############# GUI Related Functions ################## #Change the year #Sets the year
def setProjection(self,gid=None,axs=None,west=None,north=None,east=None,south=None): self.map = [None]*1 self.nx = [None] self.ny = [None] self.glons = [None]*1 self.glats = [None]*1 self.ll_lon = self.lon0[0] - 2.25 self.ll_lat = self.lat0[0] - 1.75 self.ur_lon = self.lon0[0] + 2.25 self.ur_lat = self.lat0[0] + 1.75 self.map[0] = Basemap(projection=self.projectionType, lat_0 = self.lat0[0], lon_0 = self.lon0[0], llcrnrlon = self.ll_lon, llcrnrlat = self.ll_lat, urcrnrlon = self.ur_lon, urcrnrlat = self.ur_lat, resolution=self.resolution, ax = axs) display = pyart.graph.RadarMapDisplay(self.radar) if ( self.var == 'velocity' or self.var == 'spectrum_width' ) and self.currentSweep == 0: self.currentSweep += 1 if (self.var == 'differential_reflectivity' or self.var == 'differential_phase' or self.var == 'cross_correlation_ratio') and self.currentSweep == 1: self.currentSweep -= 1 x,y = display._get_x_y(self.currentSweep,True,None) x0,y0 = self.map[0](self.lon0[0],self.lat0[0]) self.glons[0],self.glats[0] = self.map[0]((x0+x*1000.),(y0+y*1000.),inverse=True)
def Mapa(self): self.mapa = Basemap(projection='mill',lat_ts=10,llcrnrlon=self.lon0, urcrnrlon=self.lon1, llcrnrlat=self.lat0, urcrnrlat=self.lat1, resolution=self.res)
def build_map_obj(shapefile): shp = fiona.open(shapefile + '.shp') bds = shp.bounds shp.close() extra = 0.01 ll = (bds[0], bds[1]) ur = (bds[2], bds[3]) coords = list(chain(ll, ur)) w, h = coords[2] - coords[0], coords[3] - coords[1] m = Basemap( projection='tmerc', lon_0=(coords[0] + coords[2])/2, lat_0=(coords[1] + coords[3])/2, ellps='helmert', llcrnrlon=coords[0] - w * 0.01, llcrnrlat=coords[1] - extra + 0.01 * h, urcrnrlon=coords[2] + w * 0.01, urcrnrlat=coords[3] + extra + 0.01 * h, resolution='i', suppress_ticks=True ) m.readshapefile(shapefile, 'SF', color='black', zorder=2 ) return m, coords
def basemap(self, geobounds, **kwargs): """Return a :class:`matplotlib.mpl_toolkits.basemap.Basemap` object for the map projection. Args: geobounds (:class:`wrf.GeoBounds`, optional): The geobounds to get the extents. If set to None and using the *var* parameter, the geobounds will be taken from the variable. If using a file, then the geobounds will be taken from the native grid. **kwargs: Keyword arguments for creating a :class:`matplotlib.mpl_toolkits.basemap.Basemap`. By default, the domain bounds will be set to the native projection, the resolution will be set to 'l', and the other projection parameters will be set by the information in the file. Returns: :class:`matplotlib.mpl_toolkits.basemap.Basemap`: A Basemap object for the projection. See Also: :class:`matplotlib.mpl_toolkits.basemap.Basemap` """ if not basemap_enabled(): raise RuntimeError("'mpl_toolkits.basemap' is not " "installed or is disabled") return self._basemap(geobounds, **kwargs)
def _proj4(self): _proj4 = ("+proj=eqc +units=meters +a={} +b={} " "+lon_0={}".format(Constants.WRF_EARTH_RADIUS, Constants.WRF_EARTH_RADIUS, self.stand_lon)) return _proj4 # Notes (may not be correct since this projection confuses me): # Each projection system handles this differently. # 1) In WRF, if following the WPS instructions, POLE_LON is mainly used to # determine north or south hemisphere. In other words, it determines if # the globe is tipped toward or away from you. # 2) In WRF, POLE_LAT is always positive, but should be negative in the # proj4 based systems when using the southern hemisphere projections. # 3) In cartopy, pole_longitude is used to describe the dateline, which # is 180 degrees away from the normal central (standard) longitude # (e.g. center of the projection), according to the cartopy developer. # 4) In basemap, lon_0 should be set to the central (standard) longitude. # 5) In either cartopy, basemap or pyngl, I'm not sure that projections with # a pole_lon not equal to 0 or 180 can be plotted. Hopefully people # follow the WPS instructions, otherwise I need to see a sample file. # 6) For items in 3 - 4, the "longitude" (lon_0 or pole_longitude) is # determined by WRF's # STAND_LON values, with the following calculations based on hemisphere: # BASEMAP: NH: -STAND_LON; SH: 180.0 - STAND_LON # CARTOPY: NH: -STAND_LON - 180.; SH: -STAND_LON # 9) For PYNGL/NCL, you only need to set the center lat and center lon, # Center lat is the offset of the pole from +/- 90 degrees. Center # lon is -STAND_LON in NH and 180.0 - STAND_LON in SH. # 10) It also appears that NetCDF CF has no clear documentation on what # each parameter means. Going to assume it is the same as basemap, since # basemap appears to mirror the WMO way of doing things (tilt earth, then # spin globe). # 11) Basemap and cartopy produce projections that differ in their extent # calculations by either using negative values or 0-360 (basemap). For # this reason, the proj4 string for this class will use cartopy's values # to keep things in the -180 to 180, -90 to 90 range. # 12) This projection makes me sad.
def __init__(self, *args, **kwargs): """Initialisation of the map object :param *args: additional non-keyword parameters (passed to `basemap <http://matplotlib.org/basemap/api/basemap_api.html#mpl _toolkits.basemap.Basemap>`_) :param **kwargs: additional keyword parameters (passed to `basemap <http://matplotlib.org/basemap/api/basemap_api.html#mpl _toolkits.basemap.Basemap>`_) .. note:: Additional possible input in **kwargs wit to ``Basemap`` objects is "topo_data" which will be set in case it is a valid input (i.e. :class:`TopoData`) and will be used for plotting topography """ super(Map, self).__init__(*args, **kwargs) self.topo_data = None #IMPORTANT HANDLES self.contour_lines = None self.contour_filled = None self.colorbars = {} self.points = {} self.lines = {} self.polygons = {} self.texts = {} self.plotted_data = {} self.map_scale = None self.meridians = None self.parallels = None self.default_colors = {"contour_lines" : "#708090", "land" : "#FFE0B2", "water" : "#e6e6fa"}
def __init__(self, ): self.m = Basemap()
def map_world(self, res='l'): """ It creates a map of the world. :return: """ self.m = Basemap(projection='mill', resolution=res) # Draw the costal lines self.m.drawcoastlines() # Fill the continents with black self.m.fillcontinents(color='K')
def map_mediterranean(self, res='l'): """ It creates a map of the Mediterranean. :param res: Resolution of the map :return: """ self.m = Basemap(projection='mill', resolution=res, llcrnrlat=30, llcrnrlon=-13, urcrnrlat=47, urcrnrlon=38) # Draw the costal lines self.m.drawcoastlines() # Fill the continents with black self.m.fillcontinents(color='K')
def map(dataframe, title = "Map", colorbarName = None): fig, ax = plt.subplots(figsize=(80,40)) m = Basemap(resolution='l', # c, l, i, h, f or None projection='robin', lon_0=0) m.drawmapboundary(fill_color='#46bcec') m.fillcontinents(color='#f2f2f2',lake_color='#46bcec') # m.drawcoastlines() plt.title(title, fontsize=50, y=1.08) m.readshapefile('visualization/World/World', 'world',drawbounds=False) df_plot = pd.DataFrame({ 'shapes': [Polygon(np.array(shape), True) for shape in m.world], 'country': [country['ISO3'] for country in m.world_info] }) df_plot = df_plot.merge(dataframe, on='country', how='left') df_plot = df_plot.dropna() cmap = plt.get_cmap('RdYlGn') pc = PatchCollection(df_plot.shapes, zorder=2) norm = Normalize() pc.set_facecolor(cmap(norm(df_plot['value'].values))) ax.add_collection(pc) mapper = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap) mapper.set_array(df_plot['value']) cbar = plt.colorbar(mapper, shrink=0.7, label = colorbarName) fig = plt.gcf() fig.savefig("Map.jpg") plt.show()
def draw_training_points(self, filename='points.pdf', world=False, figsize=(4,3)): ''' draws training points on map ''' import matplotlib as mpl mpl.use('Agg') import matplotlib.patches as mpatches import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap, cm, maskoceans fig = plt.figure(figsize=figsize) lllat = 24.396308 lllon = -124.848974 urlat = 49.384358 urlon = -66.885444 if world: lllat = -90 lllon = -180 urlat = 90 urlon = 180 m = Basemap(llcrnrlat=lllat, urcrnrlat=urlat, llcrnrlon=lllon, urcrnrlon=urlon, resolution='c', projection='cyl') m.drawmapboundary(fill_color = 'white') m.drawcoastlines(linewidth=0.2) m.drawcountries(linewidth=0.2) ax = plt.gca() ax.xaxis.set_visible(False) ax.yaxis.set_visible(False) for spine in ax.spines.itervalues(): spine.set_visible(False) #fig = plt.figure() # figsize=(4,4.2) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) ax.spines['bottom'].set_visible(False) ax.spines['left'].set_visible(False) train_locs = self.df_train[['lat', 'lon']].values mlon, mlat = m(*(train_locs[:,1], train_locs[:,0])) #m.scatter(mlon, mlat, color='red', s=0.6) m.plot(mlon, mlat, 'r.', markersize=1) m.drawlsmask(land_color='lightgray',ocean_color="#b0c4de", lakes=True) plt.tight_layout() plt.savefig(filename) plt.close() print("the plot saved in " + filename)
def plot_map(outdict, e, savename, fig1=None, m=None): """ This function will plot the output data in a scatter plot over a map of the satellite path. Args: outdict (dict[str, obj]): Output dictionary from analyzebeacons. e (dict[str, obj]): Output dictionary from ephem_doponly. savename(:obj:`str`): Name of the file the image will be saved to. fig1(:obj:'matplotlib figure'): Figure. m (:obj:'basemap obj'): Basemap object """ t = outdict['time'] slat = e['site_latitude'] slon = e['site_longitude'] if fig1 is None: fig1 = plt.figure() plt.figure(fig1.number) latlim = [math.floor(slat-15.), math.ceil(slat+15.)] lonlim = [math.floor(slon-15.), math.ceil(slon+15.)] if m is None: m = Basemap(lat_0=slat, lon_0=slon, llcrnrlon=lonlim[0], llcrnrlat=latlim[0], urcrnrlon=lonlim[1], urcrnrlat=latlim[1]) m.drawcoastlines(color="gray") m.plot(slon, slat, "rx") scat = m.scatter(e["sublon"](t), e["sublat"](t), c=outdict['rTEC'], cmap='viridis', vmin=0, vmax=math.ceil(sp.nanmax(outdict['rTEC']))) plt.title('Map of TEC Over Satellite Path') cb = plt.colorbar(scat,label='rTEC in TECu') fig1.savefig(savename, dpi=300) #plt.draw() scat = m.scatter(e["sublon"](t), e["sublat"](t), c=outdict['rTEC_amp'], cmap='viridis', vmin=0, vmax=math.ceil(sp.nanmax(outdict['rTEC_amp']))) plt.title('Map of TEC_Amp Over Satellite Path') cb.set_clim(vmin=0, vmax=math.ceil(sp.nanmax(outdict['rTEC_amp']))) cb.draw_all() #plt.tight_layout() figpath,name = os.path.split(savename) savename = os.path.join(figpath,'amp'+name) fig1.savefig(savename, dpi=300) plt.close(fig1) #%% I/O for measurements
def Frame(self, grib_object, file_name, vmin, vmax): """ Initialize the first image of GIF grib_object: an object containing raw data to be visualized vmin: min of all data vmax: max of all data return generated figure instance """ fig = plt.figure(figsize=(8,8)) ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) lat,lon = grib_object.latlons() unit = grib_object['units'] data_type = grib_object['name'] date = self.GetTime(grib_object) data = grib_object.values m = Basemap( resolution='c', # c, l, i, h, f or None projection='cyl', lat_0=39.72, lon_0=-86.29, llcrnrlon=-87.79, llcrnrlat= 38.22, urcrnrlon=-84.79, urcrnrlat=41.22) parallels = np.arange(38.22, 41.22, 0.5) m.drawparallels(parallels,labels=[False,True,True,False]) meridians = np.arange(-87.79, -84.79, 0.5) m.drawmeridians(meridians,labels=[True,False,False,True]) x,y = m(lon, lat) im = m.pcolormesh(x,y,data, shading='flat', vmin=vmin, vmax=vmax, cmap=plt.cm.jet) # self.im = m.imshow(grib_object['values'],vmin=vmin,vmax=vmax, cmap=plt.cm.jet) cbar = plt.colorbar(location='bottom', fraction=0.046, pad=0.06) # Adjust the position of Unit cbar_ax = cbar.ax cbar_ax.text(0.0, -1.3, unit, horizontalalignment='left') m.readshapefile(self.shapeFile,'areas') title = plt.title(data_type + ' ' + date, fontsize = 'x-large') plt.savefig(file_name) plt.close(fig)
def plot_pierce_points(self,depth=410.0,ax='None',**kwargs): #################################################################################### ''' Plots pierce points for a given depth. If an axes object is supplied as an argument it will use it for the plot. Otherwise, a new axes object is created. kwargs: depth: pierce point depth proj: map projection. if none given, a default basemap axis is made, defined by the coordinates of the corners. if 'ortho' is given, a globe centered on Lat_0 and Lon_0 is made. ll_corner_lat: latitude of lower left corner of map ll_corner_lon: longitude of lower left corner of map ur_corner_lat: latitude of upper right corner of map Lat_0: latitude center of ortho Lon_0: logitude center of ortho return_ax: whether or you return the basemap axis object. default False ''' proj=kwargs.get('proj','default') ll_lat=kwargs.get('ll_corner_lat',-35.0) ll_lon=kwargs.get('ll_corner_lon',0.0) ur_lat=kwargs.get('ur_corner_lat',35.0) ur_lon=kwargs.get('ur_corner_lon',120.0) Lat_0=kwargs.get('Lat_0',0.0) Lon_0=kwargs.get('Lon_0',0.0) return_ax =kwargs.get('return_ax',False) if ax == 'None': if proj == 'default': m = Basemap(llcrnrlon=ll_lon,llcrnrlat=ll_lat,urcrnrlon=ur_lon,urcrnrlat=ur_lat) elif proj == 'ortho': m = Basemap(projection='ortho',lat_0=Lat_0,lon_0=Lon_0) m.drawmapboundary() m.drawcoastlines() m.fillcontinents() else: m = ax found_points = False for ii in self.pierce: if ii['depth'] == depth: x,y = m(ii['lon'],ii['lat']) found_points = True if found_points == True: m.scatter(x,y,100,marker='+',color='r',zorder=99) else: print "no pierce points found for the given depth" ###############################################################################
def map_stream(st,showpath=True,showplot=True,Lat_0=0.0,Lon_0=60.0): from mpl_toolkits.basemap import Basemap fig = plt.figure() plt.ion() m = Basemap(projection='ortho',lat_0=Lat_0,lon_0=Lon_0,resolution='l') m.drawmapboundary() m.drawcoastlines() m.fillcontinents(color='gray',lake_color='white') gc_lines = [] for tr in st: x1,y1 = m(tr.stats.sac['evlo'],tr.stats.sac['evla']) x2,y2 = m(tr.stats.sac['stlo'],tr.stats.sac['stla']) m.scatter(x1,y1,s=200.0,marker='*',facecolors='y',edgecolors='k',zorder=99) station_pt = m.scatter(x2,y2,s=20.0,marker='^',color='b',zorder=99,picker=1) station_pt.set_label(tr.stats.station) if showpath == True: gc_line = m.drawgreatcircle(tr.stats.sac['evlo'],tr.stats.sac['evla'], tr.stats.sac['stlo'],tr.stats.sac['stla'], linewidth=1,color='k',alpha=0.5) gc_line[0].set_label(tr.stats.station) gc_lines.append(gc_line) def onpick(event): art = event.artist picked = art.get_label() print "removing station(s) ", picked st_to_remove = st.select(station=picked) for killed_tr in st_to_remove: st.remove(killed_tr) art.remove() for gc in gc_lines: gc_line_label = gc[0].get_label() if gc_line_label == picked: gc[0].remove() fig.canvas.draw() fig.canvas.mpl_connect('pick_event', onpick) if showplot == True: plt.show() else: return m ######################################################################### # Plot multiple traces ontop of eachother, aligned on a phase and windowed #########################################################################
def find_pierce_coor(self,plot='False'): import geopy from geopy.distance import VincentyDistance ''' given an instance of the receiver function class this function returns latitude and longitude of all receiver side pierce points of Pds in a given depth range (the default range is 50 - 800 km) NOTE: be careful that ses3d typically uses colatitude, while this function returns latitude ''' depth_range = np.arange(50,800,5) #set range of pierce points #geodetic info bearing = self.az lon_s = self.ses3d_seismogram.sy lat_s = 90.0-self.ses3d_seismogram.sx lon_r = self.ses3d_seismogram.ry lat_r = 90.0-self.ses3d_seismogram.rx origin = geopy.Point(lat_s, lon_s) #find how far away the pierce point is model = TauPyModel(model='pyrolite_5km') for i in range(0,len(depth_range)): phase = 'P'+str(depth_range[i])+'s' pierce = model.get_pierce_points(self.eq_depth,self.delta_deg,phase_list=[phase]) points = pierce[0].pierce for j in range(0,len(points)): if points[j]['depth'] == depth_range[i] and points[j]['dist']*(180.0/np.pi) > 25.0: prc_dist = points[j]['dist']*(180.0/np.pi) d_km = prc_dist * ((2*np.pi*6371.0/360.0)) destination = VincentyDistance(kilometers=d_km).destination(origin,bearing) lat = destination[0] lon = destination[1] value = 0 row = {'depth':depth_range[i],'dist':prc_dist,'lat':lat,'lon':lon,'value':value} self.pierce_dict.append(row) if plot=='True': m = Basemap(projection='hammer',lon_0=0) m.drawmapboundary() m.drawcoastlines() m.drawgreatcircle(lon_s,lat_s,lon_r,lat_r,linewidth=1,color='b',alpha=0.5) for i in range(len(self.pierce_dict)): x,y = m(self.pierce_dict[i]['lon'],self.pierce_dict[i]['lat']) m.scatter(x,y,5,marker='o',color='r') plt.show()
def basicsurfmap(file,varstr,outputdir): """ Creates plots for each timestep of a file for the specified variable and outputs them to a specified directory """ import matplotlib matplotlib.use('Agg') import numpy as np import sys from netCDF4 import Dataset import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap # read in file and vars f = Dataset(file,'r') wrfeta = f.variables['ZNU'][0][:] times = f.variables['Times'][:] wrflats = f.variables['XLAT'][0][:] wrflons = f.variables['XLONG'][0][:] var = f.variables[varstr][:] # four corners of domain print wrflats.shape print wrflons.shape print wrflons[0].shape wrflat_s = wrflats[0,len(wrflats)-1] wrflat_n = wrflats[len(wrflats)-1,len(wrflons[0])-1] wrflon_w = wrflons[0,0] wrflon_e = wrflons[len(wrflats)-1,len(wrflons[0])-1] z = 0 # assuming lowest level of model # set up map map = Basemap(projection='merc',llcrnrlon=wrflon_w,urcrnrlon=wrflon_e,llcrnrlat=wrflat_s,urcrnrlat=wrflat_n,resolution='i') map.drawstates() map.drawcounties() map.drawcoastlines() x,y = map(wrflons,wrflats) # loop through times for t in range(len(times)): timestr = ''.join(times[t,:]) map.drawstates() map.drawcounties() map.drawcoastlines() plt1 = map.pcolormesh(x,y,var[t,z,:,:],vmin=np.amin(var),vmax=np.amax(var)) colorbar = map.colorbar(plt1,"right", size="5%",pad="2%") colorbar.set_label(f.variables[varstr].description+' '+f.variables[varstr].units) plt.title('WRF output valid: '+timestr) plt.savefig(outputdir+'/%03d_' % (t) +timestr+'_'+varstr+'.png') plt.clf()
def CO2nWind(file,outputdir): """ Outputs plots of CO2 in the lowest level and 10m winds """ import matplotlib matplotlib.use('Agg') import numpy as np import sys from netCDF4 import Dataset import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap # read in file and vars f = Dataset(file,'r') wrfeta = f.variables['ZNU'][0][:] times = f.variables['Times'][:] wrflats = f.variables['XLAT'][0][:] wrflons = f.variables['XLONG'][0][:] var = f.variables['CO2_ANT'][:,0,:,:] u = f.variables['U'][:,0,:,:] v = f.variables['V'][:,0,:,:] # destagger u/v u = (u[:,:,:-1] + u[:,:,1:])/2. v = (v[:,:-1,:] + v[:,1:,:])/2. # four corners of domain wrflat_s = wrflats[0,len(wrflats)-1] wrflat_n = wrflats[len(wrflats)-1,len(wrflons[0])-1] wrflon_w = wrflons[0,0] wrflon_e = wrflons[len(wrflats)-1,len(wrflons[0])-1] z = 0 # assuming lowest level of model # set up map map = Basemap(projection='merc',llcrnrlon=wrflon_w,urcrnrlon=wrflon_e,llcrnrlat=wrflat_s,urcrnrlat=wrflat_n,resolution='i') map.drawstates() map.drawcounties() map.drawcoastlines() x,y = map(wrflons,wrflats) # loop through times for t in range(len(times)): timestr = ''.join(times[t,:]) map.drawstates(color='gray',linewidth=1) map.drawcounties(color='white') map.drawcoastlines(color='gray',linewidth=1) plt1 = map.pcolormesh(x,y,var[t,:,:],vmin=380,vmax=450) #plt1 = map.pcolormesh(x,y,var[t,:,:],vmin=np.amin(var),vmax=np.amax(var)) winds = map.barbs(x[::20,::20],y[::20,::20],u[t,::20,::20]*1.94,v[t,::20,::20]*1.94,length=6,color='white') # *1.94 to convert m/s to knots (barb convention) colorbar = map.colorbar(plt1,"right", size="5%",pad="2%") colorbar.set_label(f.variables['CO2_ANT'].description+' '+f.variables['CO2_ANT'].units) plt.title('WRF output valid: '+timestr) plt.savefig(outputdir+'/%03d_' % (t) +timestr+'_CO2_wind.png') plt.clf()
def load(self): """loads shapefile onto graphical representation of data using basemap and fiona""" shape = fiona.open("data/shapefiles/chicago.shp") bounds = shape.bounds extra = 0.01 lower_left = (bounds[0], bounds[1]) upper_right = (bounds[2], bounds[3]) coords = list(chain(lower_left, upper_right)) width, height = coords[2] - coords[0], coords[3] - coords[1] self.base_map = Basemap( projection="tmerc", lon_0=-87., lat_0=41., ellps="WGS84", llcrnrlon=coords[0] - extra * width, llcrnrlat=coords[1] - extra + 0.01 * height, urcrnrlon=coords[2] + extra * width, urcrnrlat=coords[3] + extra + 0.01 * height, lat_ts=0, resolution='i', suppress_ticks=True ) self.base_map.readshapefile( "data/shapefiles/chicago", 'chicago', color='none', zorder=2 ) self.data_map = pd.DataFrame({ 'poly': [Polygon(xy) for xy in self.base_map.chicago], 'community_name': [ward['community'] for ward in self.base_map.chicago_info]}) self.data_map['area_m'] = self.data_map['poly'].map(lambda x: x.area) self.data_map['area_km'] = self.data_map['area_m'] / 100000 self.data_map['patches'] = self.data_map['poly'].map(lambda x: PolygonPatch(x, fc='#555555', ec='#787878', lw=.25, alpha=.9, zorder=4)) plt.close() self.fig = plt.figure() self.ax = self.fig.add_subplot(111, axisbg='w', frame_on=False) self.ax.add_collection(PatchCollection(self.data_map['patches'].values, match_original=True)) self.base_map.drawmapscale( coords[0] + 0.08, coords[1] + 0.015, coords[0], coords[1], 10., barstyle='fancy', labelstyle='simple', fillcolor1='w', fillcolor2='#555555', fontcolor='#555555', zorder=5)
def main(): args = get_arguments() path = load_path(args.eclipse_path_data) # The funny constants define the lower left and upper right # corners (lat and lon) of the map boundary map = Basemap(llcrnrlat=22, llcrnrlon=-119, urcrnrlat=49, urcrnrlon=-64, projection='lcc',lat_1=33,lat_2=45,lon_0=-95, resolution='l') base_color = 'white' border_color = 'lightgray' boundary_color = 'gray' map.fillcontinents(color=base_color, lake_color=border_color) map.drawstates(color=border_color) map.drawcoastlines(color=border_color) map.drawcountries(color=border_color) map.drawmapboundary(color=boundary_color) points = pickle.load(open(args.inside_filename)) lats = [] lons = [] colors = [] for point in points: lat = point[0] lon = point[1] colors.append('r') lats.append(lat) lons.append(lon) points = pickle.load(open(args.outside_filename)) for point in points: lat = point[0] lon = point[1] colors.append('g') lats.append(lat) lons.append(lon) # Draw locations as colored points (red for "inside" and green for # "outside) map.scatter(lons, lats, marker='.', c=colors, edgecolors='none', s=3, latlon=True, zorder=2, cmap=cm.plasma) fig = matplotlib.pyplot.gcf() fig.set_size_inches(18.5, 10.5) xs = [] ys = [] # Draw the eclipse path boundary as path for point in path.eclipse_boundary.exterior.coords: ys.append(point[0]) xs.append(point[1]) map.plot(xs, ys, latlon=True, alpha=0.5, zorder=3) plt.savefig(args.output, dpi=100)