Re: [Matplotlib-users] Robinson projection imshow out of bounds and different latitude limits

[Jeff W. <js...@fa...>]

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