Search code examples
pythoncartopy

Wrong arrow length using quiver with cartopy projections

I want to plot a vector field with vectors representing a displacement between one point to another on the map with cartopy.

My code works as expected when using the PlateCarree() transformation, but arrow length is several orders of magnitude off for all the other projections I tested.

Here is a MWE that should illustrate quite clearly the issue:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

import numpy

# Want to test for several different projections
projections = [
    ccrs.PlateCarree(),
    ccrs.EqualEarth(),
    ccrs.Mollweide(),
    ccrs.AzimuthalEquidistant(),
]
# ALl the coordinates will be given in the PlateCarree coordinate system.
coordinate_ccrs = ccrs.PlateCarree()

# We want N**2 points over the latitude/longitude values.
N = 5
lat, lon = numpy.meshgrid(numpy.linspace(-80, 80, N), numpy.linspace(-170, 170, N))
lat, lon = lat.flatten(), lon.flatten()

# We want arrows to appear, let make a small perturbation and try
# to do an arrow from (lon, lat) to (lon + perturbation, lat + perturbation).
rng = numpy.random.default_rng()
perturbation_amplitude = 10
lat_perturbation = perturbation_amplitude * rng.random(N * N)
lon_perturbation = perturbation_amplitude * rng.random(N * N)

# Create the matplotlib figure and axes, no projection for the moment as this
# will be changed later.
fig, axes = plt.subplots(2, 2)
axes = axes.flatten()

for i, projection in enumerate(projections):
    # Replace the existing ax with an ax with the desired projection.
    ax = axes[i]
    fig.delaxes(ax)
    ax = axes[i] = fig.add_subplot(2, 2, i + 1, projection=projection)
    # Make the plot readable.
    ax.set_global()
    ax.gridlines(draw_labels="x")

    # Non pertubed points are plotted in black.
    ax.plot(lon, lat, "k.", ms=5, transform=coordinate_ccrs)
    # Perturbed points are plotted in red.
    ax.plot(
        lon + lon_perturbation,
        lat + lat_perturbation,
        "r.",
        ms=5,
        transform=coordinate_ccrs,
    )
    # We try to draw arrows from a given black dot to its corresponding
    # red dot.
    ax.quiver(
        lon,
        lat,
        lon_perturbation,
        lat_perturbation,
        transform=coordinate_ccrs,
        # From https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.quiver.html?highlight=quiver#matplotlib.axes.Axes.quiver
        # look at the documentation of the "scale_unit" parameter.
        # The next 3 parameters are what matplotlib tell us to do. From
        # matplotlib documentation:
        #   To plot vectors in the x-y plane, with u and v having the same units
        #   as x and y, use angles='xy', scale_units='xy', scale=1.
        angles="xy",
        scale_units="xy",
        scale=1,
        # Simply make the arrows nicer, removing these last 3 parameters do not
        # solve the issue.
        minshaft=2,
        minlength=0.5,
        width=0.002,
    )

# Show everything
plt.show()

which display on the screen the following image:

enter image description here

The PlateCarree transformation is the only one displaying arrows. In fact, there are arrows in the other 3 projections, but I order to see them I need to scale them by 10000 with scale=0.00001 in the call to quiver, which gives:

enter image description here

Did I make a mistake when using cartopy API, is this expected behaviour and I missed something in the documentation, or is this a bug?

Solution

  • while there's quite some debate on github about cartopy's implementation of quiver-plot transformations GitHub-issues there is in fact a way on how to get your plot look as you want it to look...

    However, while thinking about this... I noticed that there's a thing that you might want to consider when using projected quiver-plots...

    As I see it, the re-projected arrows would would most probably need to be curved to really visualize the same direction as provided in the original data!

    (in the input-crs the arrow points as a straight line from point A to B, but if you re-project the points, the "straight line" that connected A and B is now in general a curved line, and so if the original direction was correct, I think the new direction should be indicated as a curved arrow...)

    That being said, you could achieve what you want by transforming the points manually instead of letting cartopy do the job:

    import matplotlib.pyplot as plt
import cartopy.crs as ccrs

import numpy

# Want to test for several different projections
projections = [
    ccrs.PlateCarree(),
    ccrs.EqualEarth(),
    ccrs.Mollweide(),
    ccrs.AzimuthalEquidistant(),
]
# ALl the coordinates will be given in the PlateCarree coordinate system.
coordinate_ccrs = ccrs.PlateCarree()

# We want N**2 points over the latitude/longitude values.
N = 5
lat, lon = numpy.meshgrid(numpy.linspace(-80, 80, N), numpy.linspace(-170, 170, N))
lat, lon = lat.flatten(), lon.flatten()

# We want arrows to appear, let make a small perturbation and try
# to do an arrow from (lon, lat) to (lon + perturbation, lat + perturbation).
rng = numpy.random.default_rng()
perturbation_amplitude = 10
lat_perturbation = perturbation_amplitude * rng.random(N * N)
lon_perturbation = perturbation_amplitude * rng.random(N * N)

# Create the matplotlib figure and axes, no projection for the moment as this
# will be changed later.
fig, axes = plt.subplots(2, 2)
axes = axes.flatten()

for i, projection in enumerate(projections):
    # Replace the existing ax with an ax with the desired projection.
    ax = axes[i]
    fig.delaxes(ax)
    ax = axes[i] = fig.add_subplot(2, 2, i + 1, projection=projection)
    # Make the plot readable.
    ax.set_global()
    ax.gridlines(draw_labels="x")

    # Non pertubed points are plotted in black.
    ax.plot(lon, lat, "k.", ms=5, transform=coordinate_ccrs)
    # Perturbed points are plotted in red.
    ax.plot(
        lon + lon_perturbation,
        lat + lat_perturbation,
        "r.",
        ms=5,
        transform=coordinate_ccrs,
    )
    
    
    xy_start = projection.transform_points(coordinate_ccrs, lon, lat)[:,:-1].T
    xy_end = projection.transform_points(coordinate_ccrs, lon + lon_perturbation, 
                                         lat + lat_perturbation)[:,:-1].T
    
    # We try to draw arrows from a given black dot to its corresponding
    # red dot.
    ax.quiver(
        *xy_start, 
        *(xy_end - xy_start),
        # From https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.quiver.html?highlight=quiver#matplotlib.axes.Axes.quiver
        # look at the documentation of the "scale_unit" parameter.
        # The next 3 parameters are what matplotlib tell us to do. From
        # matplotlib documentation:
        #   To plot vectors in the x-y plane, with u and v having the same units
        #   as x and y, use angles='xy', scale_units='xy', scale=1.
        angles="xy",
        scale_units="xy",
        scale=1,
        # Simply make the arrows nicer, removing these last 3 parameters do not
        # solve the issue.
        minshaft=2,
        minlength=0.5,
        width=0.002,
    )

# Show everything
plt.show()

    enter image description here