How to use Basemap (Python) to plot US with 50 states?
For anyone interested, I was able to fix it myself. The (x,y) coordinates of each segment (for Alaska and Hawaii) should be translated. I also scale down Alaska to 35% before translating it.
The second for-loop should be modified as following:
for nshape,seg in enumerate(m.states):
# skip DC and Puerto Rico.
if statenames[nshape] not in ['Puerto Rico', 'District of Columbia']:
# Offset Alaska and Hawaii to the lower-left corner.
if statenames[nshape] == 'Alaska':
# Alaska is too big. Scale it down to 35% first, then transate it.
seg = list(map(lambda (x,y): (0.35*x + 1100000, 0.35*y-1300000), seg))
if statenames[nshape] == 'Hawaii':
seg = list(map(lambda (x,y): (x + 5100000, y-900000), seg))
color = rgb2hex(colors[statenames[nshape]])
poly = Polygon(seg,facecolor=color,edgecolor=color)
ax.add_patch(poly)
Here is the new US map (using the 'Greens' colormap).
The above answer is great and was very helpful for me.
I noticed that there are many tiny islands that extend for many miles beyond the 8 main islands of Hawaii. These create little dots in Arizona, California, and Oregon, (or Nevada and Idaho) depending on how you translated Hawaii. To remove these, you need a condition on the area of the polygon. It's helpful to do one loop through the states_info
object to do this:
# Hawaii has 8 main islands but several tiny atolls that extend for many miles.
# This is the area cutoff between the 8 main islands and the tiny atolls.
ATOLL_CUTOFF = 0.005
m = Basemap(llcrnrlon=-121,llcrnrlat=20,urcrnrlon=-62,urcrnrlat=51,
projection='lcc',lat_1=32,lat_2=45,lon_0=-95)
# load the shapefile, use the name 'states'
m.readshapefile('st99_d00', name='states', drawbounds=True)
ax = plt.gca()
for i, shapedict in enumerate(m.states_info):
# Translate the noncontiguous states:
if shapedict['NAME'] in ['Alaska', 'Hawaii']:
seg = m.states[int(shapedict['SHAPENUM'] - 1)]
# Only include the 8 main islands of Hawaii so that we don't put dots in the western states.
if shapedict['NAME'] == 'Hawaii' and float(shapedict['AREA']) > ATOLL_CUTOFF:
seg = list(map(lambda (x,y): (x + 5200000, y-1400000), seg))
# Alaska is large. Rescale it.
elif shapedict['NAME'] == 'Alaska':
seg = list(map(lambda (x,y): (0.35*x + 1100000, 0.35*y-1300000), seg))
poly = Polygon(seg, facecolor='white', edgecolor='black', linewidth=.5)
ax.add_patch(poly)