Inside / outside judgment with Python ➁: About donut-shaped polygons
Thing you want to do
For the following donut-shaped polygon shp files, we have organized how to judge the inside and outside of any point.
(Source: https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A31-v2_1.html)
I will explain in the following order. Step1: Read donut-shaped polygon data Step2: Inside / outside judgment module
Reading donut-shaped polygon data
In the shp file, the donut-shaped polygon is constructed by superimposing multiple uneven polygons (parts). Therefore, first we need to extract the nodal coordinates of each part.
import numpy
import shapefile
import matplotlib.pyplot as plt
#
def outLS(pnt,fig): #Draw the line data that connects the node data of the polygon
print(*pnt)
xs,ys = zip(*pnt)
plt.plot(xs,ys)
return 0
#
fname="test.shp" #Input file
src=shapefile.Reader(fname,encoding='SHIFT-JIS') #File reading
SRS=src.shapes() #Get feature information
#
for srs in SRS:
IP=srs.parts #Acquire parts for various parts
NP=len(IP) #Get the number of parts
pnt=srs.points #Get node data
pnts=[] #Organize the nodal coordinates of parts for each part
for ip in range(NP-1):
pnts.append(pnt[IP[ip]:IP[ip+1]])
pnts.append(pnt[IP[NP-1]:])
#
fig = plt.figure() #Draw outline information for each part
for np in range(NP):
outLS(pnts[np],fig)
plt.show()
plt.clf()
Execution result
As shown below, the outline information of each part that makes up the donut-shaped polygon could be obtained separately.
Inside / outside judgment module
It's actually easy to judge inside and outside. __ You just have to check the number of intersections between the half line starting from the target point and the outline of all the parts that make up the donut-shaped polygon. If the number of intersections is odd, it will be judged inside, and if it is even, it will be judged outside __. The source code looks like the following. In the example below, the number of intersections of the half line and the outline drawn from the target point toward the east is counted.
import sys
import numpy as np
import shapefile
import matplotlib.pyplot as plt
#
def naig(long,lat,pnts,boxs):
#Inside / outside judgment:Half straight line eastward from the target point(=A)Counts the number of intersections of a straight line and the outline of a feature when
NPS=len(pnts)
ccount=0
for nps in range(NPS):
pnt=pnts[nps]
box=boxs[nps]
#
if long<=box[2] and box[1]<=lat<=box[3]: #Use rectangles to narrow down the parts that may intersect
NP=len(pnt)
for np in range(NP-1):
xx1=pnt[np][0] #Line segment start point X
yy1=pnt[np][1] #Line segment start point Y
xx2=pnt[np+1][0] #End point X of line segment
yy2=pnt[np+1][1] #End of line segment Y
miny=min([yy1,yy2]) #Minimum Y of line segment
maxy=max([yy1,yy2]) #Maximum Y of line segment
minx=min([xx1,xx2]) #Minimum line segment X
maxx=max([xx1,xx2]) #Maximum line segment X
if abs(xx1-xx2)<1.E-7: #The line segment is vertical in the north-south direction
if xx1>long and miny<=lat<maxy: #Target point Y coordinate is greater than or equal to the minimum Y of the line segment and less than the maximum Y
ccount=ccount+1
elif abs(yy1-yy2)<1.E-7: #Line segment is almost horizontal
ccount=ccount #No intersection judgment required for horizontal line segments
else:
aa=(yy2-yy1)/(xx2-xx1) #Slope of a straight line through a line segment
bb=yy2-aa*xx2 #Intercept of a straight line through a line segment
yc=lat #Intersection Y coordinate of straight line and A
xc=(yc-bb)/aa #X coordinate of the intersection of a straight line and A
#The intersection X coordinate must be greater than the target point X coordinate, and the intersection Y coordinate must be greater than or equal to the minimum Y and less than the maximum Y of the line segment.
if xc>long and miny<=yc<maxy:
ccount=ccount+1
return ccount
#
fname="test.shp" #Input file
src=shapefile.Reader(fname,encoding='SHIFT-JIS') #File reading
SRS=src.shapes() #Get feature information
SRR=src.shapeRecords() #Get feature information
#
#Get feature information
pntall=[]
for srs in SRS:
IP=srs.parts #Acquire parts for various parts
NP=len(IP) #Get the number of parts
pnt=srs.points #Get node data
pnts=[] #Organize the nodal coordinates of parts for each part
for ip in range(NP-1):
pnts.append(pnt[IP[ip]:IP[ip+1]])
pnts.append(pnt[IP[NP-1]:])
pntall.append(pnts)
#
#Get attribute information
recall=[]
for srr in SRR:
recall.append(srr.record)
if len(pntall)!=len(recall):
print("Mismatch between the number of attribute information and the number of features!!")
sys.exit()
#
#Find the rectangle that encloses each part(Rectangle information acquisition)
NPA=len(pntall) #Number of features
boxall=[]
for npa in range(NPA):
pnts=pntall[npa]
NPS=len(pnts) #Number of parts
boxs=[]
for nps in range(NPS): #Get the minimum and maximum values of latitude and longitude of the nodes of each part
pp=np.array(list(pnts[nps]),dtype=np.float32)
bbox=[pp[:,0].min(),pp[:,1].min(),pp[:,0].max(),pp[:,1].max()]
boxs.append(bbox)
boxall.append(boxs)
#
#Count the number of intersections
LON=[130.60533882626782543,130.59666405110618825,130.60918282680634661,130.60550819793903088,130.60379578346410767]
LAT=[32.76515635424413375,32.77349238606328896,32.77375748954870716,32.76751282967006063,32.77819292046771693]
NPP=len(LON) #Number of points
for npp in range(NPP): #Loop about the point
for npa in range(NPA): #Loop about features
ic=naig(LON[npp],LAT[npp],pntall[npa],boxall[npa])
if ic % 2==0: #If on the outside
print("point{0}Is a feature{1}Outside! (Number of crossings:{2})".format(npp+1,npa+1,ic))
else: #If inside
print("point{0}Is a feature{1}in! (Number of crossings:{2})".format(npp+1,npa+1,ic))
Execution result I was able to judge as follows. It seems that this self-made module is much faster than the method using sympy for internal / external judgment.
Point 1 is outside feature 1! (Number of crossings:14)
Point 2 is outside feature 1! (Number of crossings:18)
Point 3 is in feature 1! (Number of crossings:3)
Point 4 is outside feature 1! (Number of crossings:4)
Point 5 is in feature 1! (Number of crossings:1)
Recommended Posts