Re: [Matplotlib-users] contour and basemap

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

Evan Mason wrote:
> In the mpl basemap 'test.py' script, I want to add some contours to 
> the mercator projection map (test #3).  Just below line 83 (ie, below 
> the 'im  = imshow..' command I added the line:
>
> m.contour(topodat,[-1000, -2000])
>
> This returns:
>
> /home/evan/downloads/basemap-0.9.4/examples/test.py
>      82 # plot image over map.
>      83 im = m.imshow(topodat,cm.jet)
> ---> 84 m.contour(topodat,[-1000, -2000])
>      85 m.drawcoastlines ()
>      86 m.drawcountries()
>
> TypeError: contour() takes at least 4 arguments (3 given)
> WARNING: Failure executing file: <test.py>
>
> I understand the error to mean that I haven't supplied x and y values; 
> however at this point, if I close all the open figures and enter the 
> line 'contour(topodat,[-1000, -2000])' at the ipython command, it 
> gives me the plot I want.
>
> I've then tried to use meshgrid to make x and y values (see code 
> below); this is error free but the contours don't appear.  Again, 
> using contour in the command line I get the plot I want.  Is it not 
> possible to plot on top of imshow?
>
> Any help appreciated.
>
>
> Thanks, Evan
>
>
>
>
> from matplotlib.toolkits.basemap import Basemap, shiftgrid
> from scipy import interpolate as Pinterp
> from pylab import *
> import matplotlib.colors as colors
>
> def doSpline(yVec, newLatRange):
>     latRange    = range(len(yVec))
>     newLatRange = linspace(latRange[0], latRange[-1], newLatRange)
>     tck         = Pinterp.splrep(latRange, yVec)
>     yVec        = Pinterp.splev(newLatRange, tck)
>     return yVec
>
> topodatin = load('etopo20data.gz')
> lonsin = load('etopo20lons.gz')
> latsin = load('etopo20lats.gz')
> # shift data so lons go from -180 to 180 instead of 20 to 380.
> topoin,lons = shiftgrid(180.,topodatin,lonsin,start=False)
> lats = latsin
> # transform to nx x ny regularly spaced native projection grid
> nx = len(lons)
> ny = int(80.*len(lats)/90.)
> lats_ny = doSpline(lats, ny)
>
> lons_mesh, lats_mesh = meshgrid(lons, lats_ny)
>
> # setup mercator map projection (-80 to +80).
> m = Basemap(llcrnrlon=-180.,llcrnrlat=-80,urcrnrlon=180.,urcrnrlat=80.,\
>             resolution='c',area_thresh=10000.,projection='merc',\
>             lon_0=0.5*(lons[0]+lons[-1]),lat_ts=20.)
>
> topodat = m.transform_scalar(topoin,lons,lats,nx,ny)
>
> # setup figure with same aspect ratio as map.
> fig=figure()
> fig.add_axes([0.1,0.1,0.75,0.75])
> # plot image over map.
> im = m.imshow(topodat,cm.jet)
> # PLOT CONTOURS
> m.contour(lons_mesh, lats_mesh, topodat,[-1000, -2000, -3000])
> show()
>

Evan:  This works for me:

from matplotlib.toolkits.basemap import Basemap, shiftgrid
from pylab import *

topodatin = load('etopo20data.gz')
lonsin = load('etopo20lons.gz')
latsin = load('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topodatin,lonsin,start=False)
lats = latsin

# setup mercator map projection (-80 to +80).
m = Basemap(llcrnrlon=-180.,llcrnrlat=-80,urcrnrlon=180.,urcrnrlat=80.,\
            resolution='c',area_thresh=10000.,projection='merc',\
            lon_0=0.5*(lons[0]+lons[-1]),lat_ts=20.)

# transform to nx x ny regularly spaced native projection grid
nx = len(lons); ny = len(lats)
topodat = m.transform_scalar(topoin,lons,lats,nx,ny)

# setup figure with same aspect ratio as map.
# plot image over map.
im = m.imshow(topodat,cm.jet)
# PLOT CONTOURS
longrid,latgrid = meshgrid(lons,lats)
x, y  = m(longrid,latgrid)
# reverse sign of data, so contours won't be dashed.
cs = m.contour(x, y, -topoin, [1000, 2000], colors='r')
m.drawcoastlines()
show()


Unlike the pylab contour command, the basemap contour method requires 
you to supply the X and Y values.

-Jeff

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