Contextily, how to get max zoom allowable and set map extents accordingly

Question

Using Python 3.8, GeoPandas, and Contextily, I am plotting a lot of maps in different areas at different zoom levels (looping through list of points in a GeoDataFrame). Different zoom levels work for different global areas.

What is the best way to set the zoom to the max allowable, capturing the points that I am plotting? In the code below, I don't know how to:

1. Return max zoom allowable (I just used 13 in second plot because that was in error message)
2. Change the extents in the second plot-- the extents of the second plot are the same as the first; I would have expected that changing the zoom level would change the extents

I think I have correctly set the EPSG, which was the gist of the answer to a similar question here.

Code:

import pandas as pd
import geopandas as gpd
%matplotlib inline 
import matplotlib.pyplot as plt
import contextily as ctx
from shapely.geometry import Point
plt.style.use('seaborn-whitegrid')

long, lat = [(-118.02, -118.051, -118.04), (39.499, 39.512, 39.501)]
df = pd.DataFrame(list(zip(lat, long)), columns =['lat', 'long'])
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['long'], df['lat']))
gdf.crs = "EPSG:4326"

fig = plt.figure(figsize=(10,7.5), constrained_layout=True)
gs = fig.add_gridspec(1, 2)
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[0, 1])

gdf.plot(ax = ax1)
ctx.add_basemap(ax1, crs='epsg:4326', source=ctx.providers.Esri.WorldShadedRelief)
ax1.tick_params('x', labelrotation=90)
ax1.set_aspect('equal')
ax1.set_title('Auto-Zoom')

gdf.plot(ax = ax2)
ctx.add_basemap(ax2, crs='epsg:4326', source=ctx.providers.Esri.WorldShadedRelief, zoom = 13)
ax2.tick_params('x', labelrotation=90)
ax2.set_aspect('equal')
ax2.set_title('Defined Zoom')

Answer 1

As it is confirmed by the asker. Adding option reset_extent=False to ctx.add_basemap(ax2, ... ) will solve the problem of plotting on ax2.

For checking a possible zoom value before actually using it. Use this

# Test if zoom: 28 is valid or not?
retval = ctx.tile._validate_zoom(28, ctx.providers.Esri.WorldShadedRelief, auto=True)

The value return, here retval represents a valid zoom value (less or equal the check value). In this case, 28 is used to check but 13 is returned. Thus 13 represents the maximum valid zoom for that particular basemap tiles' images.