cartopy - 在 Cartopy GeoAxes 上绘制光栅栅格

Question

我已经看到了一些关于这个主题的其他问题，但是库已经发生了很大的变化，以至于这些问题的答案似乎不再适用。

栅格 used to include an example用于在 Cartopy GeoAxes 上绘制光栅栅格。这个例子大致是这样的:

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)

但是，这段代码不再绘制光栅。相反，我得到这样的东西:

它应该是这样的:

当我在 rasterio 问题跟踪器中询问这个问题时，他们告诉我这个例子已被弃用(并删除了这个例子)。不过，我想知道是否有某种方法可以做我想做的事情。谁能指出我正确的方向？

Answer 1

我认为您可能希望将数据读取到 numpy.ndarray 并使用 ax.imshow 绘制它，其中 ax 是您的cartopy.GeoAxes(你已经有了)。我在下面提供了一个例子来说明我的意思。

为此示例，我剪裁了一小块 Landsat 表面温度和一些农田。让他们在这drive link .

注意字段在 WGS 84 (epsg 4326) 中，Landsat 图像在 UTM Zone 12 (epsg 32612) 中，我希望我的 map 在 Lambert Conformal Conic 中。 Cartopy 让这一切变得简单。

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)

Cartopy 的诀窍是记住为轴对象使用 projection 关键字，因为这会以您选择的漂亮投影(在我的示例中为 LCC)渲染 map 。使用 transform 关键字来指示您的数据所在的投影系统，以便 Cartopy 知道如何渲染它。