Unable to retrieve map tiles at high zoom level (i.e. zoomed out) with Matplotlib/Cartopy/Python

Question

I am creating a mapping application with Python and Cartopy, and attempting to use open-source map tiles for backgrounds, to have more options than the default Cartopy map.

It works perfectly for maps that are zoomed in fairly close, but when I try to get a view from a higher elevation, something fails. If I have the zoom set to 11, it works. If I set it to 12, it hangs indefinitely and gives no traceback.

Same result with both the OSM and Stamen map servers.

Here is a short, self-contained example (note that a line or two may be artifacts from the various ways I've tried this)

import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt


def custom_background(source_point):

    source_point = source_point.split(" ")
    source_point = (float(source_point[0]), float(source_point[1]))
    dx = 1.5
    dy = 1.5
    lon_min, lon_max = source_point[0]-dx, source_point[0]+dx
    lat_min, lat_max = source_point[1]-dy, source_point[1]+dy
    zoom = 7
    map_url = "https://www.openstreetmap.org/#map={}/{}/{}".format(zoom,source_point[0],source_point[1])
    tile = cimgt.OSM(url=map_url)
    tile = cimgt.StamenTerrain()
    ax = plt.axes(projection=ccrs.PlateCarree())
    ax.set_extent([lat_min, lat_max, lon_min, lon_max])
    ax.add_image(tile, zoom)
    #~ ax.add_image(tile)
    return ax

custom_background("45.068466 -66.45477")
plt.savefig("tile.png")

the result, with zoom = 7:

but if I were to change zoom to say, 14, the program will not complete no matter how long I allow it to run.

The url parameter that is passed to cimgt.OSM() is optional. I get the same result with or without it. (See: https://scitools.org.uk/cartopy/docs/v0.16/cartopy/io/img_tiles.html#cartopy.io.img_tiles.OSM)

Am I missing something here? Any help would be appreciated, thank you.

Answer 1

The "zoom" levels are based on a Quadtree. Essentially increasing the "zoom" resolution increases the number of tiles by a power of four.

So:

zoom level 0: 4^0 = 1 tile(s) to cover the globe
zoom level 1: 4^1 = 4 tile(s) to cover the globe
...
zoom level 7: 4^7 = 16,384 tile(s) to cover the globe
...
zoom level 14: 4^14 = 268,435,456 tile(s) to cover the globe

So if you request tiles at a high zoom-level for a large area you can end up requesting a lot of tiles.

There is a useful but undocumented method on the Tiler object called find_images. Its implementation isn't too hairy:
https://github.com/SciTools/cartopy/blob/v0.16.0/lib/cartopy/io/img_tiles.py#L103-L122.

With this method, we can actually see the tiles that would be used for a given range. Importantly the range needs to be supplied in the coordinate system of the tiler (almost exclusively Web Mercator).

In [1]: import cartopy.io.img_tiles as cimgt

In [2]: import shapely.geometry as sgeom

In [3]: import cartopy.crs as ccrs

In [4]: tiler = cimgt.OSM()

In [5]: pt = 45.068466, -66.45477

In [6]: target = sgeom.box(pt[0] - 1.5, pt[1] - 1.5, pt[0] + 1.5, pt[1] + 1.5)

In [7]: target_mercator = tiler.crs.project_geometry(target, ccrs.Geodetic()).geoms[0]

So with all the pieces in place we can start finding out which tiles we would need to draw at particular zoom levels:

For the target you are specifying, at zoom level 0, we need the following tiles (x, y, z):

In [8]: list(tiler.find_images(target_mercator, 0))
Out[8]: [(0, 0, 0)]

For z=1, it is still only one tile:

In [9]: list(tiler.find_images(target_mercator, 1))
Out[9]: [(1, 1, 1)]

But for z=2 we clearly crossed a tile boundary, as we now need two tiles to cover the target domain:

In [10]: list(tiler.find_images(target_mercator, 2))
Out[10]: [(2, 2, 2), (2, 3, 2)]

Naturally, the list grows as we increase the zoom level:

In [11]: list(tiler.find_images(target_mercator, 3))
Out[11]: [(4, 5, 3), (5, 5, 3), (4, 6, 3), (5, 6, 3)]

In [12]: list(tiler.find_images(target_mercator, 6))
Out[12]: [(39, 47, 6), (40, 47, 6), (39, 48, 6), (40, 48, 6)]

When we hit z=7 we find that we would need 8 tiles to represent the area in question:

In [13]: list(tiler.find_images(target_mercator, 7))
Out[13]: 
[(79, 94, 7),
 (79, 95, 7),
 (80, 94, 7),
 (80, 95, 7),
 (79, 96, 7),
 (79, 97, 7),
 (80, 96, 7),
 (80, 97, 7)]

I'm sure you can see where this is going, but let's try your zoom level of 14 to find out how many tiles we would need. To save on electricity and reading time, let's just print the length of that list...

In [14]: len(list(tiler.find_images(target_mercator, 14)))
Out[14]: 47334

Right, so requesting zoom level 14 at your desired extent would need you to download ~47,000 tiles at 256x256 (8-bit colormapped PNGs). The z=0 tile (https://a.tile.openstreetmap.org/0/0/0.png) is about 8882 bytes compressed, so assuming this is typical, you'd end up downloading ~420420588 bytes (400MB). In order to hold this data in memory, you'd also need about 2.9GB of RAM. Finally, to reproject this data to PlateCarree you'd need to at least double this amount of RAM, and that is assuming a highly efficient reprojection implementation (cartopy's isn't).

Hopefully this tells you why your code seems to take a long time - you are asking it to do a lot of work. There have been a few discussions about whether cartopy should warn when an excessively large number of tiles is requested, but it always comes down to what is an unreasonable number (you might actually want to fetch that many tiles!). We have also talked about automatic zoom level selection - this is something that is really feasible if there is sufficient demand.

HTH