geospatialpython/PiP_Edge.py at master · GeospatialPython/geospatialpython · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# Improved point in polygon test which includes edge
# and vertex points

defpoint_in_poly(x,y,poly):

# check if point is a vertex
if (x,y) inpoly: return"IN"

# check if point is on a boundary
foriinrange(len(poly)):
p1=None
p2=None
ifi==0:
p1=poly[0]
p2=poly[1]
else:
p1=poly[i-1]
p2=poly[i]
ifp1[1] ==p2[1] andp1[1] ==yandx>min(p1[0], p2[0]) andx<max(p1[0], p2[0]):
return"IN"

n=len(poly)
inside=False

p1x,p1y=poly[0]
foriinrange(n+1):
p2x,p2y=poly[i%n]
ify>min(p1y,p2y):
ify<=max(p1y,p2y):
ifx<=max(p1x,p2x):
ifp1y!=p2y:
xints= (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
ifp1x==p2xorx<=xints:
inside=notinside
p1x,p1y=p2x,p2y

ifinside: return"IN"
else: return"OUT"

# Test a vertex for inclusion
poligono= [(-33.416032,-70.593016), (-33.415370,-70.589604),
(-33.417340,-70.589046), (-33.417949,-70.592351),
(-33.416032,-70.593016)]
lat=-33.416032
lon=-70.593016

printpoint_in_poly(lat, lon, poligono)

# test a boundary point for inclusion
poly2= [(1,1), (5,1), (5,5), (1,5), (1,1)]
x=3
y=1
printpoint_in_poly(x, y, poly2)