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)
python类Basemap()的实例源码
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)
