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

Question

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

Rasterio曾经包含一个在 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 问题跟踪器中询问这个问题时，他们告诉我该示例已被弃用（并删除了该示例）。不过，我想知道是否有某种方法可以做我想做的事情。谁能指出我正确的方向？

score 1 · Accepted Answer

我认为您可能希望将数据读取到 anumpy.ndarray并使用绘制它，您的ax.imshow位置在哪里（因为您已经拥有它）。我在下面提供了一个例子来说明我的意思。axcartopy.GeoAxes

我为这个例子剪下了一小块 Landsat 表面温度和一些农田。在此驱动器链接上获取它们。

注意字段在 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 知道如何渲染它。

score -1 · Accepted Answer

不需要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()

输出图：

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

2 回答 2

Related

Reference