On 4/2/10 4:27 AM, Will Hewson wrote:
Hi forum/ mailing list, When I plot in the orthographic projection I'm
getting the large artefact shown below extending away from the north
east of the globe. I'm not finding the same problem when plotting in a
full globe projection so I'm presuming the problem is with the way I'm
projecting everything rather than my data itself. I've included my
plotting code below, if anyone is able to spot some glaring omissions/
errors I'd be most grateful (I've been using python/ matplotlib for
only a couple of weeks now!).
Will: I think what's happening is that pcolormesh is having trouble
dealing with the higher curvlinear grid, which becomes nearly
pathological near the horizon of the projection. If you take a look at
the test.py file in the basemap examples directory, you'll see an
example orthographic plot that solves this problem by first
interpolating the data to a regular grid in projection coordinates (with
values over the plot horizon masked). The example uses imshow, but
pcolormesh works as well. A standalone version of the example using
pcolormesah is attached, which uses data files in the basemap examples
directory.
-Jeff
#!/usr/local/bin/python2.6
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import sys, glob
#input must be 3 col file of lons lats and data
#bins input values into half degree grid, ignores negative values
plts = glob.glob('*.plt')
x = np.arange(-180, 180, 0.5); y = np.arange(-90, 90, 0.5)
grid_lon, grid_lat = np.meshgrid(x,y) #regularly spaced 2D grid
n_vals = np.zeros((360,720)) #mean divisor
dat = np.zeros((360,720)) #2D grid of zeros
for pt in plts:
in_file = pt
data = np.loadtxt(in_file, comments = ';')
fname = in_file.split('.')[0]
lon = data[:,0] #original 1D list
lat = data[:,1] #original 1D list
slcol = data[:,2] #z data
lon = (np.around(lon*2))/2 #round to nearest .0 or 0.5
lat = (np.around(lat*2))/2 #round to nearest .0 or 0.5
##keep the below between files
j=0
for i in slcol:
if lon[j] < 0:
grid_lon_ind = 360+(lon[j]*2)
grid_lat_ind = 180+(lat[j]*2)
else:
grid_lon_ind = 360-(lon[j]*2)
grid_lat_ind = 180+(lat[j]*2)
if i > 0:
dat[grid_lat_ind, grid_lon_ind] += i #add i'th value
n_vals[grid_lat_ind, grid_lon_ind] += 1 #increase cell counter by 1
for each extra value
j+=1
dat = np.nan_to_num(dat/n_vals)
#create map object
fig = plt.figure()
m = Basemap(projection='ortho', lon_0=lon[(len(lon)/2)], lat_0=0,
resolution='l', area_thresh=10000.)
#m = Basemap(projection='moll',lon_0=0,resolution='c', area_thresh=10000.)
X,Y = m(grid_lon, grid_lat)
#pass all 2d arrays to pcolor
im = m.pcolormesh(X,Y,dat)
#add coastlines, globe boundary and colourbar
m.drawcoastlines()
m.drawmapboundary()
m.drawparallels(np.arange(-90, 90,30))
m.drawmeridians(np.arange(-180,180,30))
fig.colorbar(im)
plt.title('CH20 and ting')
plt.savefig('binplot.png')
Thanks for your help,
Will.
------------------------------------------------------------------------
View this message in context: Basemap/ orthographic projection plot
doesn't respect globe boundary
<http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28117654.html>
Sent from the matplotlib - users mailing list archive
<http://old.nabble.com/matplotlib---users-f2906.html> at Nabble.com.
------------------------------------------------------------------------------
Download Intel® Parallel Studio Eval
Try the new software tools for yourself. Speed compiling, find bugs
proactively, and fine-tune applications for parallel performance.
See why Intel Parallel Studio got high marks during beta.
http://p.sf.net/sfu/intel-sw-dev
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
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()
------------------------------------------------------------------------------
Download Intel® Parallel Studio Eval
Try the new software tools for yourself. Speed compiling, find bugs
proactively, and fine-tune applications for parallel performance.
See why Intel Parallel Studio got high marks during beta.
http://p.sf.net/sfu/intel-sw-dev
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
