在 Cartopy GeoAxes 上绘制栅格栅格

Question

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

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

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)，我希望我的地图采用 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）。使用 transform 关键字指明您的数据所在的投影系统，以便 Cartopy 知道如何渲染它。

Answer 2

不需要rasterio。获取 bluemarble 图像，然后绘制它。

这是工作代码：

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()

输出图：