quiver arrows reverse in all polar projections except orthographic. Why?

[joakimkjellsson]

Hi

I'm trying to plot vectors in the Southern ocean using basemap and matplotlib.
To illustrate my problem I've created a zonal jet at 60S.
I have U and V as 2D arrays and lon and lat as 1D arrays, all with cyclic point added and lon is from -180 to 180.
If I plot in orthographic projection the vectors correctly shows eastward motion, but if I use any of the polar projections (spstere, splaea spaeqd) the vectors reverse!
This is really strange and oddly enough it does not seem to happen if I make the same plot over the North Pole.

I've noted that in the southern hemisphere lon_0 has a different meaning for orthographic projection than for spstere, splaea etc.
So to orient the maps the same way I must pass lon_0=0 in orthographic projection and lon_0=180 in polar stereographic projection if I plot southern hemisphere.
But this is not the case in the Northern hemisphere.
Could this be part of the problem?

Here's the script that creates the results:
https://dl.dropboxusercontent.com/u/86754720/plot_quiver_bug_test.py
Here's the contourf of u and v
https://dl.dropboxusercontent.com/u/86754720/velocity_bug_contourf.png
and here are the vectors
https://dl.dropboxusercontent.com/u/86754720/velocity_bug_quiver.png

[jswhit]

For reasons not yet clear to me, you have to reverse the sign of the vectors that come back from transform_vector for SH polar plots, i.e.

ups = -ups; vps = -vps

in your test script right after the call to transform_vector. Must have something to do with how the projection coordinates are defined in proj4.

[joakimkjellsson]

Thanks. That solved it. You're probably right in that its not a basemap issue per se but rather a proj4 thing...

[efiring]

Still, it seems to me that transform_vector should take care of this transparently, so I am reopening this.

[joanpau]

I think that the reason for this behavior is that when using polar projections at the southern hemisphere the ranges of the axes are negative:

In [30]: m = Basemap(projection='spstere', boundinglat=0, lon_0=0)

In [31]: plt.gca().get_xlim()
Out[31]: (0.0, -25483987.999999996)

In [32]: plt.gca().get_ylim()
Out[32]: (0.0, -25483987.999999996)

combined with the default value of the option angles = 'uv' of Basemap's quiver (inherited from Matplotlib's quiver):

angles: [ ‘uv’ | ‘xy’ | array ]
    With the default ‘uv’, the arrow axis aspect ratio is 1, so that if U*==*V
    the orientation of the arrow on the plot is 45 degrees CCW from the horizontal axis
    (positive to the right). With ‘xy’, the arrow points from (x,y) to (x+u, y+v).
    Use this for plotting a gradient field, for example. Alternatively, arbitrary angles
    may be specified as an array of values in degrees, CCW from the horizontal axis.
    Note: inverting a data axis will correspondingly invert the arrows only with angles='xy'.

When using 'angles=xy', vectors are plotted properly.
This is consistent with the computations in rotate_vector.
However, there is something that I can not understand in that function,
namely code in lines 3137-3142:

        # Define a displacement (dlon, dlat) that moves all
        # positions (lons, lats) a small distance in the
        # direction of the original vector.
        dc = 1E-5 * np.exp(theta*1j)
        dlat = dc.imag * np.cos(np.radians(lats))
        dlon = dc.real

What is the foundation of that computation? Should not this be:

        # Define a displacement (dlon, dlat) that moves all
        # positions (lons, lats) a small distance in the
        # direction of the original vector.
        dc = 1E-5 * np.exp(theta*1j)
        dlat = dc.imag
        dlon = dc.real / np.cos(np.radians(lats))

???

[guziy]

Hi:

Yes you are right it should be divided by cos in dlon, but that might cause a division by 0 problem, so I suppose he multiplied the vector (i.e. both components) by cos. And it is OK to multiply since cos(lat) is not negative and preserves direction.

Cheers

[joanpau]

Thanks @guzly. I came to that conclusion too. I just wondered
if there were a different geometric reasoning running here I was not aware of.
Also, from this it can be deduced that the input components u and v
should be in metric units, and not in angular units.
I mean that for the case of velocities they should be in some multiple of m/s
and not deg/s. I was not able to deduce this from the documentation:

