Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Quiver with all projections: bad direction of arrows #1179

Open
lguez opened this issue Nov 13, 2018 · 18 comments · May be fixed by #1926
Open

Quiver with all projections: bad direction of arrows #1179

lguez opened this issue Nov 13, 2018 · 18 comments · May be fixed by #1926

Comments

@lguez
Copy link

lguez commented Nov 13, 2018

Description

I have a wind field : eastward component u and northward component v, in m s-1. I want to make a quiver plot of this field in north stereographic projection. I would expect the correct command to be:

ax = plt.axes(projection =ccrs.NorthPolarStereo())
ax.quiver(lon, lat, u, v, transform = ccrs.PlateCarree(), angles = "xy")

( The complete test code is below.) However, this produces arrows in the wrong direction. It appears that the way to obtain the correct direction for the arrows is to write:

ax.quiver(lon, lat, u / np.cos(lat / 180 * np.pi) , v, transform = ccrs.PlateCarree(), angles = "xy")

Note that the magnitude of the vector is not correct then, we should also renormalize it to get the correct magnitude. Compare this with the behavior of basemap:

m = basemap.Basemap(projection="npstere", boundinglat = 50,  lon_0 = 0)
u_rot, v_rot, x, y = m.rotate_vector(u, v, lon ,lat, returnxy=True)
m.quiver(x, y, u_rot, v_rot, angles = "xy")

which produces correct arrows.

I think this is a bug, or a bad feature of cartopy. The problem comes from the transform_vectors method, as shown in the test code below. In this test code, we want to plot a single vector at longitude 0°, latitude 89°N. The vector is almost westward: u < 0 and v is very small compared to u. In the north stereographic projection, the components of the projected vector should be almost the same than the original components. We see with the test code below that the rotate_vector method of basemap correctly produces this result. However, the transform_vectors of cartopy surprisingly does not produce this result. To get the correct result with transform_vectors, we have to divide u by the cosine of latitude, and renormalize the result from transform_vectors. It seems that transform_vectors expects the derivative of longitude instead of the wind. This is quite awkward to use and counter-intuitive, I think, or a bug. The test code below goes on to plot the single vector, with basemap and with cartopy. We see the red arrow with cartopy in the wrong direction.

Code to reproduce

Here is the complete test code.

#!/usr/bin/env python3

from mpl_toolkits import basemap
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np

m = basemap.Basemap(projection="npstere", boundinglat = 50, lon_0 = 0)
proj = ccrs.NorthPolarStereo()
src_crs = ccrs.PlateCarree()

lon = np.array([0])
lat = np.array([89])
u = np.array([-3])
v = np.array([0.1])
print("m.rotate_vector(u, v, lon ,lat):", m.rotate_vector(u, v, lon ,lat))
print("proj.transform_vectors(src_crs, lon, lat, u, v):",
      proj.transform_vectors(src_crs, lon, lat, u, v))
u_rot, v_rot = proj.transform_vectors(src_crs, lon, lat,
                                      u / np.cos(lat /180 * np.pi), v)
print("u_rot, v_rot:", u_rot, v_rot)
renorm = np.sqrt((u**2 + v**2) / (u_rot**2 + v_rot**2))
print("array((u_rot, v_rot)) * renorm:", np.array((u_rot, v_rot)) * renorm)

plt.figure()
u_rot, v_rot, x, y = m.rotate_vector(u, v, lon ,lat, returnxy=True)
m.quiver(x, y, u_rot, v_rot, angles = "xy")
m.drawcoastlines()
plt.title("with basemap")

plt.figure()
ax = plt.axes(projection = proj)
ax.coastlines()
ax.set_extent((-180, 180, 50, 90), src_crs)
ax.quiver(lon, lat, u, v, transform = src_crs, angles = "xy", color = "red")
ax.quiver(lon, lat, u / np.cos(lat / 180 * np.pi), v, transform = src_crs,
          angles = "xy")
plt.title("with cartopy")

plt.show()

Here is the standard output of the script:

m.rotate_vector(u, v, lon ,lat): (array([-2.99999999]), array([ 0.10000035]))
proj.transform_vectors(src_crs, lon, lat, u, v): (array([-1.3922597]), array([ 2.65925044]))
u_rot, v_rot: [-171.80062661] [ 5.72817851]
array((u_rot, v_rot)) * renorm: [[-2.99999913]
 [ 0.10002601]]

