import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.io.img_tiles import Stamen, GoogleWTS, GoogleTiles
import matplotlib.patches as mpatches
import shapely.geometry as sgeom

tiler = Stamen(style="terrain-background", cache="/tmp/tiles/")

fig = plt.figure(figsize=(12,7), dpi=120)

proj1 = ccrs.Orthographic(0, 45)
proj2 = ccrs.Orthographic(0, 45)
proj2.threshold /= 20
titles = [f"threshold = {p.threshold:,.2f}" for p in (proj1, proj2)]
for p, proj, title in zip(range(2), [proj1, proj2], titles):

    ax_pla = fig.add_subplot(2, 2, 1+p, projection=ccrs.PlateCarree(), frameon=False)
    ax_ort = fig.add_axes([0.5*p, 0.1, 0.5, 0.5], projection=proj)

    bounds_proj = sgeom.Polygon(proj.boundary)
    geom = ax_ort._get_extent_geom(tiler.crs)
    img, extent, origin = tiler.image_for_domain(geom, 4)

    for ax in [ax_pla, ax_ort]:
        ax.imshow(img, extent=extent, origin=origin, transform=tiler.crs, regrid_shape=500)
        c_geom = ax.add_geometries([geom], crs=tiler.crs, color="red", alpha=0.3, zorder=30)
        c_bounds = ax.add_geometries([bounds_proj], crs=proj, color="green", alpha=0.3, zorder=40)
    ax_pla.set_title(title)
    
geo2 = mpatches.Patch(color='green', label="ortho boundaries")
geo1 = mpatches.Patch(color='red', label='eroded boundaries')
ax_pla.legend(handles=[geo1, geo2], loc=(-0.3,-0.5))