网格重新投影#

将格栅数据绘制到地图上时 imshow ,我们需要首先设置地图的范围,以便正确完成重新投影。

在本例中,我们将一些格栅数据存储为一个引用矩形坐标系(PlateCarree)的numpy数组。Cartopy根据当前设置的地图限制重新投影数据以匹配地图的坐标系。这意味着必须设置地图范围/边界 before imshow 被称为。

✗ Adding the raster BEFORE setting the map extent, ✓ Adding the raster AFTER setting the map extent
from matplotlib.patches import Rectangle
import matplotlib.pyplot as plt
import numpy as np

import cartopy.crs as ccrs
import cartopy.feature as cf


def main():
    # Generate raster data as a numpy array
    img = np.linspace(0, 1, 10_000).reshape(100, 100)

    # Define the origin and extent of the image following matplotlib's
    # convention `(left, right, bottom, top)`. These are referenced to
    # a rectangular coordinate system.
    img_origin = "lower"
    img_extent = (-87.6, -86.4, 41.4, 42.0)
    img_proj = ccrs.PlateCarree()

    imshow_kwargs = dict(
        extent=img_extent,
        origin=img_origin,
        transform=img_proj,
        cmap="spring",
    )

    # Define extent and projection for the map
    map_extent = (-88.1, -86.1, 41.2, 42.2)
    map_proj = ccrs.RotatedPole(pole_longitude=120.0, pole_latitude=70.0)

    fig, axs = plt.subplots(
        nrows=1,
        ncols=2,
        figsize=(12, 5),
        subplot_kw={"projection": map_proj},
        sharex=True,
        sharey=True,
        layout="constrained",
    )

    # Adding the raster *before* setting the map extent
    ax = axs[0]
    ax.set_title("\u2717 Adding the raster\nBEFORE setting the map extent")

    ax.imshow(img, **imshow_kwargs)
    ax.set_extent(map_extent, crs=img_proj)

    # Adding the raster *after* setting the map extent
    ax = axs[1]
    ax.set_title("\u2713 Adding the raster\nAFTER setting the map extent")

    ax.set_extent(map_extent, crs=img_proj)
    ax.imshow(img, **imshow_kwargs)

    for ax in axs:
        # Add other map features
        ax.add_feature(cf.LAKES, alpha=0.6)
        ax.add_feature(cf.STATES)

        # Highlight raster boundaries
        xy = (img_extent[0], img_extent[2])
        width = img_extent[1] - img_extent[0]
        height = img_extent[3] - img_extent[2]
        ax.add_patch(
            Rectangle(
                xy,
                width,
                height,
                transform=img_proj,
                edgecolor="black",
                facecolor="None",
                linewidth=3,
                label="Raster data bounds",
            )
        )

        ax.legend()
        ax.gridlines(draw_labels=True, x_inline=False, dms=True)

    plt.show()


if __name__ == "__main__":
    main()

Total running time of the script: (0分4.941秒)

Gallery generated by Sphinx-Gallery _