Plotting in Different Projections#

This example shows how to overlay data and graphics in different projections, demonstrating various features of Iris, Cartopy and matplotlib.

We wish to overlay two datasets, defined on different rotated-pole grids. To display both together, we make a pseudocoloured plot of the first, overlaid with contour lines from the second. We also add some lines and text annotations drawn in various projections.

We plot these over a specified region, in two different map projections.

  • A pseudocolour plot on the Equidistant Cylindrical projection, with overlaid contours.
  • A pseudocolour plot on the North Polar Stereographic projection, with overlaid contours.
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

import iris
import iris.plot as iplt

# Define a Cartopy 'ordinary' lat-lon coordinate reference system.
crs_latlon = ccrs.PlateCarree()


def make_plot(projection_name, projection_crs):
    # Create a matplotlib Figure.
    plt.figure()

    # Add a matplotlib Axes, specifying the required display projection.
    # NOTE: specifying 'projection' (a "cartopy.crs.Projection") makes the
    # resulting Axes a "cartopy.mpl.geoaxes.GeoAxes", which supports plotting
    # in different coordinate systems.
    ax = plt.axes(projection=projection_crs)

    # Set display limits to include a set region of latitude * longitude.
    # (Note: Cartopy-specific).
    ax.set_extent((-80.0, 20.0, 10.0, 80.0), crs=crs_latlon)

    # Add coastlines and meridians/parallels (Cartopy-specific).
    ax.coastlines(linewidth=0.75, color="navy")
    ax.gridlines(crs=crs_latlon, linestyle="-")

    # Plot the first dataset as a pseudocolour filled plot.
    maindata_filepath = iris.sample_data_path("rotated_pole.nc")
    main_data = iris.load_cube(maindata_filepath)
    # NOTE: iplt.pcolormesh calls "pyplot.pcolormesh", passing in a coordinate
    # system with the 'transform' keyword:  This enables the Axes (a cartopy
    # GeoAxes) to reproject the plot into the display projection.
    iplt.pcolormesh(main_data, cmap="RdBu_r")

    # Overplot the other dataset (which has a different grid), as contours.
    overlay_filepath = iris.sample_data_path("space_weather.nc")
    overlay_data = iris.load_cube(overlay_filepath, "total electron content")
    # NOTE: as above, "iris.plot.contour" calls "pyplot.contour" with a
    # 'transform' keyword, enabling Cartopy reprojection.
    iplt.contour(overlay_data, 20, linewidths=2.0, colors="darkgreen", linestyles="-")

    # Draw a high resolution margin line, inset from the pcolormesh border.
    # First calculate rectangle corners, 7% in from each corner of the data.
    x_coord, y_coord = main_data.coord(axis="x"), main_data.coord(axis="y")
    x_start, x_end = np.min(x_coord.points), np.max(x_coord.points)
    y_start, y_end = np.min(y_coord.points), np.max(y_coord.points)
    margin = 0.07
    margin_fractions = np.array([margin, 1.0 - margin])
    x_lower, x_upper = x_start + (x_end - x_start) * margin_fractions
    y_lower, y_upper = y_start + (y_end - y_start) * margin_fractions
    steps = np.linspace(0, 1)
    zeros, ones = np.zeros(steps.size), np.ones(steps.size)
    x_delta, y_delta = (x_upper - x_lower), (y_upper - y_lower)
    x_points = x_lower + x_delta * np.concatenate((steps, ones, steps[::-1], zeros))
    y_points = y_lower + y_delta * np.concatenate((zeros, steps, ones, steps[::-1]))
    # Get the Iris coordinate system of the X coordinate (Y should be the same).
    cs_data1 = x_coord.coord_system
    # Construct an equivalent Cartopy coordinate reference system ("crs").
    crs_data1 = cs_data1.as_cartopy_crs()
    # Draw the rectangle in this crs, with matplotlib "pyplot.plot".
    # NOTE: the 'transform' keyword specifies a non-display coordinate system
    # for the plot points (as used by the "iris.plot" functions).
    plt.plot(
        x_points,
        y_points,
        transform=crs_data1,
        linewidth=2.0,
        color="white",
        linestyle="--",
    )

    # Mark some particular places with a small circle and a name label...
    # Define some test points with latitude and longitude coordinates.
    city_data = [
        ("London", 51.5072, 0.1275),
        ("Halifax, NS", 44.67, -63.61),
        ("Reykjavik", 64.1333, -21.9333),
    ]
    # Place a single marker point and a text annotation at each place.
    for name, lat, lon in city_data:
        plt.plot(
            lon,
            lat,
            marker="o",
            markersize=7.0,
            markeredgewidth=2.5,
            markerfacecolor="black",
            markeredgecolor="white",
            transform=crs_latlon,
        )
        # NOTE: the "plt.annotate call" does not have a "transform=" keyword,
        # so for this one we transform the coordinates with a Cartopy call.
        at_x, at_y = ax.projection.transform_point(lon, lat, src_crs=crs_latlon)
        plt.annotate(
            name,
            xy=(at_x, at_y),
            xytext=(30, 20),
            textcoords="offset points",
            color="black",
            backgroundcolor="white",
            size="large",
            arrowprops=dict(arrowstyle="->", color="white", linewidth=2.5),
        )

    # Add a title, and display.
    plt.title(
        "A pseudocolour plot on the {} projection,\nwith overlaid contours.".format(
            projection_name
        )
    )
    iplt.show()


def main():
    # Demonstrate with two different display projections.
    make_plot("Equidistant Cylindrical", ccrs.PlateCarree())
    make_plot("North Polar Stereographic", ccrs.NorthPolarStereo())


if __name__ == "__main__":
    main()

Total running time of the script: (0 minutes 0.609 seconds)

Gallery generated by Sphinx-Gallery