basemap and fill_between()

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

Hello, Boss.
Perhaps I am telling rubbish, but fill_between demands the coordinates within your axes, not the basemap angles. Your a,b,c contain some small numers (negative x, say -38) won't do any good, x should be, say 3001129 etc. Use m(...) to convert.

Jerzy K

···

Le 13/05/2012 11:34, David Craig a écrit :

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?

Appendix.
(and excuses for my approximate syntax in the first message; I was doing three things simultaneously).

Maestro David Craig, your filling will be lousy anyway.
If your x and y arrays come from linespace, you will surely get a triangle, and not the curvilinear area between GreatCircles. You must transform the coordinates in a more clever (but equally easy...) way.

Jerzy K.

``````  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
``````
···