import geopandas as gpd
import matplotlib.pyplot as plt
import cartopy
import contextily as ctx
import pandas as pd
from shapely.geometry.polygon import LinearRing

gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'

## cartopy format
merc_crs = cartopy.crs.Mercator()
latlon_crs = cartopy.crs.PlateCarree()

KML File

routes = gpd.read_file("data/Routes.kml", driver="KML")

Plots

USGS Tile servers are listed here: https://viewer.nationalmap.gov/services/

land_color = '#eeeeee'
coast_color = '#eeeeee'

land_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',
                                               edgecolor=coast_color,
                                               facecolor=land_color)

fig = plt.figure(figsize=(40,40))
ax = fig.add_axes([0,0,1,1], projection=merc_crs)
#ax = fig.add_axes([0,0,1,1], projection=merc_crs

xmin, xmax, ymin, ymax = -149.5, -144, 68.6, 70.6

extent = [xmin, xmax, ymin, ymax]
ax.set_extent(extent, crs=latlon_crs)

colors = ["#A1E2E6", "#E6BDA1", "#B3A16B", "#678072", "#524A4A"]
routes.to_crs(epsg=3857).plot(ax=ax, color="yellow", edgecolor="white", linewidth=5, alpha=1)
#ctx.add_basemap(ax, url=ctx.providers.Stamen.TonerLite)
#ctx.add_basemap(ax, source=ctx.providers.Esri.WorldTopoMap)
#ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery)
#ctx.add_basemap(ax, source=ctx.providers.Esri.NatGeoWorldMap)

#ctx.add_basemap(ax, source=ctx.providers.OpenTopoMap)
#ctx.add_basemap(ax, source='https://a.tile.opentopomap.org/{z}/{x}/{y}.png')

#"USGSShadedReliefOnly"
#"USGSTopo"
#"USGSImageryTopo"
ctx.add_basemap(ax, source="https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryTopo/MapServer/tile/{z}/{y}/{x}")

ax.axis("off")

## Context plot
# ax_sub = fig.add_axes([0.55, 0.8, 0.2, 0.2], projection=cartopy.crs.Mollweide(central_longitude=180))
ax_sub = fig.add_axes(
    [0.79, 0.81, 0.18, 0.18],
    projection=cartopy.crs.NearsidePerspective(central_longitude=-145, central_latitude=55)
)
ax_sub.add_feature(land_10m)
ax_sub.add_geometries([LinearRing([(xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin)])], crs=latlon_crs,
                      color="black")
ax_sub.spines['geo'].set_visible(False)

fig.savefig("output_imagery.png")
plt.show()