From: <gruben@bi...>  20070601 04:19:12

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/basemapstereographicprojectionboundingboxestf2170166.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 Wulffnet 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() 
From: Jeff Whitaker <jswhit@fa...>  20070601 11:49:05

gruben@... wrote: > 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 ab= le to use the basemap toolkit's "stere" plot, centred at lat_0=3D0, lon_0= =3D0, with limits set at +/90deg lat and +/180deg lon, but we keep gett= ing tracebacks and I wondered whether this is possible, based on a commen= t from Jeff Whitaker <http://www.nabble.com/basemapstereographicproje= ctionboundingboxestf2170166.html#a6000978> which implies that basemap'= s stereographic projection code can't handle these default limits. Is thi= s the case? Would it be asking too much to request a small sample of gene= rating some Wulffnet 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 =3D figure(figsize=3D(4,4)) > # setup map projection > m =3D Basemap(projection=3D'stere',lat_0=3D0.,lon_0=3D0.,llcrnrlat=3D5= 0.,llcrnrlon=3D120., urcrnrlat=3D90., urcrnrlon=3D90.) > =20 > ax =3D 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=3D25483988.0 map =3D Basemap(lon_0=3D0,lat_0=3D0,projection=3D'stere',width=3Dw,height= =3Dw) print map.llcrnrlat,map.llcrnrlon,map.urcrnrlat,map.urcrnrlon map.drawparallels(p.arange(80,81,10)) map.drawmeridians(p.arange(90,91,10)) ax =3D p.gca() circ =3D 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=20 circular region by filling the region between the circle and the=20 bounding square. HTH, Jeff =20 Jeffrey S. Whitaker Phone : (303)4976313 NOAA/OAR/CDC R/PSD1 FAX : (303)4976449 325 Broadway Boulder, CO, USA 803053328 
From: Gary Ruben <gruben@bi...>  20070601 14:03:41

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@... wrote: >> 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/basemapstereographicprojectionboundingboxestf2170166.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 Wulffnet 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 