Innen- / Außenbeurteilung mit Python ➁: Über Donut-förmige Polygone

Was du machen willst

Für die folgenden Donut-förmigen Polygon-Shp-Dateien haben wir organisiert, wie das Innere und Äußere eines beliebigen Punkts beurteilt werden kann. (Quelle: https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A31-v2_1.html)

Ich werde in der folgenden Reihenfolge erklären. Schritt 1: Lesen Sie Donut-förmige Polygondaten Schritt 2: Innen- / Außenbewertungsmodul

Donutförmige Polygondaten lesen

In der shp-Datei wird das Donut-förmige Polygon durch Überlagerung mehrerer ungleichmäßiger Polygone (Teile) konstruiert. Daher ist es zunächst erforderlich, die Knotenkoordinaten jedes Teils zu extrahieren.

import numpy
import shapefile
import matplotlib.pyplot as plt
#
def outLS(pnt,fig):    #Zeichnen Sie die Liniendaten, die die Knotendaten des Polygons verbinden
    print(*pnt)
    xs,ys = zip(*pnt)
    plt.plot(xs,ys) 
    return 0
#
fname="test.shp"    #Eingabedatei
src=shapefile.Reader(fname,encoding='SHIFT-JIS')  #Datei lesen
SRS=src.shapes()    #Funktionsinformationen abrufen
#
for srs in SRS:
    IP=srs.parts    #Teile für verschiedene Teile erwerben
    NP=len(IP)      #Holen Sie sich die Anzahl der Teile
    pnt=srs.points  #Knotendaten abrufen
    pnts=[]         #Organisieren Sie die Knotenkoordinaten der Teile für jedes Teil
    for ip in range(NP-1):
        pnts.append(pnt[IP[ip]:IP[ip+1]])
    pnts.append(pnt[IP[NP-1]:])
    #
    fig = plt.figure() #Zeichnen Sie Umrissinformationen für jedes Teil
    for np in range(NP):
        outLS(pnts[np],fig)
    plt.show()
    plt.clf()

Ausführungsergebnis Wie unten gezeigt, konnten wir die Umrissinformationen jedes Teils unterscheiden, aus dem das donutförmige Polygon besteht.

Innen- / Außenbeurteilungsmodul

Es ist eigentlich leicht, drinnen und draußen zu beurteilen. Sie können die Anzahl der Schnittpunkte zwischen der Halblinie ab dem Zielpunkt und den Umrissen aller Teile überprüfen, aus denen das Donut-förmige Polygon besteht. Wenn die Anzahl der Kreuzungen ungerade ist, wird sie innen beurteilt, und wenn sie gerade ist, wird sie außerhalb __ beurteilt. Der Quellcode sieht wie folgt aus. Im folgenden Beispiel wird die Anzahl der Schnittpunkte zwischen der halben Linie und dem vom Zielpunkt nach Osten gezeichneten Umriss gezählt.

import sys
import numpy  as np
import shapefile
import matplotlib.pyplot as plt
#
def naig(long,lat,pnts,boxs):
　　#Innen / Außen Urteil:Halb geradeaus nach Osten vom Zielpunkt(=A)Zählt die Anzahl der Schnittpunkte der geraden Linie und den Umriss des Features, wenn
    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]:   #Verwenden Sie ein Rechteck, um die Teile einzugrenzen, die sich schneiden können
            NP=len(pnt)
            for np in range(NP-1):
                xx1=pnt[np][0]  #Liniensegmentstartpunkt X.
                yy1=pnt[np][1]  #Liniensegmentstartpunkt Y.
                xx2=pnt[np+1][0]  #Linienendpunkt X.
                yy2=pnt[np+1][1]  #Endpunkt Y des Liniensegments
                miny=min([yy1,yy2]) #Minimum Y des Liniensegments
                maxy=max([yy1,yy2]) #Maximales Y des Liniensegments
                minx=min([xx1,xx2]) #Minimum X des Liniensegments
                maxx=max([xx1,xx2]) #Maximales X für Liniensegment
                if abs(xx1-xx2)<1.E-7: #Das Liniensegment ist vertikal in Nord-Süd-Richtung
                    if xx1>long and miny<=lat<maxy: #Die Y-Koordinate des Zielpunkts ist größer oder gleich dem Minimum Y des Liniensegments und kleiner als das Maximum Y.
                        ccount=ccount+1
                elif abs(yy1-yy2)<1.E-7: #Linien sind fast horizontal
                    ccount=ccount        #Für horizontale Linien ist keine Schnittbewertung erforderlich
                else:
                    aa=(yy2-yy1)/(xx2-xx1) #Die Steigung der geraden Linie, die durch das Liniensegment verläuft
                    bb=yy2-aa*xx2          #Ein Abschnitt einer geraden Linie, der durch ein Liniensegment verläuft
                    yc=lat                 #Y-Schnittpunktkoordinate von Gerader und A.
                    xc=(yc-bb)/aa          #X-Koordinate des Schnittpunkts von Gerader und A.
                    #Die Schnittpunkt-X-Koordinate muss größer als die Zielpunkt-X-Koordinate sein, und die Schnittpunkt-Y-Koordinate muss größer oder gleich dem minimalen Y und kleiner als das maximale Y des Liniensegments sein.
                    if xc>long and miny<=yc<maxy:   
                        ccount=ccount+1
    return ccount
