Re: [Matplotlib-users] basemap and fill_between()

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

On 5/13/12 3:34 AM, David Craig wrote:
> Hi, I'm having a problem usinf fill_between() with basemap. I plot two 
> great circles and want to shade the region between them. My code is 
> below, it doesnt give any error just creates the plot without filling 
> the area. Does anyone know if it's possible to do this or should I try 
> a different method?
> Thanks,
> David
>

David Try replacing

m.drawgreatcircle(x1, y1, x2, y2, del_s=10, color='gray', lw=1.)
m.drawgreatcircle(x1, y1, x3, y3, del_s=10, color='gray', lw=1.)
a=linspace(x3,x1)
b=linspace(y2,y1)
c=linspace(y3,y1)
fill_between(a, b, c, where=None, alpha=0.2)

with

xx1,yy1 = m.gcpoints(x1,y1,x2,y2,100)
xx2,yy2 = m.gcpoints(x1,y1,x3,y3,100)
fill_between(xx1, yy1, yy2)

-Jeff


> from mpl_toolkits.basemap import Basemap
> from pylab import *
>
> ### PARAMETERS FOR MATPLOTLIB :
> import matplotlib as mpl
> rcParams['font.size'] = 10.
> rcParams['font.family'] = 'Comic Sans MS'
> rcParams['axes.labelsize'] = 8.
> rcParams['xtick.labelsize'] = 6.
> rcParams['ytick.labelsize'] = 6.
>
> def shoot(lon, lat, azimuth, maxdist=None):
>     """Shooter Function
>     Original javascript on http://williams.best.vwh.net/gccalc.htm
>     Translated to python by Thomas Lecocq
>     """
>     glat1 = lat * pi / 180.
>     glon1 = lon * pi / 180.
>     s = maxdist / 1.852
>     faz = azimuth * pi / 180.
>
>     EPS= 0.00000000005
>     if ((abs(cos(glat1))<EPS) and not (abs(sin(faz))<EPS)):
>         alert("Only N-S courses are meaningful, starting at a pole!")
>
>     a=6378.13/1.852
>     f=1/298.257223563
>     r = 1 - f
>     tu = r * tan(glat1)
>     sf = sin(faz)
>     cf = cos(faz)
>     if (cf==0):
>         b=0.
>     else:
>         b=2. * arctan2 (tu, cf)
>
>     cu = 1. / sqrt(1 + tu * tu)
>     su = tu * cu
>     sa = cu * sf
>     c2a = 1 - sa * sa
>     x = 1. + sqrt(1. + c2a * (1. / (r * r) - 1.))
>     x = (x - 2.) / x
>     c = 1. - x
>     c = (x * x / 4. + 1.) / c
>     d = (0.375 * x * x - 1.) * x
>     tu = s / (r * a * c)
>     y = tu
>     c = y + 1
>     while (abs (y - c) > EPS):
>
>         sy = sin(y)
>         cy = cos(y)
>         cz = cos(b + y)
>         e = 2. * cz * cz - 1.
>         c = y
>         x = e * cy
>         y = e + e - 1.
>         y = (((sy * sy * 4. - 3.) * y * cz * d / 6. + x) *
>               d / 4. - cz) * sy * d + tu
>
>     b = cu * cy * cf - su * sy
>     c = r * sqrt(sa * sa + b * b)
>     d = su * cy + cu * sy * cf
>     glat2 = (arctan2(d, c) + pi) % (2*pi) - pi
>     c = cu * cy - su * sy * cf
>     x = arctan2(sy * sf, c)
>     c = ((-3. * c2a + 4.) * f + 4.) * c2a * f / 16.
>     d = ((e * cy * c + cz) * sy * c + y) * sa
>     glon2 = ((glon1 + x - (1. - c) * d * f + pi) % (2*pi)) - pi
>
>     baz = (arctan2(sa, b) + pi) % (2 * pi)
>
>     glon2 *= 180./pi
>     glat2 *= 180./pi
>     baz *= 180./pi
>
>     return (glon2, glat2, baz)
>
> #Create a basemap around N. Atlantic
> m = Basemap(llcrnrlon=-45.0,llcrnrlat=30.0,urcrnrlon=15.0,urcrnrlat=75.0,
>             resolution='i',projection='merc',lon_0=-17.5,lat_0=60.0)
>
>
> m.drawcountries(linewidth=0.5)
> m.drawcoastlines(linewidth=0.5)
> m.bluemarble()
> m.drawparallels(arange(40.,75.,10.),labels=[1,0,0,0],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) 
> # draw parallels
> m.drawmeridians(arange(-45.,15.,10.),labels=[0,0,0,1],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) 
> # draw meridians
>
> # Shade region defined by great circles.
> x1, y1 = -9.1676613, 51.6029999
> az1 = 270.
> az2 = 290.
> maxdist = 2000
> x2, y2, baz = shoot(x1, y1, az1, maxdist)
> x3, y3, baz = shoot(x1, y1, az2, maxdist)
>
> m.drawgreatcircle(x1, y1, x2, y2, del_s=10, color='gray', lw=1.)
> m.drawgreatcircle(x1, y1, x3, y3, del_s=10, color='gray', lw=1.)
> a=linspace(x3,x1)
> b=linspace(y2,y1)
> c=linspace(y3,y1)
> fill_between(a, b, c, where=None, alpha=0.2)
> show()
>
>
> ------------------------------------------------------------------------------
> Live Security Virtual Conference
> Exclusive live event will cover all the ways today's security and
> threat landscape has changed and how IT managers can respond. Discussions
> will include endpoint security, mobile security and the latest in malware
> threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
>
>
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users