The input vector field is defined in spherical coordinates (it
has eastward and northward components) while the output
vector field is rotated to map projection coordinates (relative
to x and y). The magnitude of the vector is preserved.

The wording does not specify it, and other mapping systems use angular units
instead of metric units for vector component inputs
(e.g. quiverm from MATLAB Mapping Toolbox)
I think that this should be more clearly stated in the docstring.
Do you think this deserves a new issue and/or pull request?

[guziy]

I think this function would work for any units. Spherical coordinates, as I
understand, here means that velocity components are along meridians and
parallels... But let's see what others think, maybe the documentation would
benefit from additional clarification...

Cheers

2015-09-14 7:28 GMT-04:00 Joan Pau Beltran notifications@github.com:

Thanks @guzly. I came to that conclusion too. I just wondered
if there were a different geometric reasoning running here I was not aware
of.
Also, from this it can be deduced that the input components u and v
should be in metric units, and not in angular units.
I mean that for the case of velocities they should be in some multiple of
m/s
and not deg/s. I was not able to deduce this from the documentation:

The input vector field is defined in spherical coordinates (it
has eastward and northward components) while the output
vector field is rotated to map projection coordinates (relative
to x and y). The magnitude of the vector is preserved.

The wording does not specify it, and other mapping systems use angular
units
instead of metric units for vector component inputs
(e.g. quiverm from MATLAB Mapping Toolbox
http://www.mathworks.com/help/map/ref/quiverm.html)
I think that this should be more clearly stated in the docstring.
Do you think this deserves a new issue and/or pull request?

—
Reply to this email directly or view it on GitHub
#157 (comment).

Sasha

[joanpau]

Not actually. If the velocity components are eastward and northward but already expressed in deg/s
the factor np.cos(np.radians(lats)) in the above code should not be there, leaving instead:

        dc = 1E-5 * np.exp(theta*1j)
        dlat = dc.imag
        dlon = dc.real

Since only one component is affected, this is not a scaling: the direction of the
approximated differential increment is different.

I am not asking to change the current behavior to this one.
I am just saying that the expected unit convention of the input vector components (metric vs angular)
could be more explicitly stated in the docstring.

[guziy]

Maybe it is a bit tricky to understand... I'll try to explain later or
maybe someone who knows better could help.... I have some deadlines for
tomorrow...

Cheers

2015-09-14 9:00 GMT-04:00 Joan Pau Beltran notifications@github.com:

Not actually. If the velocity components are eastward and northward but
already expressed in deg/s:
the factor np.cos(np.radians(lats)) in the above code should not be
there, leaving instead:

    dc = 1E-5 * np.exp(theta*1j)
    dlat = dc.imag
    dlon = dc.real

Since only one component is affected, this is not a scaling: the direction
of the
approximated differential increment is different.

—
Reply to this email directly or view it on GitHub
#157 (comment).

Sasha

[guziy]

Actually ...

You are right, sorry)) I hope nobody will try to use it this way....
Because in this case units for the components are different though they
have the same name and in order to calculate the module one of them should
be divided or multiplied by cos ....

Confusing ...))

2015-09-14 9:26 GMT-04:00 Oleksandr Huziy guziy.sasha@gmail.com:

Maybe it is a bit tricky to understand... I'll try to explain later or
maybe someone who knows better could help.... I have some deadlines for
tomorrow...

Cheers

2015-09-14 9:00 GMT-04:00 Joan Pau Beltran notifications@github.com:

Not actually. If the velocity components are eastward and northward but
already expressed in deg/s:
the factor np.cos(np.radians(lats)) in the above code should not be
there, leaving instead:

    dc = 1E-5 * np.exp(theta*1j)
    dlat = dc.imag
    dlon = dc.real

Since only one component is affected, this is not a scaling: the
direction of the
approximated differential increment is different.

—
Reply to this email directly or view it on GitHub
http://t.sidekickopen18.com/e1t/c/5/f18dQhb0S7lC8dDMPbW2n0x6l2B9nMJN7t5XZsf6M1jN2B8hqKdnGrCVd_V9H56dHzydWxrcl02?t=https%3A%2F%2Fgithub.com%2Fmatplotlib%2Fbasemap%2Fissues%2F157%23issuecomment-140063177&si=5057221580292096&pi=9bb774f1-28b6-4766-d4a4-2dcf8506bc70
.

Sasha

Sasha