Robinson projection imshow out of bounds and different latitude limits

Dear all!

I would like to ask two questions: one concerning imshow with
the Robinson projection and the second about the latitude limits for
the same map projection.

First, I am trying to make a Robinson projection map using imshow
instead of contourf as described in the Basemap
documentation. I modified the script as follows:

#----------------------------------------------------------------------

from matplotlib.toolkits.basemap import Basemap, shiftgrid

from pylab import load, meshgrid, title, arange, show, cm, pi

#

# read in topo data (on a regular lat/lon grid)

etopo = load('etopo20data.gz')

lons = load('etopo20lons.gz')

lats = load('etopo20lats.gz')

#

# create Basemap instance for Robinson projection.

m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))

#

# compute native map projection coordinates for lat/lon grid.

etopo, lons = shiftgrid(180., etopo, lons, start=False)

x, y = m(*meshgrid(lons,lats))

dx = 2.*pi*m.rmajor/len(lons)

nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1

dat, x, y = m.transform_scalar(etopo, lons, lats, nx, ny, returnxy=True)

#

# make filled contour plot.

im = m.imshow(dat, cmap=cm.jet)

m.drawcoastlines() # draw coastlines

m.drawmapboundary() # draw a line around the map region

m.drawparallels(arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels

m.drawmeridians(arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians

title('Robinson Projection') # add a title

show()

#----------------------------------------------------------------------

The result looks as I would expect, except that some data is plot
outside of the map boundaries. Using another projection as the
Orthographic for example, this problem doesn’t happen. Am I doing
something wrong?

Second, I work with Topex/POSEIDON and Jason-1 sea surface height
anomalies datasets where the latitudes range from about 67.5S to 67.5N.
Outside these limits hopefully its blank, as anyone would expect.
Aesthetically I find it more appealing if I could limit my map
boundaries to these limits, or even lower limits if I zoom the
equatorial region. Has anyone ever tried to do this?

Thank’s in advance and kind regards,

Sebastian

Sebastian Krieger wrote:

Dear all!

I would like to ask two questions: one concerning imshow with the Robinson projection and the second about the latitude limits for the same map projection.

First, I am trying to make a Robinson projection map using imshow instead of contourf as described in the Basemap documentation. I modified the script as follows:

    #----------------------------------------------------------------------
    from matplotlib.toolkits.basemap import Basemap, shiftgrid
    from pylab import load, meshgrid, title, arange, show, cm, pi
    #
    # read in topo data (on a regular lat/lon grid)
    etopo = load('etopo20data.gz')
    lons = load('etopo20lons.gz')
    lats = load('etopo20lats.gz')
    #
    # create Basemap instance for Robinson projection.
    m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
    #
    # compute native map projection coordinates for lat/lon grid.
    etopo, lons = shiftgrid(180., etopo, lons, start=False)
    x, y = m(*meshgrid(lons,lats))
    dx = 2.*pi*m.rmajor/len(lons)
    nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1
    dat, x, y = m.transform_scalar(etopo, lons, lats, nx, ny,
    returnxy=True)
    #
    # make filled contour plot.
    im = m.imshow(dat, cmap=cm.jet)
    m.drawcoastlines() # draw coastlines
    m.drawmapboundary() # draw a line around the map region
    m.drawparallels(arange(-90.,120.,30.),labels=[1,0,0,0]) # draw
    parallels
    m.drawmeridians(arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians
    title('Robinson Projection') # add a title
    show()
    #----------------------------------------------------------------------

The result looks as I would expect, except that some data is plot outside of the map boundaries. Using another projection as the Orthographic for example, this problem doesn't happen. Am I doing something wrong?

Sebastian: The Robinson projection is non-rectangular, so it's not straightforward to plot an image on it. You would somehow have to mask the points outside the projection limb (this is what happens for the 'ortho' projection). It's much easier just to use pcolor, for example

from matplotlib.toolkits.basemap import Basemap, shiftgrid
from pylab import load, meshgrid, title, arange, show, cm, pi
# read in topo data (on a regular lat/lon grid)
etopo = load('etopo20data.gz')
lons = load('etopo20lons.gz')
lats = load('etopo20lats.gz')
# create Basemap instance for Robinson projection.
m = Basemap(projection='robin',lon_0=0.5*(lons[0]+lons[-1]))
x,y = m(*meshgrid(lons,lats))
p = m.pcolormesh(x,y,etopo,shading='flat')
m.drawcoastlines() # draw coastlines
m.drawmapboundary() # draw a line around the map region
m.drawparallels(arange(-90.,120.,30.),labels=[1,0,0,0]) # draw parallels
m.drawmeridians(arange(0.,420.,60.),labels=[0,0,0,1]) # draw meridians
title('Robinson Projection') # add a title
show()

Second, I work with Topex/POSEIDON and Jason-1 sea surface height anomalies datasets where the latitudes range from about 67.5S to 67.5N. Outside these limits hopefully its blank, as anyone would expect. Aesthetically I find it more appealing if I could limit my map boundaries to these limits, or even lower limits if I zoom the equatorial region. Has anyone ever tried to do this?

You can't with the Robinson projection, it's only defined globally. You could do it with the Mercator ('merc'), Miller ('mill') or Cylindrical Equidistant ('cyl') projections by specifying the lat/lon values of the upper right and lower left corners.

Thank's in advance and kind regards,
Sebastian

HTH,

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
NOAA/OAR/CDC R/PSD1 FAX : (303)497-6449
325 Broadway Boulder, CO, USA 80305-3328