wrapping of longitude

[petrushy]

I am trying to scatter plot on maps with nonzero middle longitude, but the areas where the map is "wrapped" are not plotted.

My test script:

---

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

m = Basemap(lon_0=160)
m.fillcontinents()

for x in range(-179, 179):
px1, py1 = m(x, x/4.0)
m.scatter(px1,py1, zorder=3)

plt.show()

---

If lon_0 is zero it works fine but with non-zero lat_0 some areas are not plotted.

It seems like the m(lon,lat) is not wrapping the points, the points I expect to wrap are negative in px1?

Or am I missing some transformation step?

[jswhit]

You need to make sure the longitudes you are trying to plot are within the map region (lon_0-180 to lon_0+180). In your example, this means adding 180 to each value of x. Basemap does not automatically do this for you, although I suppose it could. I've shied away from building this in, because I'm not sure that it is what users expect in all cases.

[petrushy]

Ok thanks, I'll just wrap it manually to that range.

I am not an expert at all of how the library is used, but my spontanous feeling is that the mapping from a lat,lon to a px,py would correspond to the same pixel location as the underlying map for that lat lon.

Thanks for the clarification!

[jswhit]

I see your point. Right now, for better or worse, Basemap.call is a noop for the 'cyl' projection (it just returns the same longitudes you give it). I'll investigate the possibility of having it transform the longitudes to the plot region.

[WeatherGod]

As an additional data point, I have a co-worker who recently ran into similar issues with pcolormesh(). The coordinate data was monotonic (0, 360), but the display limits was (-180, 180), so half of the pcolor image looked smeared. Because it is basemap that is performing an implicit setting of the limits, I would agree that basemap should also be smart about handling user inputs. At the extreme least, it should be documented in each of the functions how to deal with it (introducing users to numpy's unwrap, as well as possible needs to resort the data). Something else that might be useful is if the problem could be detected somehow and a warning issued.

But, ultimately, I think the best goal would be to either have basemap either provide the tools for rectifying input data, or to just simply do it for the user.

[jswhit]

Here's a prototype function to do what's needed:

import numpy as np

def shiftlon(lon,lon_0):
    """returns original sequence of longitudes (in degrees) recentered
    in the interval [lon_0-180,lon_0+180]"""
    lon_shift = np.asarray(lon)
    lon_shift = np.where(lon_shift > lon_0+180, lon_shift-360 ,lon_shift)
    lon_shift = np.where(lon_shift < lon_0-180, lon_shift+360 ,lon_shift)
    itemindex = len(lon)-np.where(lon_shift[0:-1]-lon_shift[1:] >= 180)[0]
    return np.roll(lon_shift,itemindex-1)

if __name__ == "__main__":
    lon_0 = 90
    lon = -180+5*np.arange(72)
    lon_shift = shiftlon(lon,lon_0)
    print "original lon",lon
    print " shifted lon",lon_shift

I'll modify my fork so that when longitudes are transformed to map projection coordinates for 'cyl','merc' or 'mill' this transformation is applied.

[jswhit]

I just added this feature to my fork (https://github.com/jswhit/basemap). The OP's test script now works.

[jswhit]

This is more subtle than I though at first. The longitudes are now only transformed to [lon_0-180,lon_0+180] before the forward transform. numpy.roll is not applied, since it only makes sense for ordered, 1-d sequences. This will not address the problem mentioned by Ben (trying to contour a field defined in 0 to 360 on a map region defined on -180 to 180). To accomplish this automatically requires shifting not only the input longitudes, but the data to be plotted. This can be done manually using the shiftgrid function, but I hesitate to make this automatic since it just seems to be a bit too much under-the-hood magic.

[jswhit]

OK, I backed off the change that automatically shifts the longitudes to the map region in Basemap.call. The convenience function shiftlons should be used to do this manually. There are just too many cases where automatically applying the transform produces surprising results. I also created a Basemap.shiftdata method that addresses the problem Ben mentioned. It's use is demonstrated by examples/shiftdata.py.

[jswhit]

Since both the data and the longitudes need to be shifted together, it makes more sense to do this in the plotting methods themselves, instead of in the transformation method (call). I've added a latlon keyword to pcolor,pcolormesh,contour,contourf,quiver and barb. The default is False, but it is set to True, the x,y arguments are interpreted as longitudes and latitudes instead of map projection coordinates. The data and longitudes are then automatically shifted to match the map projection region (using the shiftdata method) and the lons/lats are then converted to map projection coordinates. The shiftdata.py example exercises this functionality.

[WeatherGod]

I agree that making the object call function perform some sort of magic isn't the way to go. That would impact users who may be using their basemap object as an quick-n-dirty interface to pyproj4 and may not be interested in creating plots. The approach of having it in the plotting functions makes much more sense.

Just thinking out loud here, it would be nice if mpl's plotting function had a decorator system that would enable 3rd party tools to add special actions before and after a plotting action happens. This would reduce the amount of code basemap would need to maintain as many of basemap's plotting functions have repeated code, and users could potentially write their basemap using code in a manner very similar to non-basemap using code. Furthermore, updates to mpl's plotting functions would be easier to track. Of course, not all of mpl's plotting functions make sense to basemap, such as pie().

[jswhit]

Ben - can you provide a example of how this might work and what it might be used for?

[petrushy]

Just wanted to say thanks for investigating this issue, to me right now it seems clean to have the latlong option in the plot function themselves.

[jswhit]

issued a pull request for this (#56). Following Ben's suggestion, the code to shift the data and longitudes was refactored into a decorator.