#
fname="test.shp"    #Eingabedatei
src=shapefile.Reader(fname,encoding='SHIFT-JIS')  #Datei lesen
SRS=src.shapes()       #Funktionsinformationen abrufen
SRR=src.shapeRecords() #Funktionsinformationen abrufen
#
#Funktionsinformationen abrufen
pntall=[]
for srs in SRS:
    IP=srs.parts    #Teile für verschiedene Teile erwerben
    NP=len(IP)      #Holen Sie sich die Anzahl der Teile
    pnt=srs.points  #Knotendaten abrufen
    pnts=[]         #Organisieren Sie die Knotenkoordinaten der Teile für jedes Teil
    for ip in range(NP-1):
        pnts.append(pnt[IP[ip]:IP[ip+1]])
    pnts.append(pnt[IP[NP-1]:])
    pntall.append(pnts)
#
#Attributinformationen abrufen
recall=[]
for srr in SRR:
    recall.append(srr.record)
if len(pntall)!=len(recall):
    print("Nichtübereinstimmung zwischen der Anzahl der Attributinformationen und der Anzahl der Features!!")
    sys.exit()
#
#Suchen Sie das Rechteck, das jedes Teil umschließt(Informationen zu Rechtecken abrufen)
NPA=len(pntall)     #Anzahl der Funktionen
boxall=[]
for npa in range(NPA):
    pnts=pntall[npa]    
    NPS=len(pnts)   #Anzahl der Teile
    boxs=[]
    for nps in range(NPS):  #Ermitteln Sie die minimalen und maximalen Werte für den Breiten- und Längengrad der Knoten jedes Teils
        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)
#
#Zählen Sie die Anzahl der Kreuzungen
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)   #Anzahl der Punkte
for npp in range(NPP):      #Schleife über den Punkt
    for npa in range(NPA):  #Schleife über Funktionen
        ic=naig(LON[npp],LAT[npp],pntall[npa],boxall[npa])
        if ic % 2==0:  #Wenn außen
            print("Punkt{0}Ist eine Funktion{1}Draußen! (Anzahl der Kreuzungen:{2})".format(npp+1,npa+1,ic))
        else:       #Wenn drinnen
            print("Punkt{0}Ist eine Funktion{1}im! (Anzahl der Kreuzungen:{2})".format(npp+1,npa+1,ic))

Ausführungsergebnis Ich konnte wie folgt beurteilen. Es scheint, dass dieses selbst erstellte Modul viel schneller ist als die Methode, bei der Sympy für die interne / externe Beurteilung verwendet wird.

Punkt 1 liegt außerhalb von Merkmal 1! (Anzahl der Kreuzungen:14)
Punkt 2 liegt außerhalb von Merkmal 1! (Anzahl der Kreuzungen:18)
Punkt 3 ist in Merkmal 1! (Anzahl der Kreuzungen:3)
Punkt 4 liegt außerhalb von Merkmal 1! (Anzahl der Kreuzungen:4)
Punkt 5 ist in Merkmal 1! (Anzahl der Kreuzungen:1)