edsaac
commented
6 days ago 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).
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
```py 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
Description
imshowworks 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 toimshow(Figure 3)To reproduce
Code
```py 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.