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
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
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
