I'm using Basemap to display the East Coast of the US and the Atlantic Ocean. Shapefile is used to read a shapefile the contents of which is loaded into a PyPlot subplot hosted Line Collection which overlays the ocean with a series of closed polygons. PyPlot text is used to label each polygon with it's designator.
I have two questions both of which involve the polygons I'm loading from the shapefile.
1.) In order to simplify polygon\shoreline intersections the sides of the polygons were brought onto shore to close the polygon. How do I make these lines not show up on the land portion of the basemap? I tried zorder but it didn't seem to make a difference.
2.) Is there a way to color fill the polygons being created in the LineCollection?
Thank you for the assist. PNG available by PM if needed.
Jim
import shapefile
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib import cm
from mpl_toolkits.basemap import Basemap
fig = plt.figure(figsize=(11.7,8.3))
#Custom adjust of the subplots
plt.subplots_adjust(left=0.05,right=0.95,top=0.90,bottom=0.05,wspace=0.15,hspace=0.05)
ax = plt.subplot(111)
#Let's create a basemap of the Northeast Region
x1 = -78.
x2 = -64.
y1 = 34.
y2 = 46.
m = Basemap(resolution='i',projection='merc', llcrnrlat=y1,urcrnrlat=y2,llcrnrlon=x1,urcrnrlon=x2,lat_ts=(x1+x2)/2)
m.drawcountries(linewidth=0.5,zorder = 10)
m.drawcoastlines(linewidth=0.5,zorder = 11)
m.drawstates(linewidth=0.5,zorder = 12)
##m.drawrivers(linewidth=0.5)
m.fillcontinents(color='coral')
sf = shapefile.Reader("c:\\temp\\US_Stat_Areas_no_Coastline\\US_Stat_Areas_no_Coastline")
shapeRecs = sf.shapeRecords()
centroids = []
for area in shapeRecs:
shpsegs = []
shpinfo = []
points = []
stat_area = area.record[0]
if int(stat_area) > 700:
continue
else:
if int(stat_area) <= 463:
continue
else:
pass
print "Processing Stat Area: %s" % stat_area
vertices = area.shape.points
lons = []
lats = []
for each_vertice in vertices:
lons.append(each_vertice[0])
lats.append(each_vertice[1])
x, y = m(lons, lats)
shpsegs.append(zip(x,y))
shpseg = zip(x,y)
for point in shpseg:
points.append(Point(point[0],point[1]))
lines = LineCollection(shpsegs,antialiaseds=(1,),zorder = 5)
lines.set_edgecolors('k')
lines.set_linewidth(0.3)
ax.add_collection(lines)
plt.savefig('test_plot.png',dpi=300)