Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Jan 2, 2018
Merged

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

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
cf8a3cb
Updated fillstates.py to add Alaska and Hawaii as map insets
Aug 7, 2017
cfff23e
Some changed according to suggestions
Aug 9, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
86 changes: 65 additions & 21 deletions examples/fillstates.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,18 +2,29 @@
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

fig, ax = plt.subplots()

# 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,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can this be read in only once and the drawbounds be a part of drawing the polygons?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not really. Alaska and Hawaii are plotted using Mercator projection, while the other contiguous states are plotted using Lambert projection. The reason for choosing Mercator projection to plot Alaska and Hawaii is purely an aesthetic one, and I have noticed that this map-making choice (i.e., Alaska and Hawaii being plotted by Mercator projection) was adopted from as early as half a century ago.

I am not sure what you meant by "drawbounds be a part of drawing the polygons". Line 20 is in the original example script.

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,
Expand Down Expand Up @@ -65,13 +76,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.
Expand All @@ -80,18 +91,51 @@
# 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.
ax = plt.gca() # get current axes instance

#%% --------- cycle through state names, color each one. --------------------
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.19 # scale down Alaska to show as a map inset
HI_OFFSET_X = -1900000 # X coordinate offset amount to move Hawaii "beneath" Texas
HI_OFFSET_Y = 250000 # similar to above: Y offset for Hawaii
AK_OFFSET_X = -250000 # X offset for Alaska (These four values are obtained
AK_OFFSET_Y = -750000 # via manual trial and error, thus changing them is not recommended.)

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 = [(x + HI_OFFSET_X, y + HI_OFFSET_Y) for x, y in seg]
color = rgb2hex(colors[statenames[nshape]])
elif shapedict['NAME'] == 'Alaska' and float(shapedict['AREA']) > AREA_2:
seg = [(x*AK_SCALE + AK_OFFSET_X, y*AK_SCALE + AK_OFFSET_Y)\
for x, y in 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
x1,y1 = m_([-190,-183,-180,-180,-175,-171,-171],[29,29,26,26,26,22,20])
x2,y2 = m_([-180,-180,-177],[26,23,20]) # these numbers are fine-tuned manually
m_.plot(x1,y1,color=light_gray,linewidth=0.8) # do not change them drastically
m_.plot(x2,y2,color=light_gray,linewidth=0.8)

#%% --------- 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}$]')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not a regular colorbar creation approach instead of manually creating axes and such? An example like this is already getting pretty complicated, so we should aim to reduce complexity where-ever possible.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried the regular m.colorbar() way, but I couldn't figure out what I should use for mappable. Please provide advice if you can. Thanks.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @jsh9:

I have this file around where I was coloring countries based on population. I used it for teaching purposes, so there is an intentional bug somewhere, but you can check out how I did the colorbar there at the end:

https://github.com/guziy/CMOS-python-tutorial/blob/master/fiona_test.py

Cheers


plt.show()