ValueError: could not broadcast input array from shape (256,256,4) into shape (512,512,4)

[galenseilis]

I'm at a loss to what has changed, but this Jupyter notebook has stopped working. Troubleshooting feedback is welcome. The issue seems to be a mismatch in array shape. Not sure why, but one array appears to be twice the size of the other.

import time
import warnings

warnings.filterwarnings("ignore")

import arviz as az
import contextily as cx
from IPython.display import YouTubeVideo
import graphviz
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
from scipy import stats
import seaborn as sns
import xarray as xr
from watermark import watermark

%load_ext watermark

# Details about the system running the report and when the report was generated.
%watermark

Here is the watermark:

Last updated: 2023-10-03T15:11:30.607964-07:00

Python implementation: CPython
Python version       : 3.10.4
IPython version      : 8.14.0

Compiler    : MSC v.1929 64 bit (AMD64)
OS          : Windows
Release     : 10
Machine     : AMD64
Processor   : Intel64 Family 6 Model 85 Stepping 0, GenuineIntel
CPU cores   : 8
Architecture: 64bit

The is the block in the Jupyter notebook which is breaking.

# Plot Annual Exam Volume
axes = merged_df.plot(
    "annual_exam_volume",
    cmap="cool",
    scheme="NaturalBreaks",
    missing_kwds={
        "color": "lightgrey",
        "edgecolor": "red",
        "hatch": "///",
        "label": "Missing values",
    },
    legend=True,
    legend_kwds={"fmt": "{:.0f}", "interval": True},
    figsize=(25, 15),
    alpha=0.7,
)

plt.axis("off")

axes.set_title("Annual MRI Breast Cancer Exam Volume By LHA")

cx.add_basemap(axes, crs=merged_df.crs)

merged_df.apply(
    lambda x: axes.annotate(
        text=x["LHA"], xy=x.lha_boundary.centroid.coords[0], ha="left"
    ),
    axis=1,
);

Here is the traceback:

