matplotlib-users Mailing List for matplotlib

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

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