This is great Jeff, thanks for the help - I'll give it a try over the weekend
(it's bank holiday here in the UK!) and get back to you, if I'm still having
trouble I'll stick up the plotting data too... thanks again.
Will
Jeff Whitaker wrote:
>
> 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
>
> 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()
>
>
--
View this message in context:
http://old.nabble.com/Basemap--orthographic-projection-plot-doesn%27t-respect-globe-boundary-tp28117654p28118555.html
Sent from the matplotlib - users mailing list archive 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
