Re: [Matplotlib-users] Basemap/ orthographic projection plot doesn't respect globe boundary

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

On 4/2/10 6:32 AM, Will Hewson wrote:
> 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
>    

Will:  I forgot to mention that contourf will work on your data without 
having to interpolate to projection coordinates.

-Jeff
>
>
> 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()
>>
>>
>>      
>