basemap stereographic projection problem

I've been trying to help a friend who wants to plot directional data on a "Wulff net" <http://en.wikipedia.org/wiki/Pole_figure#Geometry_in_the_pole_figure>, which is a stereographic projection plot. She wants to plot points specified by latitude and longitude in degrees. We hoped to be able to use the basemap toolkit's "stere" plot, centred at lat_0=0, lon_0=0, with limits set at +/-90deg lat and +/-180deg lon, but we keep getting tracebacks and I wondered whether this is possible, based on a comment from Jeff Whitaker <http://www.nabble.com/-basemap--stereographic-projection-bounding-boxes-tf2170166.html#a6000978> which implies that basemap's stereographic projection code can't handle these default limits. Is this the case? Would it be asking too much to request a small sample of generating some Wulff-net axes? Plotting points on these seems simple enough. We started with the polarmaps.py example and the code below is as close as we could get. Any suggestions would be welcome.

Gary Ruben

···

--

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

# loop over projections, one for each panel of the figure.
fig = figure(figsize=(4,4))
# setup map projection
m = Basemap(projection='stere',lat_0=0.,lon_0=0.,llcrnrlat=-50.,llcrnrlon=-120., urcrnrlat=90., urcrnrlon=90.)
           
ax = fig.add_subplot(1,1,1)
# draw parallels and meridians.
m.drawparallels(arange(-180.,180.,10.))
m.drawmeridians(arange(-90.,90.,10.))
# draw boundary around map region.
m.drawmapboundary()
show()

gruben@...636... wrote:

I've been trying to help a friend who wants to plot directional data on a "Wulff net" <Pole figure - Wikipedia, which is a stereographic projection plot. She wants to plot points specified by latitude and longitude in degrees. We hoped to be able to use the basemap toolkit's "stere" plot, centred at lat_0=0, lon_0=0, with limits set at +/-90deg lat and +/-180deg lon, but we keep getting tracebacks and I wondered whether this is possible, based on a comment from Jeff Whitaker <http://www.nabble.com/-basemap--stereographic-projection-bounding-boxes-tf2170166.html#a6000978&gt; which implies that basemap's stereographic projection code can't handle these default limits. Is this the case? Would it be asking too much to request a small sample of generating some Wulff-net axes? Plotting points on these seems simple enough. We started with the polarmaps.py example and the code below is as close as we could get. Any suggestions would be welco
me.

Gary Ruben

--

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

# loop over projections, one for each panel of the figure.
fig = figure(figsize=(4,4))
# setup map projection
m = Basemap(projection='stere',lat_0=0.,lon_0=0.,llcrnrlat=-50.,llcrnrlon=-120., urcrnrlat=90., urcrnrlon=90.)
           ax = fig.add_subplot(1,1,1)
# draw parallels and meridians.
m.drawparallels(arange(-180.,180.,10.))
m.drawmeridians(arange(-90.,90.,10.))
# draw boundary around map region.
m.drawmapboundary()
show()

-

Gary:

How about this?

import pylab as p
from matplotlib.toolkits.basemap import Basemap
from matplotlib.patches import Circle, Polygon
w=25483988.0
map = Basemap(lon_0=0,lat_0=0,projection='stere',width=w,height=w)
print map.llcrnrlat,map.llcrnrlon,map.urcrnrlat,map.urcrnrlon
map.drawparallels(p.arange(-80,81,10))
map.drawmeridians(p.arange(-90,91,10))
ax = p.gca()
circ = Circle(map(0,0),0.5*w)
circ.set_fill(False)
circ.set_edgecolor('k')
circ.set_linewidth(1.0)
circ.set_clip_on(True)
ax.add_patch(circ)
#ax.set_frame_on(False)
p.show()

It should be possible to get rid of the lines extending outside the circular region by filling the region between the circle and the bounding square.

HTH,

-Jeff

···

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

Many thanks for the special case code Jeff,
This appears to work well.
I see you picked up on our confusion about the -180 to 180 longitude range. I'll pass this on and look at clipping the lines outside the circle later,
regards,
Gary

Jeff Whitaker wrote:

···

gruben@...636... wrote:

I've been trying to help a friend who wants to plot directional data on a "Wulff net" <Pole figure - Wikipedia, which is a stereographic projection plot. She wants to plot points specified by latitude and longitude in degrees. We hoped to be able to use the basemap toolkit's "stere" plot, centred at lat_0=0, lon_0=0, with limits set at +/-90deg lat and +/-180deg lon, but we keep getting tracebacks and I wondered whether this is possible, based on a comment from Jeff Whitaker <http://www.nabble.com/-basemap--stereographic-projection-bounding-boxes-tf2170166.html#a6000978&gt; which implies that basemap's stereographic projection code can't handle these default limits. Is this the case? Would it be asking too much to request a small sample of generating some Wulff-net axes? Plotting points on these seems simple enough. We started with the polarmaps.py example and the code below is as close as we could get. Any suggestions would be welco
me.

Gary Ruben

--

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

# loop over projections, one for each panel of the figure.
fig = figure(figsize=(4,4))
# setup map projection
m = Basemap(projection='stere',lat_0=0.,lon_0=0.,llcrnrlat=-50.,llcrnrlon=-120., urcrnrlat=90., urcrnrlon=90.)
           ax = fig.add_subplot(1,1,1)
# draw parallels and meridians.
m.drawparallels(arange(-180.,180.,10.))
m.drawmeridians(arange(-90.,90.,10.))
# draw boundary around map region.
m.drawmapboundary()
show()

-

Gary:

How about this?

import pylab as p
from matplotlib.toolkits.basemap import Basemap
from matplotlib.patches import Circle, Polygon
w=25483988.0
map = Basemap(lon_0=0,lat_0=0,projection='stere',width=w,height=w)
print map.llcrnrlat,map.llcrnrlon,map.urcrnrlat,map.urcrnrlon
map.drawparallels(p.arange(-80,81,10))
map.drawmeridians(p.arange(-90,91,10))
ax = p.gca()
circ = Circle(map(0,0),0.5*w)
circ.set_fill(False)
circ.set_edgecolor('k')
circ.set_linewidth(1.0)
circ.set_clip_on(True)
ax.add_patch(circ)
#ax.set_frame_on(False)
p.show()

It should be possible to get rid of the lines extending outside the circular region by filling the region between the circle and the bounding square.

HTH,

-Jeff