Unable to get the correct tiles using `bounds2img` via contextily and geopandas for USA

The main idea is to get a png/jpg image file as a basemap for the complete bounds of a geodataframe.

The geodataframe here happens to be a map of the USA which I took from the topojson/us-atlas repository – this shapefile -> US States

I do some plotting on the map, but that is not important. We will use the simple raw map itself. Ideally I plot it in "albersUsa" projection, but I have been stuck at not getting the correct tiles itself to worry about the projection.

Attempt
Loading the data and assigning it a proper CRS since the original file does not seem to have any information about its CRS and finnaly converting it to the CRS that contextily expects –

import pandas as pd
import contextily as ctx
import matplotlib.pyplot as plt # mostly for debugging

usa_gdf = gpd.read_file('https://cdn.jsdelivr.net/npm/us-atlas@3/states-10m.json')
empty = usa_gdf.geometry.is_empty
usa_gdf = usa_gdf[~empty]
print(usa_gdf.crs)
usa_gdf.crs = 'esri:102003'
print(usa_gdf.crs)
usa_gdf = usa_gdf.to_crs(epsg=4326)  #EPSG:3857 also seems to be okay
print(usa_gdf.crs)

Output:

None
esri:102003
epsg:4326

Plotting the data to see what it looks like (ideally I want to keep everything in albersUsa projection) –

ax = usa_gdf.plot()

Output:

Just for checking, I use the add_basemap method of contextily which seems to work seamlessly with a matplotlib plot as the documentation shows –

ctx.add_basemap(ax, crs=usa_gdf.crs.to_string())

But it gives me the following error –

UserWarning: The inferred zoom level of 20 is not valid for the current tile provider (valid zooms: 0 - 18).
  warnings.warn(msg)

---------------------------------------------------------------------------

HTTPError                                 Traceback (most recent call last)

/usr/local/lib/python3.6/dist-packages/contextily/tile.py in _retryer(tile_url, wait, max_retries)
    446         request = requests.get(tile_url, headers={"user-agent": USER_AGENT})
--> 447         request.raise_for_status()
    448     except requests.HTTPError:

8 frames

HTTPError: 404 Client Error: Not Found for url: https://stamen-tiles-a.a.ssl.fastly.net/terrain/18/61167/101576.png


During handling of the above exception, another exception occurred:

HTTPError                                 Traceback (most recent call last)

/usr/local/lib/python3.6/dist-packages/contextily/tile.py in _retryer(tile_url, wait, max_retries)
    450             raise requests.HTTPError(
    451                 "Tile URL resulted in a 404 error. "
--> 452                 "Double-check your tile url:n{}".format(tile_url)
    453             )
    454         elif request.status_code == 104:

HTTPError: Tile URL resulted in a 404 error. Double-check your tile url:
https://stamen-tiles-a.a.ssl.fastly.net/terrain/18/61167/101576.png

When I provide a source I get a text output with warning but no image –

ctx.add_basemap(ax, crs=usa_gdf.crs.to_string(), source=ctx.providers.Stamen.Toner)

Output:

UserWarning: The inferred zoom level of 20 is not valid for the current tile provider (valid zooms: 0 - 19).
  warnings.warn(msg)

<Figure size 432x288 with 0 Axes>

Then I try some lower zoom numbers zoom=3 with a few tile providers but everytime I just get a text output, no map or the previous error – the errors are themselves not consistent –

<Figure size 432x288 with 0 Axes>

Getting the tiles using inbuilt methods
To get the basemap image I looked at the documentation of contextily and found that bound2img method may help me get what I want so I tried it this way –

bounds = usa_gdf.total_bounds
us_img, us_ext = ctx.bounds2img(bounds[0],bounds[1],bounds[2],bounds[3],ll=True, zoom=4, source=ctx.providers.CartoDB.Positron)

This will output us_img as an RGB array. So we will convert and store it into a jpg image and then read it using PIL OR plot it directly using Matplotlib. Since both give the exact same result, I will show the PIL method here –

from PIL import Image
im = Image.fromarray(us_img)
im.save("us.jpeg")

Reading the file and plotting it –

f, ax = plt.subplots(1)
ax.imshow(plt.imread('us.jpeg'), extent=us_ext)

But all I get is this –

Playing with different zoom levels just gives me the world or zooms into Kansas…never the proper whole USA.

Any idea what is it that I am not doing right? And how to see the map of USA in albersUsa projection along with a basemap?

NOTE: It is important for me to get the tiles without using Matplotlib, as I said, here I used Matplotlib just for debugging. Ideally i’d just use PIL, contextily and geopandas. Because of this constraint I decided to use total_bounds and bounds2img , if there are other ways to get this done then I am open to those too. Also, if I work on a different geodataframe then I’d like to have the basemap for that full geodataframe, so the solution should take that into consideration too.

Resources:

• Intro to Contextily
• Data/CRS/Projections that contextily works well with
• Saving tiles as images and working locally