Stamen terrain map not working in Cartopy

Question

I am using Cartopy (v0.21.1) to make an inset map in matplotlib (v3.7.1) with cartopy.io.img_tiles.Stamen('terrain-background'). This method worked seamlessly up until about a week ago.

Now I get this error that only occurs sometimes:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/anaconda3/envs/xuscd/lib/python3.8/site-packages/IPython/core/formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/IPython/core/pylabtools.py in <lambda>(fig)
    246 
    247     if 'png' in formats:
--> 248         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    249     if 'retina' in formats or 'png2x' in formats:
    250         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, **kwargs)
    130         FigureCanvasBase(fig)
    131 
--> 132     fig.canvas.print_figure(bytes_io, **kw)
    133     data = bytes_io.getvalue()
    134     if fmt == 'svg':

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2340                 )
   2341                 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2342                     self.figure.draw(renderer)
   2343 
   2344             if bbox_inches:

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     93     @wraps(draw)
     94     def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95         result = draw(artist, renderer, *args, **kwargs)
     96         if renderer._rasterizing:
     97             renderer.stop_rasterizing()

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer)
     70                 renderer.start_filter()
     71 
---> 72             return draw(artist, renderer)
     73         finally:
     74             if artist.get_agg_filter() is not None:

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/matplotlib/figure.py in draw(self, renderer)
   3138 
   3139             self.patch.draw(renderer)
-> 3140             mimage._draw_list_compositing_images(
   3141                 renderer, self, artists, self.suppressComposite)
   3142 

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    129     if not_composite or not has_images:
    130         for a in artists:
--> 131             a.draw(renderer)
    132     else:
    133         # Composite any adjacent images together

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer)
     70                 renderer.start_filter()
     71 
---> 72             return draw(artist, renderer)
     73         finally:
     74             if artist.get_agg_filter() is not None:

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py in draw(self, renderer, **kwargs)
    529         if not self._done_img_factory:
    530             for factory, factory_args, factory_kwargs in self.img_factories:
--> 531                 img, extent, origin = factory.image_for_domain(
    532                     self._get_extent_geom(factory.crs), factory_args[0])
    533                 self.imshow(img, extent=extent, origin=origin,

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/cartopy/io/img_tiles.py in image_for_domain(self, target_domain, target_z)
     97 
     98         print(f'img shape = {img.shape}')
---> 99         img, extent, origin = _merge_tiles(tiles)
    100         return img, extent, origin
    101 

~/anaconda3/envs/xuscd/lib/python3.8/site-packages/cartopy/io/img_tiles.py in _merge_tiles(tiles)
    652         print(img_slice)
    653 
--> 654         img[img_slice] = tile_img
    655 
    656     return img, [min(xs), max(xs), min(ys), max(ys)], 'lower'

ValueError: could not broadcast input array from shape (512,512,3) into shape (766,766,3)

Did the Stamen map update? Could this be breaking the existing Cartopy code by modifying the size of the map image?

Answer 1

Stamen stopped providing a base map service last week - this was announced some time ago, their maps are now provided by Stadia.

Answer 2

As described in other answers, Stamen has stopped providing a base map service. As a workaround to continue using your Cartopy code (almost) as-is:

1. Create a free account (200k credits/month) and get an API key from Stadia Maps here: https://stadiamaps.com/stamen/onboarding/create-account/

2. Override the Cartopy Stamen map class as follows, using your API key:

import cartopy.io.img_tiles as cimgt

class StadiaStamen(cimgt.Stamen):
    def _image_url(self, tile):
         x,y,z = tile
         url = f"https://tiles.stadiamaps.com/tiles/stamen_terrain_background/{z}/{x}/{y}.png?api_key={API_KEY}"
         return url

3. Replace calls to the old Stamen class with the new class, e.g.

stamen_terrain = StadiaStamen("terrain-background")