in <module>:21                                                                                   │
│                                                                                                  │
│   18                                                                                             │
│   19 axes.set_title("Annual MRI Breast Cancer Exam Volume By LHA")                               │
│   20                                                                                             │
│ ❱ 21 cx.add_basemap(axes, crs=merged_df.crs)                                                     │
│   22                                                                                             │
│   23 merged_df.apply(                                                                            │
│   24 │   lambda x: axes.annotate(                                                                │
│                                                                                                  │
│ E:\Breast_Cancer_MRI_Screen_Demand_Forecast\venv\lib\site-packages\contextily\plotting.py:129 in │
│ add_basemap                                                                                      │
│                                                                                                  │
│   126 │   │   │   │   left, right, bottom, top, crs, {"init": "epsg:3857"}                       │
│   127 │   │   │   )                                                                              │
│   128 │   │   # Download image                                                                   │
│ ❱ 129 │   │   image, extent = bounds2img(                                                        │
│   130 │   │   │   left, bottom, right, top, zoom=zoom, source=source, ll=False                   │
│   131 │   │   )                                                                                  │
│   132 │   │   # Warping                                                                          │
│                                                                                                  │
│ E:\Breast_Cancer_MRI_Screen_Demand_Forecast\venv\lib\site-packages\contextily\tile.py:225 in     │
│ bounds2img                                                                                       │
│                                                                                                  │
│   222 │   │   image = _fetch_tile(tile_url, wait, max_retries)                                   │
│   223 │   │   tiles.append(t)                                                                    │
│   224 │   │   arrays.append(image)                                                               │
│ ❱ 225 │   merged, extent = _merge_tiles(tiles, arrays)                                           │
│   226 │   # lon/lat extent --> Spheric Mercator                                                  │
│   227 │   west, south, east, north = extent                                                      │
│   228 │   left, bottom = mt.xy(west, south)                                                      │
│                                                                                                  │
│ E:\Breast_Cancer_MRI_Screen_Demand_Forecast\venv\lib\site-packages\contextily\tile.py:628 in     │
│ _merge_tiles                                                                                     │
│                                                                                                  │
│   625 │                                                                                          │
│   626 │   for ind, arr in zip(indices, arrays):                                                  │
│   627 │   │   x, y = ind                                                                         │
│ ❱ 628 │   │   img[y * h : (y + 1) * h, x * w : (x + 1) * w, :] = arr                             │
│   629 │                                                                                          │
│   630 │   bounds = np.array([mt.bounds(t) for t in tiles])                                       │
│   631 │   west, south, east, north = (  

[chudlerk]

Having the same issue outside of Jupyter notebook (using Spyder, personally). Same issue pops up even with the most basic example from the documentation.

import contextily as cx
import geopandas

data_url = "https://ndownloader.figshare.com/files/20232174"

db = geopandas.read_file(data_url)

ax = db.plot(color="red", figsize=(9, 9))
cx.add_basemap(ax, crs=db.crs.to_string())

  File C:\ProgramData\Miniconda3\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\users\kchudler\onedrive - research triangle institute\documents\hydromet_tools\storm_analysis_system\untitled1.py:18
    cx.add_basemap(ax, crs=db.crs.to_string())

  File C:\ProgramData\Miniconda3\lib\site-packages\contextily\plotting.py:129 in add_basemap
    image, extent = bounds2img(

  File C:\ProgramData\Miniconda3\lib\site-packages\contextily\tile.py:225 in bounds2img
    merged, extent = _merge_tiles(tiles, arrays)

  File C:\ProgramData\Miniconda3\lib\site-packages\contextily\tile.py:628 in _merge_tiles
    img[y * h : (y + 1) * h, x * w : (x + 1) * w, :] = arr

ValueError: could not broadcast input array from shape (512,512,4) into shape (256,256,4)

[dmaldona]

I am having the same issue.

[martinfleis]

This seems to be an issue with Stamen Terrain tiles. The code from the @chudlerk's post fetches 9 tiles, but their shape are not consistent

(256, 256, 4), (256, 256, 4), (256, 256, 4), (256, 256, 4), (256, 256, 4), (256, 256, 4), (256, 256, 4), (256, 256, 4), (512, 512, 4)

See that the last one has double the resolution. I wonder if that is intention or some artifact of the transition to Stadia tile service. Maybe @ianthetechie will know?

In the meantime, you can use different tiles which do not seem to have the same issue, like cx.providers.OpenStreetMap.HOT or any other.

[jorisvandenbossche]

Can you try with a different provider source? For example cx.add_basemap(ax, source=cx.providers.OpenStreetMap.Mapnik)

I think this error is due because the default of Stamen is going to be removed (see https://github.com/geopandas/contextily/issues/220), and it seems that they now changed the tiles to have a message, like:

for https://stamen-tiles-a.a.ssl.fastly.net/terrain/6/28/23.png. And this image has a size of 512x512 instead of 256x256 like the normal tiles.

[jorisvandenbossche]

I suppose they give you such a tile once every few or with some logic, so that we see only some of the tiles have this other size.

In any case, time to do a release so we have the changed default out?

[martinfleis]

In any case, time to do a release so we have the changed default out?

Ideally with #222 in and a fix for tests (not sure how, yet, probably by relaxing the checks or depending on more stable tiles).

[ianthetechie]

Hey, thanks for the summons @martinfleis.

See that the last one has double the resolution. I wonder if that is intention or some artifact of the transition to Stadia tile service.

That'd be a bit of an accident. We currently don't discriminate on the 1x/2x resolution on these error tiles since most maps display them correctly. As @jorisvandenbossche noticed, we started serving these "brownout" warning tiles yesterday for every 16th tile. We will increase the frequency as the month goes on. If you access the tiles using a Stadia Maps API key, you'll get the correct tiles.

Stamen announced back in July (https://maps.stamen.com/stadia-partnership/) that they were shutting down their Fastly CDN as it became too expensive. You will need an API key to continue access Stamen tiles. They remain free for non-commercial use, but commercial use and higher volumes require a paid plan. You can sign up for a key here: https://stadiamaps.com/stamen/onboarding/create-account.

[martinfleis]

Thanks! I guess that we just cut a release with the new default and point users to new Stadia service if they want to keep using Stamen tiles.

[ianthetechie]

Sounds good. You know where to shout if you need some help with that (or just a quick once-over). You can find the updated URLs here: https://docs.stadiamaps.com/guides/migrating-from-stamen-map-tiles/#switch-your-tile-urls.

[galenseilis]

Can you try with a different provider source? For example cx.add_basemap(ax, source=cx.providers.OpenStreetMap.Mapnik)

I think this error is due because the default of Stamen is going to be removed (see #220), and it seems that they now changed the tiles to have a message, like:

for https://stamen-tiles-a.a.ssl.fastly.net/terrain/6/28/23.png. And this image has a size of 512x512 instead of 256x256 like the normal tiles.

Yes, using a different provider source provided an acceptable workaround.

[galenseilis]

Thanks! I guess that we just cut a release with the new default and point users to new Stadia service if they want to keep using Stamen tiles.

Sounds good. Thanks.

[peanutfun]

I still experience the original issue when switching to the new Stamen URLs at the Stadia service:

>>> contextily.add_basemap(plt.gca(), source="https://tiles.stadiamaps.com/tiles/stamen_toner_lite/{z}/{x}/{y}{r}.png")

ValueError: could not broadcast input array from shape (256,256,4) into shape (512,512,4)

It seems like the resolution of the tiles changed and contextily cannot handle that. Is there any workaround? I would like to keep using the Stamen maps with contextily.

[jorisvandenbossche]

It seems that new link is still giving the same "warning" tiles every x requests (eg I tried https://tiles.stadiamaps.com/tiles/stamen_toner_lite/3/1/1.png several times). Maybe that's because you didn't provide an API key? (if that's the case)

[peanutfun]

Thank you, I did not realize I needed to register! Worked fine once I added the API key to the URL 🙌

For reference:

• https://tiles.stadiamaps.com/tiles/stamen_toner/{z}/{x}/{y}{r}.png gives the ValueError
• https://tiles.stadiamaps.com/tiles/stamen_toner/{z}/{x}/{y}.png gives the ValueError
• https://tiles.stadiamaps.com/tiles/stamen_toner/{z}/{x}/{y}@2x.png gives a plot with the warning tiles

Adding an API key from the Stadia Maps registration to any of these URLs provided the map as expected.

[ianthetechie]

@peanutfun Good to gear you got it sorted :) For anyone finding this issue wanting to know the full resolution, here are some details that I think should help.

We have been basically letting everything slide through since the deprecation announcement in July, but these tiles DO require an API key going forward. I'm not up on whether contextily has a way to manage API keys, but you can add it as a query string parameter like ?api_key=YOUR-API-KEY-HERE to your URLs.

@jorisvandenbossche has the right idea. Stamen never collected any information on who was using their tiles, so there was no way to contact everyone using the old Fastly URLs. The solution we settled on was a "warning" tile. Since ~1 week ago, all unauthenticated requests will now get these warning tiles pseudo-randomly. Currently it's ~every 16th tile on the map, and will increase as the old URL shutoff date approaches.

[jorisvandenbossche]

@ianthetechie would it be possible to generate a version of the warning tile at 256x256 ?

[ianthetechie]

Yeah, that's a good question @jorisvandenbossche. I'll look into that. Our current tile is 512 because all the major web renderers handle it fine, but we probably should offer a 256 version as well.

[martinfleis]

Because Stadia allows both domain whitelisting and API key authentication, passing the key requires a bit more than usual. See the snippet below.

import geopandas
import geodatasets
import contextily

nybb = geopandas.read_file(geodatasets.get_path("nybb")).to_crs(3857)

# get the provider object to adapt
tiles = contextily.providers.Stadia.StamenToner
tiles["url"] = tiles["url"] + "?api_key={api_key}"  # adding API key placeholder

ax = nybb.plot(alpha=.5)
contextily.add_basemap(ax=ax, source=tiles(api_key="my_key_here"))