And here are the two resulting figures:
figure_1
figure_2

Full environment definition

Operating system

Linux Mint 19 Cinnamon

Cartopy version

0.16.0

pip list

apt-clone (0.2.1)
apturl (0.5.2)
asn1crypto (0.24.0)
basemap (1.1.0)
beautifulsoup4 (4.6.0)
Brlapi (0.6.6)
bsddb3 (6.1.0)
Cartopy (0.16.0)
certifi (2018.1.18)
chardet (3.0.4)
command-not-found (0.3)
configobj (5.0.6)
cryptography (2.1.4)
cupshelpers (1.0)
cycler (0.10.0)
Cython (0.26.1)
decorator (4.1.2)
defer (1.0.6)
gpg (1.10.0)
gramps (4.2.8)
graphviz (0.5.2)
httplib2 (0.9.2)
idna (2.6)
ipython (5.5.0)
ipython-genutils (0.2.0)
louis (3.5.0)
lxml (4.2.1)
macaroonbakery (1.1.3)
Mako (1.0.7)
MarkupSafe (1.0)
matplotlib (2.1.1)
mutagen (1.38)
netCDF4 (1.3.1)
netifaces (0.10.4)
numpy (1.15.4)
onboard (1.4.1)
openshot-qt (2.4.1)
OWSLib (0.16.0)
PAM (0.4.2)
paramiko (2.0.0)
pexpect (4.2.1)
pickleshare (0.7.4)
Pillow (5.1.0)
pip (9.0.1)
prompt-toolkit (1.0.15)
protobuf (3.0.0)
psutil (5.4.2)
pyasn1 (0.4.2)
pycairo (1.16.2)
pycrypto (2.6.1)
pycups (1.9.73)
pycurl (7.43.0.1)
Pygments (2.2.0)
pygobject (3.26.1)
PyICU (1.9.8)
pyinotify (0.9.6)
pymacaroons (0.13.0)
PyNaCl (1.1.2)
pyparsing (2.2.0)
pyproj (1.9.5.1)
pyRFC3339 (1.0)
pyshp (2.0.1)
python-apt (1.6.3)
python-dateutil (2.6.1)
python-debian (0.1.32)
python-xapp (1.2.0)
python-xlib (0.20)
pytz (2018.3)
pyxdg (0.25)
PyYAML (3.12)
pyzmq (16.0.2)
reportlab (3.4.0)
requests (2.18.4)
requests-unixsocket (0.1.5)
scipy (0.19.1)
sessioninstaller (0.0.0)
setproctitle (1.1.10)
setuptools (40.6.0)
Shapely (1.6.4.post2)
simplegeneric (0.8.1)
six (1.11.0)
system-service (0.3)
systemd-python (234)
traitlets (4.3.2)
ubuntu-drivers-common (0.0.0)
ufw (0.35)
urllib3 (1.22)
wcwidth (0.1.7)
wheel (0.30.0)
xdot (0.9)
xkit (0.0.0)
youtube-dl (2018.11.7)
@lguez
Copy link
Author

lguez commented Nov 13, 2018

Made a test with the Mercator projection instead of north polar stereographic and the problem is the same. So I guess the problem is for all projections.

@lguez lguez changed the title Quiver with north stereographic projection: bad direction of arrows Quiver with all projections: bad direction of arrows Nov 13, 2018
@pelson
Copy link
Member

pelson commented Jan 15, 2019

Confirming that the transform_vectors algorithm is based on "derivatives of source coordinate system". This is documented at https://scitools.org.uk/cartopy/docs/latest/crs/index.html#cartopy.crs.CRS.transform_vectors, though there is no doubt that it could be made even more explicit (your help would be greatly appreciated there).

