From: Sebastian Krieger <sebastian@io...>  20080531 02:27:42
Attachments:
Message as HTML

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.xmaxm.xmin)/dx)+1; ny = int((m.ymaxm.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 Jason1 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 
From: Jeff Whitaker <jswhit@fa...>  20080531 11:50:12

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.xmaxm.xmin)/dx)+1; ny = int((m.ymaxm.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 nonrectangular, 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 Jason1 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)4976313 NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 325 Broadway Boulder, CO, USA 803053328 