basemap icon indicating copy to clipboard operation
basemap copied to clipboard

CEA projection: plotting all the data works, plotting part of the data fails.

Open jakevdp opened this issue 10 years ago • 9 comments

Using the following setup, Python 2.7/matplotlib 1.2.1/basemap 1.0.6

from mpl_toolkits import basemap

m1 = basemap.Basemap(projection='cea', lon_0=0)
lon = [-135, -45, 45, 135]
lat = [45, 45, 45, 45]

If I then run the following command, the line plots as expected:

m1.plot(lon, lat, latlon=True)

But plotting just the first two points alone gives a failure:

m1.plot(lon[:2], lat[:2], latlon=True)

Here's the result:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-273-c6f75eb64ab8> in <module>()
      8 lat = [45, 45, 45, 45]
      9 #m1.plot(lon, lat, latlon=True)
---> 10 m1.plot(lon[:2], lat[:2], latlon=True)

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in with_transform(self, x, y, *args, **kwargs)
    525             if self.projection in _cylproj or self.projection in _pseudocyl:
    526                 if x.ndim == 1:
--> 527                     x = self.shiftdata(x)
    528                 elif x.ndim == 0:
    529                     if x > 180:

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in shiftdata(self, lonsin, datain, lon_0)
   4700             londiff = np.abs(lonsin[0:-1]-lonsin[1:])
   4701             londiff_sort = np.sort(londiff)
-> 4702             thresh = 360.-londiff_sort[-2]
   4703             itemindex = len(lonsin)-np.where(londiff>=thresh)[0]
   4704             if itemindex:

IndexError: index -2 is out of bounds for axis 0 with size 1

Edit: here's an error from a similar call, but I believe it's unrelated to the first:

lat = [-45, 45, 45, 45, 45, 0]
lon = [135, -135, -45, 45, 135, -90]
m1.plot(lon, lat, latlon=True)

error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-295-d815d338a0b6> in <module>()
     11 lon = [135, -135, -45, 45, 135, -90]
     12 
---> 13 m1.plot(lon, lat, latlon=True)
     14 

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in with_transform(self, x, y, *args, **kwargs)
    525             if self.projection in _cylproj or self.projection in _pseudocyl:
    526                 if x.ndim == 1:
--> 527                     x = self.shiftdata(x)
    528                 elif x.ndim == 0:
    529                     if x > 180:

/Users/jakevdp/anaconda/python.app/Contents/lib/python2.7/site-packages/mpl_toolkits/basemap/__init__.pyc in shiftdata(self, lonsin, datain, lon_0)
   4702             thresh = 360.-londiff_sort[-2]
   4703             itemindex = len(lonsin)-np.where(londiff>=thresh)[0]
-> 4704             if itemindex:
   4705                 # check to see if cyclic (wraparound) point included
   4706                 # if so, remove it.

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

jakevdp avatar Sep 25 '13 17:09 jakevdp

Here is a small workaround for you:

from mpl_toolkits import basemap
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt

m1 = basemap.Basemap(projection='cea', lon_0=0)
m1.drawcoastlines()
lon = [-135, -45, 45, 135]
lat = [45, 45, 45, 45]

x, y = m1(lon[:2], lat[:2])
ax = plt.gca()
ax.plot(*m1(lon[:2], lat[:2]))

This way (without latlon=True) the problematic function is not called... Probably latlon=True is there only for large arrays, and I suppose ideally even uniform..

Cheers

guziy avatar Aug 27 '15 15:08 guziy

Actually this also works fine:

from mpl_toolkits import basemap

m1 = basemap.Basemap(projection='cea', lon_0=0)
m1.drawcoastlines()
lon = [-135, -45, 45, 135]
lat = [45, 45, 45, 45]

m1.plot(*m1(lon[:2], lat[:2]))

guziy avatar Aug 27 '15 15:08 guziy

I think a good fix would be to check the latlon=True and the length of the arrays and throw an exception if the length is less than 3, with an explanation like: latlon=True is only for uniform grids with lots of points...

guziy avatar Aug 27 '15 15:08 guziy

Is latlon=True really only for a uniform grid? I've been using it pretty regularly for any input – I assumed that it basically means "call m(x, y) automatically", and it seems to act that way for arrays longer than 2 elements.

jakevdp avatar Aug 27 '15 17:08 jakevdp

Actually it is also calling shiftdata and the shiftdata is kind of limited, here is the doc to shiftdata:

def shiftdata(self,lonsin,datain=None,lon_0=None):
    """ Shift longitudes
    (and optionally data) so that they match map projection region. Only valid
    for cylindrical/pseudo-cylindrical global projections and data on regular
    lat/lon grids. longitudes and data can be 1-d or 2-d, if 2-d it is assumed
    longitudes are 2nd (rightmost) dimension.
    """

It actually might work for nonuniform grids, but I am not sure if it is bulletproof...

guziy avatar Aug 27 '15 19:08 guziy

Duplicating the code here, I could not make it highlighted above, probably because I sent it using email...

    def shiftdata(self,lonsin,datain=None,lon_0=None):
        """ Shift longitudes
        (and optionally data) so that they match map projection region. Only valid
        for cylindrical/pseudo-cylindrical global projections and data on regular
        lat/lon grids. longitudes and data can be 1-d or 2-d, if 2-d it is assumed
        longitudes are 2nd (rightmost) dimension.
        """

guziy avatar Aug 27 '15 19:08 guziy

I am pretty sure that shiftdata is not needed for scatter, so I've created a PR for that, but I am not sure about plot though...

guziy avatar Aug 27 '15 19:08 guziy

You need it for plot so that a path drawn around over the pacific doesn't look like it crosses the Atlantic, for example.

On Thu, Aug 27, 2015 at 3:20 PM, Huziy Oleksandr (Sasha) < [email protected]> wrote:

I am pretty sure that shiftdata is not needed for scatter, so I've created a PR for that, but I am not too sure about plot though...

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/126#issuecomment-135528647.

WeatherGod avatar Aug 27 '15 19:08 WeatherGod

We have also encountered issues using latlon=True and I'm wondering if it wouldn't hurt to tighten up the documentation to note some of these complexities with its use. It appears a better practice to not use it if possible. However, for someone just starting to use basemap, it is very appealing and further instructions on its usage would be helpful to avoid the issues we were having in https://github.com/MPAS-Dev/geometric_features/issues/14.

pwolfram avatar Feb 08 '16 14:02 pwolfram