imshow fails when transform projection is different that the axes projection

[edsaac]

Description

imshow works fine when both axes and raster data projections are the same (Figure 1), but it fails when the axes has a projection that is different from the data (Figure 2). As a workaround, the data extent has to be projected independently to the axes projection and pass that as the extent to imshow (Figure 3)

Figure	Map
1
2
3

To reproduce

Code

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs  
import cartopy.feature as cf
from shapely import Point

map_extent = (-88, -86, 41.2, 42.2)

# Some dummy raster data
raster = np.array([[1,2], [3,4]])  
raster_extent = (-87.6, -86.4, 41.4, 42.0) # (left, right, bottom, top)
raster_proj = ccrs.PlateCarree()

# Corners of the raster data
points = (raster_extent[:2], raster_extent[2:])
points_proj = ccrs.PlateCarree()


#----- Figure 1 ------------
## Map configuration
axes_proj = ccrs.PlateCarree()
fig, ax = plt.subplots(subplot_kw={"projection": axes_proj}, figsize=(6, 6))

## Add some map features
ax.add_feature(cf.LAKES, alpha=0.3)
ax.add_feature(cf.STATES)

## Add my raster data
ax.imshow(
    raster,
    extent=raster_extent,
    transform=raster_proj,
    origin="lower",
    cmap="spring"
)

## Add my point data
ax.scatter(
    *points, 
    transform=ccrs.PlateCarree(),
    s=100, c='k',
    marker="x",
    label="imshow extent"
)

## Final touches
ax.set_extent(map_extent, crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)
ax.set_title("PlateCarree")
ax.legend()
plt.show()


#----- Figure 2 ------------
## Map configuration
axes_proj = ccrs.LambertConformal(central_longitude=-87, central_latitude=41.7)

fig, ax = plt.subplots(subplot_kw={"projection": axes_proj}, figsize=(6, 6))

## Add some map features
ax.add_feature(cf.LAKES, alpha=0.3)
ax.add_feature(cf.STATES)

## Add my raster data
ax.imshow(
    raster,
    extent=raster_extent,
    transform=raster_proj,
    origin="lower",
    cmap="spring"
)

## Add my point data
ax.scatter(
    *points, 
    transform=ccrs.PlateCarree(),
    s=100, c='k',
    marker="x",
    label="imshow extent"
)

## Final touches
ax.set_extent(map_extent, crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)
ax.set_title("LambertConformal")
ax.legend()
plt.show()


#----- Figure 3 ------------
## Map configuration
axes_proj = ccrs.LambertConformal(central_longitude=-87, central_latitude=41.7)


fig, ax = plt.subplots(subplot_kw={"projection": axes_proj}, figsize=(6, 6))

## Add some map features
ax.add_feature(cf.LAKES, alpha=0.3)
ax.add_feature(cf.STATES)

## Add my raster data

### Create shapely points
lower_left = Point(raster_extent[0], raster_extent[2])
upper_right = Point(raster_extent[1], raster_extent[3])

### Project from the original raster projection to the axes projection
projected_lower_left = axes_proj.project_geometry(lower_left, raster_proj)
projected_upper_right = axes_proj.project_geometry(upper_right, raster_proj)

### Assemble to pass to extent
projected_raster_extent = (
    projected_lower_left.x,
    projected_upper_right.x,
    projected_lower_left.y,
    projected_upper_right.y
)

ax.imshow(
    raster,
    extent=projected_raster_extent,
    transform=axes_proj,
    origin="lower",
    cmap="spring"
)

## Add my point data
ax.scatter(
    *points, 
    transform=ccrs.PlateCarree(),
    s=100, c='k',
    marker="x",
    label="imshow extent"
)

## Final touches
ax.set_extent(map_extent, crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)
ax.set_title("LambertConformal [workaround]")
ax.legend()
plt.show()

Environment

Python==3.12.7
Cartopy==0.24.1 
shapely==2.0.6

PD:
#2459 seemed related but I am afraid this is a separate issue.

[edsaac]

