How to draw multiple GeoDataFrames (overlapping) in one map figure?

Question

I'm trying to draw a few different layers onto a single map in GeoPandas, but only the last one is ever displayed.

Now, my eventual goal is draw some raster-like polygon 'pixels' based on another numeric column, and I want that overlaid over a map of my area of interest.

Here's what an example with some large, manually created squares:

# modules you'd need to run this
import geopandas as gpd
from matplotlib import pyplot as plt
from shapely.geometry import polygon
%matplotlib inline

# a few example squares
square1 = {
    'geometry':polygon.Polygon([
        (-4164911.2311, 2834480.454299999),
        (-3514002.14019091, 2834480.454299999),
        (-3514002.14019091, 3394480.454299999),
        (-4164911.2311, 3394480.454299999),
        (-4164911.2311, 2834480.454299999)
    ]),
    'count':95
}
square2 = {
    'geometry':polygon.Polygon([
        (-4164911.2311, 3394480.454299999), 
        (-3514002.14019091, 3394480.454299999), 
        (-3514002.14019091, 3954480.454299999), 
        (-4164911.2311, 3954480.454299999), 
        (-4164911.2311, 3394480.454299999)
    ]),
    'count':65
}
square3 = {
    'geometry': polygon.Polygon([
        (-4164911.2311, 3954480.454299999),
        (-3514002.14019091, 3954480.454299999),
        (-3514002.14019091, 4514480.454299999),
        (-4164911.2311, 4514480.454299999),
        (-4164911.2311, 3954480.454299999)
    ]),
    'count':0
}

# squares put into a GeoDataFrame
squares = gpd.GeoDataFrame([square1,square2,square3],
crs={"proj-string: +proj=aea"})

# world country outlines, included with GeoPandas library
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# selecting only for the western hemisphere
westhem = world[(world['continent'] == 'North America') | 
                (world['continent'] == 'South America')]

# ensuring crs info aligns; squares relate to data written onto 
# Albers Equal Area projection
westhem.crs = {"proj-string: +proj=aea"}

# somewhere I read that establishing your pyplot objects is always a 
# good idea
fig,ax = plt.subplots(1,1,sharex=True,sharey=True,figsize=(11,11))

# this should draw some black outlines for countries
westhem.plot(ax=ax,color='white', edgecolor='black');
# and this should map my squares onto the SAME axis.
squares.plot(ax=ax,cmap='Reds',column='count',legend=True);

Instead of drawing both maps, however, I just get this:

How do I make sure both maps appear together so that my squares have context?

Answer 1

You need appropriate CRS conversion to bring two data sources into a common reference system. In the code below, the naturalearth_lowres basemap data is converted to AlbersEqualArea for plotting. Some parameters of AlbersEqualArea I use might be wrong, but they can be corrected and rerun to get the output you need.

import geopandas as gpd
from matplotlib import pyplot as plt
from shapely.geometry import polygon
from cartopy import crs as ccrs

# proj4str  = ccrs.AlbersEqualArea().proj4_init
# '+ellps=WGS84 +proj=aea +lon_0=0.0 +lat_0=0.0 +x_0=0.0 +y_0=0.0 +lat_1=20.0 +lat_2=50.0 +no_defs'

# a few example squares
square1 = {
    'geometry':polygon.Polygon([
        (-4164911.2311, 2834480.454299999),
        (-3514002.14019091, 2834480.454299999),
        (-3514002.14019091, 3394480.454299999),
        (-4164911.2311, 3394480.454299999),
        (-4164911.2311, 2834480.454299999)
    ]),
    'count':95
}
square2 = {
    'geometry':polygon.Polygon([
        (-4164911.2311, 3394480.454299999), 
        (-3514002.14019091, 3394480.454299999), 
        (-3514002.14019091, 3954480.454299999), 
        (-4164911.2311, 3954480.454299999), 
        (-4164911.2311, 3394480.454299999)
    ]),
    'count':65
}
square3 = {
    'geometry': polygon.Polygon([
        (-4164911.2311, 3954480.454299999),
        (-3514002.14019091, 3954480.454299999),
        (-3514002.14019091, 4514480.454299999),
        (-4164911.2311, 4514480.454299999),
        (-4164911.2311, 3954480.454299999)
    ]),
    'count':0
}

# squares put into a GeoDataFrame
squares = gpd.GeoDataFrame([square1,square2,square3],
                           crs={"proj-string: +proj=aea"})

# world country outlines, included with GeoPandas library
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# selecting only for the western hemisphere
westhem = world[(world['continent'] == 'North America') | 
                (world['continent'] == 'South America')]

# prep AlbersEqualArea projection
# default parameters:
#  central_longitude=0.0  <-- not good for USA
#  standard_parallels=(20.0, 50.0)

# for demo purposes: parameters used here may not match the ..
#  intended CRS you want
crs_aea = ccrs.AlbersEqualArea(central_longitude=0, \
                               standard_parallels=(20, 50))

# crs_aea = ...
# '+ellps=WGS84 +proj=aea +lon_0=-75 +lat_0=0.0 +x_0=0.0 +y_0=0.0 +lat_1=20 +lat_2=50 +no_defs'

# convert the geo dataframe to the projection (crs_aea) just created
westhem_crs_aea = westhem.to_crs(crs_aea.proj4_init)

fig,ax = plt.subplots(1,1,sharex=True,sharey=True,figsize=(11,11))

# this should draw some black outlines for countries
westhem_crs_aea.plot(ax=ax, color='lightgray', edgecolor='black');

# and this should map my squares onto the SAME axis
squares.plot(ax=ax,cmap='Reds',column='count',legend=True);

The resulting plot:

Hope this is useful.