Strange behaviour on plotting data on Ronbinson projection using Basemap

Dear all,

I am trying to plot some data (see attached data.txt) on global coverage with 1-degree resolution on the Robinson projection using Basemap. However I get some strange band on the high latitude, and imshow function by matplotlib does not show similar thing. Please refer to the two attached figures.

Could anyone give me some tips? thanks!!! Below is a working example:

import numy as np

import mpl_toolkits.basemap as bmp

import matplotlib.pyplot as plt

#The example file data.txt could be downloaded from dropbox:
https://www.dropbox.com/s/xma4w540qa83sa6/data.txt

RonbinsonB.jpg

imgB.jpg

···

#############
data=np.genfromtxt(‘data.txt’,usemask=True,missing_values=‘0.000000000000000000e+00’)

print np.ma.unique(data)
lev = np.arange(0.5,8.6,1)
print lev

here we have to first build the equal-distance grid, this is inspired from here:

http://matplotlib.1069221.n5.nabble.com/Basemap-plotting-data-on-projection-td40973.html

cyl_basemap = bmp.Basemap(projection=‘cyl’, llcrnrlat=-90, urcrnrlat=90,
llcrnrlon=-180, urcrnrlon=180, resolution=‘l’)
lon, lat = cyl_basemap.makegrid(360, 180)

fig,ax=plt.subplots(1,1)
m=bmp.Basemap(projection=‘robin’,lon_0=0,resolution=‘c’,ax=ax)
m.drawcoastlines()
x, y = m(lon, np.flipud(lat))
m.contourf(x,y,data=data,levels=lev)

#############################

Cheers,

Chao


please visit:
http://www.globalcarbonatlas.org/


Chao YUE
Laboratoire des Sciences du Climat et de l’Environnement (LSCE-IPSL)
UMR 1572 CEA-CNRS-UVSQ
Batiment 712 - Pe 119
91191 GIF Sur YVETTE Cedex

Tel: (33) 01 69 08 29 02; Fax:01.69.08.77.16


Hi Chao,

The warning you are getting:

WARNING: x coordinate not monotonically increasing - contour plot
may not be what you expect. If it looks odd, your can either
adjust the map projection region to be consistent with your data, or
(if your data is on a global lat/lon grid) use the shiftgrid
function to adjust the data to be consistent with the map projection
region (see examples/contour_demo.py).

Is important here. It looks like the x coordinate is not in appropriate
longitudes. Printing the first 5 and last 5 longitudes gives us our first
clue:

First 5: [-180. -178.99720764 -177.99443054 -176.99163818 -175.98886108]
Last 5 : [ 175.98886108 176.9916687 177.9944458 178.9972229 180.00003052]

Notice that the last longitude wraps around beyond 180. So if we were to
clip these numbers to -180 and +180 we will see that the warning disappears
and the contour is correct. This can be achieved with:

lon = np.clip(lon, -180, 180)

Alternatively, we can just construct the latitudes and longitudes ourselves
directly with:

lon, lat = np.meshgrid(np.linspace(-180, 180, 360), np.linspace(-90, 90,
180))

Incidentally, I tried these numbers with cartopy which has been designed to
handle dateline wrapping automatically, and the contour worked with the
unmodified longitudes (http://nbviewer.ipython.org/gist/pelson/10830039).

···

---------------------------

@JeffWhitaker - This looks like a bug with float tolerances in the makegrid
function. It currently does:

    def makegrid(self,nx,ny,returnxy=False):
        dx = (self.urcrnrx-self.llcrnrx)/(nx-1)
        dy = (self.urcrnry-self.llcrnry)/(ny-1)

But might be better if it did:

   def makegrid(self,nx,ny,returnxy=False):
        x = np.linspace(self.llcrnrx, self.urcrnrx, nx)
        y = np.linspace(self.llcrnry, self.urcrnry, ny)

To avoid the multiplicative floating point drift that is currently being
seen.

HTH,

Phil

Dear Phil,

Thank you. This solves my problem. So the title of my mail is wrong, the behaviour is reasonable but I am using wrong coordinates.

And also thanks to Jeff.

Cheers,

Chao

···

On Wed, Apr 16, 2014 at 10:14 AM, Phil Elson <pelson.pub@…287…> wrote:


please visit:
http://www.globalcarbonatlas.org/


Chao YUE
Laboratoire des Sciences du Climat et de l’Environnement (LSCE-IPSL)
UMR 1572 CEA-CNRS-UVSQ
Batiment 712 - Pe 119
91191 GIF Sur YVETTE Cedex

Tel: (33) 01 69 08 29 02; Fax:01.69.08.77.16


Hi Chao,

The warning you are getting:

WARNING: x coordinate not monotonically increasing - contour plot

may not be what you expect. If it looks odd, your can either
adjust the map projection region to be consistent with your data, or
(if your data is on a global lat/lon grid) use the shiftgrid
function to adjust the data to be consistent with the map projection

region (see examples/contour_demo.py).

Is important here. It looks like the x coordinate is not in appropriate longitudes. Printing the first 5 and last 5 longitudes gives us our first clue:

First 5: [-180. -178.99720764 -177.99443054 -176.99163818 -175.98886108]
Last 5 : [ 175.98886108 176.9916687 177.9944458 178.9972229 180.00003052]

Notice that the last longitude wraps around beyond 180. So if we were to clip these numbers to -180 and +180 we will see that the warning disappears and the contour is correct. This can be achieved with:

lon = np.clip(lon, -180, 180)

Alternatively, we can just construct the latitudes and longitudes ourselves directly with:

lon, lat = np.meshgrid(np.linspace(-180, 180, 360), np.linspace(-90, 90, 180))

Incidentally, I tried these numbers with cartopy which has been designed to handle dateline wrapping automatically, and the contour worked with the unmodified longitudes (http://nbviewer.ipython.org/gist/pelson/10830039).


@JeffWhitaker - This looks like a bug with float tolerances in the makegrid function. It currently does:

def makegrid(self,nx,ny,returnxy=False):

    dx = (self.urcrnrx-self.llcrnrx)/(nx-1)
    dy = (self.urcrnry-self.llcrnry)/(ny-1)

But might be better if it did:

def makegrid(self,nx,ny,returnxy=False):

    x = np.linspace(self.llcrnrx, self.urcrnrx, nx)
    y = np.linspace(self.llcrnry, self.urcrnry, ny)

To avoid the multiplicative floating point drift that is currently being seen.

HTH,

Phil