Updated fillstates.py to add Alaska and Hawaii as map insets

Author: jsh9

examples/fillstates.py

@@ -2,18 +2,27 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from mpl_toolkits.basemap import Basemap as Basemap
-from matplotlib.colors import rgb2hex
+from matplotlib.colors import rgb2hex, Normalize
 from matplotlib.patches import Polygon
+from matplotlib.colorbar import ColorbarBase
 
 # Lambert Conformal map of lower 48 states.
-m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
+m = Basemap(llcrnrlon=-119,llcrnrlat=20,urcrnrlon=-64,urcrnrlat=49,
             projection='lcc',lat_1=33,lat_2=45,lon_0=-95)
-# draw state boundaries.
-# data from U.S Census Bureau
-# http://www.census.gov/geo/www/cob/st2000.html
-shp_info = m.readshapefile('st99_d00','states',drawbounds=True)
-# population density by state from
-# http://en.wikipedia.org/wiki/List_of_U.S._states_by_population_density
+
+# Mercator projection, for Alaska and Hawaii
+m_ = Basemap(llcrnrlon=-190,llcrnrlat=20,urcrnrlon=-143,urcrnrlat=46,
+            projection='merc',lat_ts=20)  # do not change these numbers
+
+#%% ---------   draw state boundaries  ----------------------------------------
+## data from U.S Census Bureau
+## http://www.census.gov/geo/www/cob/st2000.html
+shp_info = m.readshapefile('st99_d00','states',drawbounds=True,
+                           linewidth=0.45,color='gray')
+shp_info_ = m_.readshapefile('st99_d00','states',drawbounds=False)
+
+## population density by state from
+## http://en.wikipedia.org/wiki/List_of_U.S._states_by_population_density
 popdensity = {
 'New Jersey':  438.00,
 'Rhode Island':   387.35,
@@ -65,13 +74,13 @@
 'Montana':     2.39,
 'Wyoming':      1.96,
 'Alaska':     0.42}
-print(shp_info)
-# choose a color for each state based on population density.
+
+#%% -------- choose a color for each state based on population density. -------
 colors={}
 statenames=[]
-cmap = plt.cm.hot # use 'hot' colormap
+cmap = plt.cm.hot_r # use 'reversed hot' colormap
 vmin = 0; vmax = 450 # set range.
-print(m.states_info[0].keys())
+norm = Normalize(vmin=vmin, vmax=vmax)
 for shapedict in m.states_info:
     statename = shapedict['NAME']
     # skip DC and Puerto Rico.
@@ -80,18 +89,59 @@
         # calling colormap with value between 0 and 1 returns
         # rgba value.  Invert color range (hot colors are high
         # population), take sqrt root to spread out colors more.
-        colors[statename] = cmap(1.-np.sqrt((pop-vmin)/(vmax-vmin)))[:3]
+        colors[statename] = cmap(np.sqrt((pop-vmin)/(vmax-vmin)))[:3]
     statenames.append(statename)
-# cycle through state names, color each one.
+
+#%% ---------  cycle through state names, color each one.  --------------------
+fig = plt.gcf()  # get current figure instance
 ax = plt.gca() # get current axes instance
+
 for nshape,seg in enumerate(m.states):
     # skip DC and Puerto Rico.
-    if statenames[nshape] not in ['District of Columbia','Puerto Rico']:
-        color = rgb2hex(colors[statenames[nshape]]) 
+    if statenames[nshape] not in ['Puerto Rico', 'District of Columbia']:
+        color = rgb2hex(colors[statenames[nshape]])
         poly = Polygon(seg,facecolor=color,edgecolor=color)
         ax.add_patch(poly)
-# draw meridians and parallels.
-m.drawparallels(np.arange(25,65,20),labels=[1,0,0,0])
-m.drawmeridians(np.arange(-120,-40,20),labels=[0,0,0,1])
-plt.title('Filling State Polygons by Population Density')
+
+AREA_1 = 0.005  # exclude small Hawaiian islands that are smaller than AREA_1
+AREA_2 = AREA_1 * 30.0  # exclude Alaskan islands that are smaller than AREA_2
+AK_SCALE = 0.18  # scale down Alaska to show as a map inset
+
+for nshape, shapedict in enumerate(m_.states_info):  # plot Alaska and Hawaii as map insets
+    if shapedict['NAME'] in ['Alaska', 'Hawaii']:
+        seg = m_.states[int(shapedict['SHAPENUM'] - 1)]
+        if shapedict['NAME'] == 'Hawaii' and float(shapedict['AREA']) > AREA_1:
+            seg = list(map(lambda (x,y): (x-1900000, y+250000), seg))
+            color = rgb2hex(colors[statenames[nshape]])
+        elif shapedict['NAME'] == 'Alaska' and float(shapedict['AREA']) > AREA_2:
+            seg = list(map(lambda (x,y): (AK_SCALE*x-200000, AK_SCALE*y-650000), seg))
+            color = rgb2hex(colors[statenames[nshape]])
+        poly = Polygon(seg, facecolor=color, edgecolor='gray', linewidth=.45)
+        ax.add_patch(poly)
+
+ax.set_title('United states population density by state')
+
+#%% ---------  Plot bounding boxes for Alaska and Hawaii insets  --------------
+light_gray = [0.8]*3  # define light gray color RGB
+m_.plot(np.linspace(170,177),np.linspace(29,29),linewidth=1.,
+        color=light_gray,latlon=True)
+m_.plot(np.linspace(177,180),np.linspace(29,26),linewidth=1.,
+        color=light_gray,latlon=True)
+m_.plot(np.linspace(180,180),np.linspace(26,23),linewidth=1.,
+        color=light_gray,latlon=True)
+m_.plot(np.linspace(-180,-177),np.linspace(23,20),linewidth=1.,
+        color=light_gray,latlon=True)
+m_.plot(np.linspace(-180,-175),np.linspace(26,26),linewidth=1.,
+        color=light_gray,latlon=True)
+m_.plot(np.linspace(-175,-171),np.linspace(26,22),linewidth=1.,
+        color=light_gray,latlon=True)
+m_.plot(np.linspace(-171,-171),np.linspace(22,20),linewidth=1.,
+        color=light_gray,latlon=True)
+
+#%% ---------   Show color bar  ---------------------------------------
+ax_c = fig.add_axes([0.9, 0.1, 0.03, 0.8])
+cb = ColorbarBase(ax_c,cmap=cmap,norm=norm,orientation='vertical',
+                  label=r'[population per $\mathregular{km^2}$]')
+
 plt.show()
+