Plotting a rasterio raster on a Cartopy GeoAxes

Question

I've seen a few other questions on this topic, but the library has changed enough that the answers to those no longer seem to apply.

Rasterio used to include an example for plotting a rasterio raster on a Cartopy GeoAxes. The example went roughly like this:

import matplotlib.pyplot as plt
import rasterio
from rasterio import plot

import cartopy
import cartopy.crs as ccrs

world = rasterio.open(r"../tests/data/world.rgb.tif")

fig = plt.figure(figsize=(20, 12))
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())
ax.set_global()
plot.show(world, origin='upper', transform=ccrs.PlateCarree(), interpolation=None, ax=ax)

ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS)

However, this code no longer draws the raster. Instead, I get something like this:

It should look like this:

When I asked about this in the rasterio issues tracker, they told me the example was deprecated (and deleted the example). Still, I wonder if there's some way to do what I'm trying to do. Can anyone point me in the right direction?

dgketchum · Answer 1 · 2021-12-29T21:19:04.167

I think you may want to read the data to a numpy.ndarray and plot it using ax.imshow, where ax is your cartopy.GeoAxes (as you have it already). I offer an example of what I mean, below.

I clipped a small chunk of Landsat surface temperature and some agricultural fields for this example. Get them on this drive link.

Note fields are in WGS 84 (epsg 4326), Landsat image is in UTM Zone 12 (epsg 32612), and I want my map in Lambert Conformal Conic. Cartopy makes this easy.

import numpy as np
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
import rasterio
import matplotlib.pyplot as plt


def cartopy_example(raster, shapefile):
    with rasterio.open(raster, 'r') as src:
        raster_crs = src.crs
        left, bottom, right, top = src.bounds
        landsat = src.read()[0, :, :]
        landsat = np.ma.masked_where(landsat <= 0,
                                     landsat,
                                     copy=True)
        landsat = (landsat - np.min(landsat)) / (np.max(landsat) - np.min(landsat))

    proj = ccrs.LambertConformal(central_latitude=40,
                                 central_longitude=-110)

    fig = plt.figure(figsize=(20, 16))
    ax = plt.axes(projection=proj)
    ax.set_extent([-110.8, -110.4, 45.3, 45.6], crs=ccrs.PlateCarree())

    shape_feature = ShapelyFeature(Reader(shapefile).geometries(),
                                   ccrs.PlateCarree(), edgecolor='blue')
    ax.add_feature(shape_feature, facecolor='none')
    ax.imshow(landsat, transform=ccrs.UTM(raster_crs['zone']),
              cmap='inferno',
              extent=(left, right, bottom, top))
    plt.savefig('surface_temp.png')


feature_source = 'fields.shp'
raster_source = 'surface_temperature_32612.tif'

cartopy_example(raster_source, feature_source)

The trick with Cartopy is to remember to use the projection keyword for your axes object, as this renders the map in a nice projection of your choice (LCC in my case). Use transform keyword to indicate what projection system your data is in, so Cartopy knows how to render it.

this is a great soln! Do you know if there is a way that the raster still shows up if the `facecolor='white'` or some other color for the cartopy vector features? I find that the cartopy features are drawn at the top — user308827, Feb 06 '22 at 18:33
This worked well but I had to change raster_crs['zone'] to just 12, the utm zone number for it to work. Has there been some change that causes raster_crs['zone'] to no longer work? Is there any other way we can get the zone number? — Casivio, Jun 16 '22 at 23:14

score -1 · Answer 2 · answered Jul 15 '19 at 16:22

-1

No need of rasterio. Get a bluemarble image, then plot it.

Here is the working code:

import cartopy
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

fig = plt.figure(figsize=(10, 5))
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())
# source of the image:
# https://eoimages.gsfc.nasa.gov/images/imagerecords/73000/73909/world.topo.bathy.200412.3x5400x2700.jpg
fname = "./world.topo.bathy.200412.3x5400x2700.jpg"
img_origin = 'lower'
img = plt.imread(fname)
img = img[::-1]
ax.imshow(img, origin=img_origin, transform=ccrs.PlateCarree(), extent=[-180, 180, -90, 90])

ax.coastlines()
ax.add_feature(cartopy.feature.BORDERS)

ax.set_global()
plt.show()

The output plot:

answered Jul 15 '19 at 16:22

swatchai

17,400
3
39
58

Yeah, but I need this to work on extremely high-resolution topography for the moon (as well as the earth example in the notebook), and your solution won't do that. I didn't just decide to use rasterio because I like cutting and pasting code. :) – Translunar Jul 26 '19 at 14:31
@DoctorMohawk You didn't mention the required resolution and the moon may be ++. I will wait and see the accepted answer. – swatchai Jul 26 '19 at 15:06
While you're correct that I didn't mention the required resolution, I *did* specifically say I wanted to render a rasterio raster on a Cartopy GeoAxes. – Translunar Jul 26 '19 at 18:36
How about "...I wonder if there's some way to do what I'm trying to do ..." – swatchai Jul 26 '19 at 18:41

Plotting a rasterio raster on a Cartopy GeoAxes

2 Answers2