I believe there are three cases:

  • 1: We have x and y coords, and some u&v in a system invariant form
    • e.g. UTM x&y, with u and v in m/s in northward and eastward directions (which categorically are not aligned to the x and y axes of the projection)
    • we only want to transform x and y, and compute the new direction as a simple "new-Northward" rotation
  • 2: We have x and y coords, and u&v giving a direction defined in the same coordinate system as x&y, but not the magnitude
    • e.g. longs & lats, with u and v in m/s in the direction of the projection (not necessarily northward and eastward)
    • we want to use the existing algorithm - it converts the direction, and preserves the original magnitude
  • 3: We have x and y coords, and u&v giving a direction and magnitude defined in the same coordinate system as x&y
    • e.g. longs & lats, with u and v in degrees
    • we want to use the existing algorithm + compute the new length of the vector

cc @ajdawson as the original transform_vectors author.

@pelson
Copy link
Member

pelson commented Jan 15, 2019

In putting together the summary above, I had the following example:

src = ccrs.PlateCarree()
nx, ny = 6, 5
xs = np.linspace(*src.x_limits, nx+2)[1:-1]
ys = np.linspace(*src.y_limits, ny+2)[1:-1]

# u and v are in the x and y direction of src.
u = 5  # m/s, or some other non-src CRS bound system.
v = 0  # m/s

xs, ys = np.meshgrid(xs, ys)
us = np.tile(u, xs.shape)
vs = np.tile(v, ys.shape)

Which when visualised in the source CRS:

plt.figure(figsize=(10, 10))
ax = plt.axes(projection=src)
ax.coastlines()
ax.quiver(xs, ys, us, vs, angles="xy", color="red")

img1

And in a projected CRS:

plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0))
ax.coastlines()
ax.quiver(xs, ys, us, vs, angles="xy", color="red", transform=src)

img2

Reprojecting that case:

plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.NorthPolarStereo(central_longitude=0))
ax.coastlines()
ax.quiver(xs, ys, us, vs, angles="xy", color="red", transform=src, regrid_shape=20)

img3

It is quite clear that the magnitude of the vectors is preserved ref and that direction is based on the source grid.

@pelson
Copy link
Member

pelson commented Jan 16, 2019

I've thought over this some more, and now wonder if cases 1 and 2 are actually the same thing, and we are looking at a bug that is acute for this projection. In particular, I think the transform_vectors algorithm uses a forward difference to compute the rotation, but a central difference would provide a much better estimation.

@lguez
Copy link
Author

lguez commented Jan 21, 2019

Thank you for your answer. Maybe I am missing something in your line of thought but it seems to me that the example you have chosen does not illustrate the problem. The problem comes from the difference between the directions of (u /cos (latitude), v) and (u, v). So if you choose v = 0, there is no difference in direction. That is why I chose : u = -3, v =0.1 and latitude = 89°, in my post above.

Second, thanks for confirming that the transform_vectors algorithm is based on derivatives of source coordinate system. My point is then that I think the public interface of quiver with cartopy is very inconvenient. The majority of use cases of quiver in Cartopy is probably to plot a wind field in some projection, with the wind given by its eastward and northward components. In order to plot this wind field now in Cartopy, I have to do something like:

u_src_crs = u / np.cos(lat[:, np.newaxis] / 180 * np.pi)
v_src_crs = v
magnitude = ma.sqrt(u**2 + v**2)
magn_src_crs = ma.sqrt(u_src_crs**2 + v_src_crs**2)
ax.quiver(lon, lat, u_src_crs * magnitude / magn_src_crs, v_src_crs * magnitude / magn_src_crs,
                transform = src_crs)

I should not have to do all this just to plot a wind field. There should be a one-line command for this, included with cartopy.

@iljah
Copy link

iljah commented Mar 8, 2019

+1 for simple way of plotting vectors with starting point given by lat & lon and direction and length given by geographic north & east components in some physical units. In my case this would be geomagnetic and geoelectric fields with north and east components of nT and mV/km respectively.

@iljah
Copy link

iljah commented Mar 8, 2019

Ideally there would also be a boolean toggle to choose between making vector lengths depend only on the physical quantity (in addition automatic or user defined global scaling) or also depend on how far they seem to be located from the camera.

@akpetty
Copy link

akpetty commented Jul 20, 2020

Did anyone ever resolve this or find a way of doing something like rotate_vector without Basemap? That was a pretty useful function but it's hard to use this without needing to define the projection with Basemap too...

@arnea
Copy link

