For the conterminous states, using geoDataFrame.plot() is done as usual. But Alaska and Hawaii must be plotted individually as inset maps. In the complete and runnable code below there are many inserted comments that explain the important steps to help the readers. About the thematic data, I prefer to use more realistic population data (rather than the makeup data in the question) to demonstrate the thematic mapping on the plot. Package mapclassify is used to classify data for proper thematic data processing.
import matplotlib.pyplot as plt
import matplotlib
import cartopy
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import shapely.geometry as sgeom #for box drawing
import geopandas as gpd
import numpy as np
import mapclassify as mc
import pandas as pd
import requests
import json
# population per sq-km
# include `Puerto Rico` but not used
popdensity = {
'Alaska': 0.8,
'District of Columbia': 4251,
'Hawaii': 86,
'Puerto Rico': 360,
'New Jersey': 438.00,
'Rhode Island': 387.35,
'Massachusetts': 312.68,
'Connecticut': 271.40,
'Maryland': 209.23,
'New York': 155.18,
'Delaware': 154.87,
'Florida': 114.43,
'Ohio': 107.05,
'Pennsylvania': 105.80,
'Illinois': 86.27,
'California': 83.85,
'Virginia': 69.03,
'Michigan': 67.55,
'Indiana': 65.46,
'North Carolina': 63.80,
'Georgia': 54.59,
'Tennessee': 53.29,
'New Hampshire': 53.20,
'South Carolina': 51.45,
'Louisiana': 39.61,
'Kentucky': 39.28,
'Wisconsin': 38.13,
'Washington': 34.20,
'Alabama': 33.84,
'Missouri': 31.36,
'Texas': 30.75,
'West Virginia': 29.00,
'Vermont': 25.41,
'Minnesota': 23.86,
'Mississippi': 23.42,
'Iowa': 20.22,
'Arkansas': 19.82,
'Oklahoma': 19.40,
'Arizona': 17.43,
'Colorado': 16.01,
'Maine': 15.95,
'Oregon': 13.76,
'Kansas': 12.69,
'Utah': 10.50,
'Nebraska': 8.60,
'Nevada': 7.03,
'Idaho': 6.04,
'New Mexico': 5.79,
'South Dakota': 3.84,
'North Dakota': 3.59,
'Montana': 2.39,
'Wyoming': 1.96}
# use this simple colormap
my_colormap = matplotlib.cm.Reds
# some settings
edgecolor = "gray"
# use this column for thematic mapping
theme_value = "pop_per_sqkm"
# A function that draws inset map
# ===============================
def add_insetmap(axes_extent, map_extent, state_name, facecolor, edgecolor, geometry):
# create new axes, set its projection
use_projection = ccrs.Mercator() # preserves shape
#use_projection = ccrs.PlateCarree() # large distortion in E-W, bad for for Alaska
geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
sub_ax = plt.axes(axes_extent, projection=use_projection) # normal units
sub_ax.set_extent(map_extent, geodetic) # map extents
# option to add basic land, coastlines of the map
# can comment out if you don't need them
sub_ax.add_feature(cartopy.feature.LAND)
sub_ax.coastlines()
sub_ax.set_title(state_name)
# add map `geometry`
sub_ax.add_geometries([geometry], ccrs.PlateCarree(), \
facecolor=facecolor, edgecolor=edgecolor, lw=0.3)
# +++ more features can be added here +++
# plot box around the map
extent_box = sgeom.box(map_extent[0], map_extent[2], map_extent[1], map_extent[3])
sub_ax.add_geometries([extent_box], ccrs.PlateCarree(), color='none')
# access USA shapefile
# use data from internet
URL = 'https://github.com/kjhealy/us-county/raw/master/data/geojson/gz_2010_us_040_00_500k.json'
# more hard-code for settings
state_name = "NAME" #defined column header
fname = "usa_json.json"
# Only on the first-run, set True in the next statement
if True:
us_json = requests.get(URL).json()
with open(fname, 'w') as file:
file.write(json.dumps(us_json))
# saved it to disk as fname
# then load from file
usa = gpd.read_file(fname)
usa.sort_values(state_name, inplace=True)
# make dataframe from population data `popdensity` dict object
usa_popden = pd.DataFrame.from_dict(popdensity, orient='index',
columns=['pop_per_sqkm'])
# merge `usa_popden` to main dataframe, `usa`
newusa = usa.merge(usa_popden, how='left', left_on='NAME', right_index=True)
# take only some columns in `newusa` for our operation
newusa = newusa[['NAME', 'pop_per_sqkm', 'geometry']]
# Data classification for thematic mapping
# choose number of classes of population density
# classes --> assigned colors in thematic mapping
num_classes = 7 #will get class-values: 0,1,2,3,4,5,6
sclass = mc.Quantiles(newusa["pop_per_sqkm"].values, k=num_classes)
#print(sclass)
# add new column, "sclass", for raw class values, sclass.yb
# its values will be used to assign color for polygon's facecolor
newusa["sclass"] = sclass.yb
# extract parts of the whole 'newusa' geodataframe for separate plotting/manipulation
# 'usa_main': excluding non-conterminous states
usa_main = newusa[~newusa[state_name].isin(["Alaska", "Hawaii", "Puerto Rico"])] # exclude these
# re-project usa_main to equal-area conic projection "EPSG:2163"
usa_main.crs = {'init': 'epsg:4326'}
usa_main = usa_main.to_crs(epsg=2163)
# 'usa_more': non-conterminous states, namely, Alaska and Hawaii
usa_more = newusa[newusa[state_name].isin(["Alaska", "Hawaii"])] # include these
# ------------ Plot --------------
# plot 1st part, using usa_main and grab its axis as 'ax2'
ax2 = usa_main.plot(column="sclass", legend=False,
cmap=matplotlib.cm.Reds, ec=edgecolor, lw=0.4)
# manipulate colorbar/legend
fig = ax2.get_figure()
cax = fig.add_axes([0.9, .25, 0.02, 0.5]) #[left,bottom,width,height]
sm = plt.cm.ScalarMappable(cmap=my_colormap,
norm=plt.Normalize(vmin=min(newusa["sclass"]),vmax=max(newusa["sclass"])))
# clear the array of the scalar mappable
sm._A = []
cb = fig.colorbar(sm, cax=cax)
cb.set_label("Pop-density class")
# manipulate the axis seetings
ax2.set_frame_on(False)
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_xticklabels([])
ax2.set_yticklabels([])
ax2.set_title("Population Density Plot")
# add more features on ax2
# plot Alaska, Hawaii as inset maps
for index,state in usa_more.iterrows():
if state[state_name] in ("Alaska", "Hawaii"):
st_name = state[state_name]
# set fill color, using normalized `sclass` on `my_colormap`
facecolor = my_colormap( state["sclass"] / max(newusa["sclass"] ))
if st_name == "Alaska":
# (1) Alaska
# Custom extent, relative size
map_extent = (-178, -135, 46, 73) # degrees: (lonmin,lonmax,latmin,latmax)
axes_extent = (0.04, 0.06, 0.29, 0.275) # axes units: 0 to 1, (LLx,LLy,width,height)
if st_name == "Hawaii":
# (2) Hawaii
# Custom extent, relative size
map_extent = (-162, -152, 15, 25)
axes_extent = (0.27, 0.06, 0.15, 0.15)
# add inset maps
add_insetmap(axes_extent, map_extent, st_name, \
facecolor, \
edgecolor, \
state["geometry"])
plt.show()
Data classification:
In the code above, sclass carries the classification details of the population density. To get the list of classes, run the code below:
low_val = 0
max_cval = len(sclass.bins)-1
print("Population density (per square_km)")
print("Class Value_ranges")
for ix,val in enumerate(sclass.bins):
#print(low_val, "< x <=", val)
print("{ix:} {low_val:.2f} < density <= {val:.2f}".format(ix=ix, low_val=low_val, val=val))
#print(my_colormap(ix/max_cval))
low_val = val
and you should get the result below that should be accompanied with the map above.
Population density (per square_km)
Class Value_ranges
0 0.00 < density <= 7.48
1 7.48 < density <= 18.56
2 18.56 < density <= 30.50
3 30.50 < density <= 51.70
4 51.70 < density <= 75.38
5 75.38 < density <= 155.09
6 155.09 < density <= 4251.00
Extra colorbar manipulation
It is possible to manipulate the colorbar to change the label texts from simple class number 0,1,2,...6 to the classified value ranges. Just insert these lines of code before the line plt.show() of the code above and rerun.
class_txts = []
low_val = 0
for ix,val in enumerate(sclass.bins):
class_txts.append("{low_val:.1f}, {val:.1f}".format(ix=ix, low_val=low_val, val=val))
low_val = val
cb.ax.set_yticklabels( class_txts )
cb.set_label("Pop-density range")
New plot will look like this:
