[Bug]: Issue overlaying geographic vector data on geographic raster data

[Robinlovelace]

Bug summary

I have 2 datasets in the same CRS. They plot fine individually. But when you try to plot one on top of another it fails.

Code for reproduction

import matplotlib as mpl
from matplotlib import pyplot
import geopandas as gpd
import rasterio
import rasterio.plot
nz_elev = rasterio.open('https://github.com/geocompr/py/blob/main/data/nz_elev.tif?raw=true')
nz_elev_read = nz_elev.read(1)
nz_transformed = gpd.read_file('https://github.com/geocompr/py/releases/download/0.1/nz_transformed.gpkg')

nz_transformed.plot(facecolor='none', edgecolor='r'); # works

rasterio.plot.show(nz_elev_read);

fig, ax = pyplot.subplots(figsize=(15, 15))
rasterio.plot.show(nz_elev_read, ax=ax);
nz_transformed.plot(ax=ax, facecolor='none', edgecolor='r'); # fails

Actual outcome

Expected outcome

The overlay between the two images below

Additional information

Illustration of this issue built on reproducible CI on GitHub actions: https://geocompr.github.io/py/09-mapping.html

Operating system

Ubuntu

Matplotlib Version

3.1.2

Matplotlib Backend

module://matplotlib_inline.backend_inline

Python version

3.8.10 and 3.9

Jupyter version

?

Installation

conda

[Robinlovelace]

A bit more context and working example that may be of use. First question I have is: can anyone else reproduce this behaviour?

Here's an issue that puts in context of open source teaching materials: geocompx/geocompy#83
And here's the perhaps not perfect demo of this approach working: https://gis.stackexchange.com/questions/294072/how-can-i-superimpose-a-geopandas-dataframe-on-a-raster-plot

[oscargus]

I can reproduce it.

Although I do not have any insights into the problem, I can confirm that the nz_transformed call is adding a collection and the rasterio call is adding an image to ax. The same thing happens to the first two calls (pyplot.gca()). (I do only get the image, not the borders, on the first two calls, btw, but that is maybe not what you meant with "working"?)

(I used the qt-backend and current main-branch, btw.)

[Robinlovelace]

I can reproduce it.

Thanks for reproducing it. Just to clarify, do you get the upside-down New Zealand result?

[oscargus]

Ohh, it is actually doing the "right" thing. The image has a range of about 1000, while the outline has a range of about 1000000, with an offset. Also, the yaxis is inverted somewhere.

So zoom in to the top left corner and the image will appear.

[Rabscuttler]

Pass the rasterio object transform as an argument to the rasterio.plot.show:
transform=nz_elev.transform

From rasterio docs.

Clue is originally the axis labels are going from 0 to (max image height) - ie plotting as pixels rather than the desired coordinates, so need the transform.

[Robinlovelace]

Many thanks @Rabscuttler this is awesome! Demo of the fix:

# Plot combined
fig, ax = pyplot.subplots(figsize=(15, 15))
rasterio.plot.show(nz_elev_read, ax=ax, transform=nz_elev.transform);
nz_transformed.plot(ax=ax, facecolor='none', edgecolor='r'); # fails