arnea commented Sep 23, 2020

Using the "angles" keyword in call to quiver seems to work fine.
Edit: using "angles" keyword, with array of angles in degrees seems to work fine.

@lguez
Copy link
Author

lguez commented Sep 23, 2020

I do not think it does in a convenient way : that is what I tried to show in submitting this issue. See the example in my post at the top of this issue.

@liuzheng-arctic
Copy link

I total agree with @lguez , the current user interface for vector rotation in Cartopy is inconvenient at the least. I don't think this issue is clearly documented so the first time users may not realize they need to rescale the u/v components (like @lguez's code above) before feeding the vectors into transform_vectors . This will end up with wrong vector maps. Can someone add a short example to showcase the correct way to rotate vectors to the documentation before there is a fix for this? I am sure I am not the only one trying to figure out why the wind maps look funny in the polar regions.

@iljah
Copy link

iljah commented Jan 29, 2021

+1 for documentation on (optionally simple) way of plotting vectors with starting point given by lat & lon and direction and length given by geographic north & east components in some physical units

@jabadge
Copy link

jabadge commented Oct 28, 2021

I definitely agree that this needs to be better documented at the very least. I would have continued to plot wind vectors in Antarctica incorrectly if I had not accidentally run across this post. In case anyone needs it, here's another example of the issue:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
from matplotlib import pyplot as plt

lat = np.array([[-87, -87, -87, -87, -87], [-82, -82, -82, -82, -82], [-77, -77, -77, -77, -77]])
lon = np.array([[0, 90, 180, -90, -180], [0, 90, 180, -90, -180], [0, 90, 180, -90, -180]])
u = np.array([[-3, -3, -3, -3, -3], [-3, -3, -3, -3, -3], [-3, -3, -3, -3, -3]]) # wind component pointing East
v = np.array([[0, .1, 1, 3, 1], [0, .1, 1, 3, 1], [0, .1, 1, 3, 1]]) # wind component pointing North

def plot_quiver(lon, lat, u, v):
    plt.figure()

    # Plot in South Polar Stereographic projection
    ax = plt.axes(projection = ccrs.SouthPolarStereo())
   
    # Set the extent of the plot
    ax.set_extent([-180, 180, -90, -65], crs = ccrs.PlateCarree())
    
    # Plot Antartic coastline for reference
    ax.add_feature(cfeature.COASTLINE)
    
    # Plot the arrows
    ax.quiver(lon, lat, u, v, transform = ccrs.PlateCarree(), angles = "xy")
    
    return

# Simply plot the data above without doing anything special
plot_quiver(lon, lat, u, v)

# Rotate and rescale the wind vectors using code from this GitHub issue
u_src_crs = u / np.cos(lat / 180 * np.pi)
v_src_crs = v
magnitude = np.sqrt(u**2 + v**2)
magn_src_crs = np.sqrt(u_src_crs**2 + v_src_crs**2)

plot_quiver(lon, lat, u_src_crs * magnitude / magn_src_crs, v_src_crs * magnitude / magn_src_crs)

image

image

The arrows at longitude 0 should point West.
The arrows at longitude 90 should point West and imperceptibly North.
The arrows at longitude 180 should point West and a little North.
The arrows at longitude 270 should point perfectly Northwest.

As you can see in the first figure, without the rotating and rescaling provided by Iguez, the arrows point in the wrong directions. With the rotation and rescaling, the arrows are correct (second figure).

@greglucas
Copy link
Contributor

I agree the current situation is less than ideal and there should be a way to handle this easier within Cartopy.

I opened #1920 with a proposed fix for this.

@greglucas greglucas linked a pull request Nov 7, 2021 that will close this issue
@leeyupeng119
Copy link

It has been 3 years and this problem still exists, I would like to know if there are plans to improve it in the future. Thank you in advance.

@greglucas
Copy link
Contributor

There is a linked PR with a proposed fix that no one has commented on: #1926

@axelschweiger
Copy link

Glad I found this after beating my ahead against this for half a day. It looks like the fix is still pending, in the meantime the cartopy quiver doc could probably use an update. Not sure what the mechanism for this is.

@generalpeng26
Copy link

I have the same problem. Have you resolved this problem?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet