|
From: Jeff W. <js...@fa...> - 2010-04-02 12:16:56
|
from mpl_toolkits.basemap import Basemap, shiftgrid
import numpy as np
import matplotlib.pyplot as plt
# read in topo data (on a regular lat/lon grid)
# longitudes go from 20 to 380.
topoin = np.loadtxt('etopo20data.gz')
lons = np.loadtxt('etopo20lons.gz')
lats = np.loadtxt('etopo20lats.gz')
# shift data so lons go from -180 to 180 instead of 20 to 380.
topoin,lons = shiftgrid(180.,topoin,lons,start=False)
m = Basemap(projection='ortho',lon_0=-105,lat_0=40,resolution='l')
# transform to nx x ny regularly spaced native projection grid
nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
topodat,x,y =\
m.transform_scalar(topoin,lons,lats,nx,ny,returnxy=True,masked=True,order=1)
# create the figure.
fig=plt.figure(figsize=(8,8))
im = m.pcolormesh(x,y,topodat,cmap=plt.cm.jet)
m.drawcoastlines()
m.drawparallels(np.arange(0.,80,20.))
m.drawmeridians(np.arange(10.,360.,30.))
m.drawmapboundary()
plt.show()
|