Digging a little in geoaxis.py I noticed that the image reprojection is dependent on the plot extent. Setting the map extent before adding the raster data with imshow plots the map correctly (except for the top/bottom pixels that get clipped as described in #2459).

cartopy/lib/cartopy/mpl/geoaxes.py

Lines 1290 to 1305 in 0d4e39f

	target_extent = self.get_extent(self.projection)
	regrid_shape = kwargs.pop('regrid_shape', 750)
	regrid_shape = self._regrid_shape_aspect(regrid_shape,
	target_extent)
	# Lazy import because scipy/pykdtree in img_transform are only
	# optional dependencies
	from cartopy.img_transform import warp_array
	original_extent = extent
	img, extent = warp_array(img,
	source_proj=transform,
	source_extent=original_extent,
	target_proj=self.projection,
	target_res=regrid_shape,
	target_extent=target_extent,
	mask_extrapolated=True,
	)

Below is an example of this, adapting from the Reprojecting Images example in the docs

map extent set after imshow	map extent set before imshow

Code

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf


def dummy_image():
    """
    Return a 2x2 dummy image with a projection and extent.

    Returns
    -------
    img : numpy array
        The pixels of the image in a numpy array.
    img_proj : cartopy CRS
        The rectangular coordinate system of the image.
    img_extent : tuple[float, float, float, float]
        The extent of the image ``(left, right, bottom, top)`` referenced in
        the ``img_proj`` coordinate system.
    origin : str
        The origin of the image to be passed through to matplotlib's imshow.

    """
    img = np.linspace(0, 1, 12).reshape(6, 2)
    img_proj = ccrs.PlateCarree()
    img_extent = (-87.6, -86.4, 41.4, 42.0)
    return img, img_proj, img_extent, "lower"


def buggy_map(axes_proj):
    fig, ax = plt.subplots(subplot_kw={"projection": axes_proj})

    # Add some map features
    ax.add_feature(cf.LAKES, alpha=0.3)
    ax.add_feature(cf.STATES)

    # Add the image before the extent of the map is set
    img, crs, extent, origin = dummy_image()
    ax.imshow(img, transform=crs, extent=extent, origin=origin, cmap="spring")

    ## Add raster corners
    ax.scatter(
        extent[:2],
        extent[2:],
        transform=ccrs.PlateCarree(),
        s=100,
        c="k",
        marker="x",
        label="imshow extent",
    )

    # Set the map extexts after imshow
    map_extent = (-88, -86, 41.2, 42.2)
    ax.set_extent(map_extent, crs=ccrs.PlateCarree())

    # Final touches
    ax.gridlines(draw_labels=True)
    ax.set_title(f"\u2717 Using {type(axes_proj)}")
    ax.legend()
    fig.savefig("buggy_map.png")


def not_buggy_map(axes_proj):
    fig, ax = plt.subplots(subplot_kw={"projection": axes_proj})

    # Set the map extexts before imshow
    map_extent = (-88, -86, 41.2, 42.2)
    ax.set_extent(map_extent, crs=ccrs.PlateCarree())

    # Add the image *after* the extent of the map is set
    img, crs, extent, origin = dummy_image()
    ax.imshow(img, transform=crs, extent=extent, origin=origin, cmap="spring")

    # Add some map features
    ax.add_feature(cf.LAKES, alpha=0.3)
    ax.add_feature(cf.STATES)

    ## Add raster corners
    ax.scatter(
        extent[:2],
        extent[2:],
        transform=ccrs.PlateCarree(),
        s=100,
        c="k",
        marker="x",
        label="imshow extent",
    )

    # Final touches
    ax.gridlines(draw_labels=True)
    ax.set_title(f"\u2713 Using {type(axes_proj)}")
    ax.legend()
    fig.savefig("not_buggy_map.png")


def main():
    axes_proj = ccrs.LambertConformal(central_longitude=-87, central_latitude=41.7)
    buggy_map(axes_proj)
    not_buggy_map(axes_proj)


if __name__ == "__main__":
    main()

I have not found this being mentioned in the documentation though. Should a note be added?

[greglucas]

@edsaac, documentation improvements are always welcome! I think we have a few issues open relating to this same issue of ordering calls. Here is one saying also that documentation would help.
#1468