Hi again,
still trying to accomplish the filling of certain paths only inside U.S. It ocurred to me that set_clip_path could work but so far, I have not been able to show the curve in the map.
Could someone point out what is missing? Thanks
#!/usr/bin/env python
"""http://www.scipy.org/Cookbook/Finding_Convex_Hull"""
import numpy as n, pylab as p, time
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.path as mpath
from matplotlib.patches import Polygon
from mpl_toolkits.basemap import Basemap
import sys
import numpy as np
import matplotlib.cm as cm
def _angle_to_point(point, centre):
'''calculate angle in 2-D between points and x axis'''
delta = point - centre
res = n.arctan(delta[1] / delta[0])
if delta[0] < 0:
res += n.pi
return res
def _draw_triangle(p1, p2, p3, **kwargs):
tmp = n.vstack((p1,p2,p3))
x,y = [x[0] for x in zip(tmp.transpose())]
p.fill(x,y, **kwargs) # commented out to see result
def area_of_triangle(p1, p2, p3):
'''calculate area of any triangle given co-ordinates of the corners'''
return n.linalg.norm(n.cross((p2 - p1), (p3 - p1)))/2.
def convex_hull(points, graphic=True, smidgen=0.0075):
'''Calculate subset of points that make a convex hull around points
Recursively eliminates points that lie inside two neighbouring points until only convex hull is remaining.
:Parameters:
points : ndarray (2 x m)
array of points for which to find hull
graphic : bool
use pylab to show progress?
smidgen : float
offset for graphic number labels - useful values depend on your data range
:Returns:
hull_points : ndarray (2 x n)
convex hull surrounding points
'''
# Commenting these 2 lines out shows the plotting of the country border
# if graphic:
# p.clf()
n_pts = points.shape[1]
print "POINTS shape[1]:%d" %(n_pts)
assert(n_pts > 5)
centre = points.mean(1)
angles = n.apply_along_axis(_angle_to_point, 0, points, centre)
pts_ord = points[:,angles.argsort()]
pts = [x[0] for x in zip(pts_ord.transpose())]
prev_pts = len(pts) + 1
k = 0
while prev_pts > n_pts:
prev_pts = n_pts
n_pts = len(pts)
if graphic: p.gca().patches = []
i = -2
while i < (n_pts - 2):
Aij = area_of_triangle(centre, pts[i], pts[(i + 1) % n_pts])
Ajk = area_of_triangle(centre, pts[(i + 1) % n_pts], \
pts[(i + 2) % n_pts])
Aik = area_of_triangle(centre, pts[i], pts[(i + 2) % n_pts])
if Aij + Ajk < Aik:
if graphic: p.plot((pts[i + 1][0],),(pts[i + 1][1],)) # without green circle marker
del pts[i+1]
i += 1
n_pts = len(pts)
k += 1
return n.asarray(pts)
if __name__ == "__main__":
x1 = -124.
x2 = -57.
y1 = 18.
y2 = 51.
m = Basemap(resolution='i',projection='stere', lat_0=39., lon_0=-98.,llcrnrlat=y1,urcrnrlat=y2,llcrnrlon=x1,urcrnrlon=x2)
point0 = []
point1 = []
cpaths = (m.drawcountries()).get_paths()
for cpath in cpaths:
vs = cpath.vertices
for v in vs:
x=v[0]
y=v[1]
point0.append(x)
point1.append(y)
nonumpy = []
nonumpy.append(point0)
nonumpy.append(point1)
points = n.array(nonumpy)
hull_pts = convex_hull(points)
poly = Polygon(hull_pts)
USpath = poly.get_path()
USpatch = mpatches.PathPatch(USpath,facecolor='none',edgecolor='none')
ax = plt.gca()
ax.add_patch(USpatch)
# Add my curve
Path = mpath.Path
pathdata = []
pathdata.append((Path.MOVETO,m(-120.,35.)))
pathdata.append((Path.LINETO,m(-100.,35.)))
pathdata.append((Path.LINETO,m(-100.,30.)))
pathdata.append((Path.LINETO,m(-90.,28.)))
pathdata.append((Path.LINETO,m(-90,25)))
pathdata.append((Path.LINETO,m(-100,25)))
pathdata.append((Path.CURVE3,m(-120.,35.)))
codes,verts = zip(*pathdata)
im = ax.imshow(verts,cmap=cm.gray,clip_path=USpatch,clip_on=True)
im.set_clip_path(USpatch)
m.drawparallels(np.arange(20,60,5),labels=[1,0,0,0],color='black')
m.drawmeridians(np.arange(-140,60,10),labels=[0,0,0,1],color='black')
plt.gcf().savefig('./ch.png')