Jugement intérieur / extérieur avec Python ➁: À propos des polygones en forme d'anneau

Chose que tu veux faire

Pour les fichiers shp de polygones en forme d'anneau suivants, nous avons organisé comment juger de l'intérieur et de l'extérieur d'un point arbitraire. (Source: https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A31-v2_1.html)

J'expliquerai dans l'ordre suivant. Étape 1: Lire les données de polygones en forme d'anneau Étape 2: module de jugement intérieur / extérieur

Lecture de données de polygones en forme d'anneau

Dans le fichier shp, le polygone en forme d'anneau est construit en superposant plusieurs polygones inégaux (parties). Par conséquent, tout d'abord, il est nécessaire d'extraire les coordonnées des nœuds de chaque pièce.

import numpy
import shapefile
import matplotlib.pyplot as plt
#
def outLS(pnt,fig):    #Tracez les données de ligne qui connectent les données de nœud du polygone
    print(*pnt)
    xs,ys = zip(*pnt)
    plt.plot(xs,ys) 
    return 0
#
fname="test.shp"    #Fichier d'entrée
src=shapefile.Reader(fname,encoding='SHIFT-JIS')  #Lecture de fichiers
SRS=src.shapes()    #Obtenir des informations sur les fonctionnalités
#
for srs in SRS:
    IP=srs.parts    #Acquérir des pièces pour diverses pièces
    NP=len(IP)      #Obtenez le nombre de pièces
    pnt=srs.points  #Obtenir les données des nœuds
    pnts=[]         #Organisez les coordonnées des nœuds des pièces pour chaque pièce
    for ip in range(NP-1):
        pnts.append(pnt[IP[ip]:IP[ip+1]])
    pnts.append(pnt[IP[NP-1]:])
    #
    fig = plt.figure() #Dessinez des informations de contour pour chaque pièce
    for np in range(NP):
        outLS(pnts[np],fig)
    plt.show()
    plt.clf()

__ Résultat de l'exécution__ Comme indiqué ci-dessous, nous avons pu distinguer et acquérir les informations de contour de chaque partie qui compose le polygone en forme d'anneau.

Module de jugement intérieur / extérieur

Il est en fait facile de juger à l'intérieur et à l'extérieur. __ Il suffit de vérifier le nombre d'intersections entre la demi-droite partant du point cible et le contour de toutes les parties qui composent le polygone en forme d'anneau. Si le nombre d'intersections est impair, il sera jugé à l'intérieur, et s'il est pair, il sera jugé à l'extérieur de __. Le code source ressemble à ce qui suit. Dans l'exemple ci-dessous, le nombre d'intersections entre la demi-ligne et le contour tracé vers l'est à partir du point cible est compté.

import sys
import numpy  as np
import shapefile
import matplotlib.pyplot as plt
#
def naig(long,lat,pnts,boxs):
　　#Jugement intérieur / extérieur:À moitié droit vers l'est à partir du point cible(=A)Compte le nombre d'intersections de la ligne droite et du contour de l'entité lorsque
    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]:   #Utilisez un rectangle pour réduire les parties qui peuvent se croiser
            NP=len(pnt)
            for np in range(NP-1):
                xx1=pnt[np][0]  #Point de départ du segment de ligne X
                yy1=pnt[np][1]  #Point de départ du segment de ligne Y
                xx2=pnt[np+1][0]  #Point de fin de ligne X
                yy2=pnt[np+1][1]  #Point final Y du segment de ligne
                miny=min([yy1,yy2]) #Y minimum du segment de ligne
                maxy=max([yy1,yy2]) #Y maximum du segment de ligne
                minx=min([xx1,xx2]) #X minimum du segment de ligne
                maxx=max([xx1,xx2]) #X maximum pour le segment de ligne
                if abs(xx1-xx2)<1.E-7: #Le segment de ligne est vertical dans la direction nord-sud
                    if xx1>long and miny<=lat<maxy: #La coordonnée Y du point cible est supérieure ou égale au minimum Y du segment de ligne et inférieure au maximum Y
                        ccount=ccount+1
                elif abs(yy1-yy2)<1.E-7: #Les lignes sont presque horizontales
                    ccount=ccount        #Aucun jugement d'intersection requis pour les lignes horizontales
                else:
                    aa=(yy2-yy1)/(xx2-xx1) #La pente de la droite passant par le segment de ligne
                    bb=yy2-aa*xx2          #Une section de ligne droite passant par un segment de ligne
                    yc=lat                 #Coordonnée Y de l'intersection de la ligne droite et A
                    xc=(yc-bb)/aa          #Coordonnée X de l'intersection de la ligne droite et A
                    #La coordonnée X de l'intersection doit être supérieure à la coordonnée X du point cible et la coordonnée Y de l'intersection doit être supérieure ou égale au minimum Y et inférieur au maximum Y du segment de ligne.
                    if xc>long and miny<=yc<maxy:   
                        ccount=ccount+1
    return ccount
#
fname="test.shp"    #Fichier d'entrée
src=shapefile.Reader(fname,encoding='SHIFT-JIS')  #Lecture de fichiers
SRS=src.shapes()       #Obtenir des informations sur les fonctionnalités
SRR=src.shapeRecords() #Obtenir des informations sur les fonctionnalités
#
#Obtenir des informations sur les fonctionnalités
pntall=[]
for srs in SRS:
    IP=srs.parts    #Acquérir des pièces pour diverses pièces
    NP=len(IP)      #Obtenez le nombre de pièces
    pnt=srs.points  #Obtenir les données des nœuds
    pnts=[]         #Organisez les coordonnées des nœuds des pièces pour chaque pièce
    for ip in range(NP-1):
        pnts.append(pnt[IP[ip]:IP[ip+1]])
    pnts.append(pnt[IP[NP-1]:])
    pntall.append(pnts)
#
#Obtenir des informations sur les attributs
recall=[]
for srr in SRR:
    recall.append(srr.record)
if len(pntall)!=len(recall):
    print("Inadéquation entre le nombre d'informations attributaires et le nombre d'entités!!")
    sys.exit()
#
#Trouvez le rectangle qui entoure chaque partie(Obtenir des informations sur le rectangle)
NPA=len(pntall)     #Nombre de fonctionnalités
boxall=[]
for npa in range(NPA):
    pnts=pntall[npa]    
    NPS=len(pnts)   #Nombre de pièces
    boxs=[]
    for nps in range(NPS):  #Obtenez les valeurs minimales et maximales de la latitude et de la longitude des nœuds de chaque partie
        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)
#
#Comptez le nombre d'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)   #Nombre de points
for npp in range(NPP):      #Boucle sur le point
    for npa in range(NPA):  #Boucle sur les fonctionnalités
        ic=naig(LON[npp],LAT[npp],pntall[npa],boxall[npa])
        if ic % 2==0:  #Si à l'extérieur
            print("point{0}Est une fonctionnalité{1}À l'extérieur! (Nombre d'intersections:{2})".format(npp+1,npa+1,ic))
        else:       #Si à l'intérieur
            print("point{0}Est une fonctionnalité{1}dans! (Nombre d'intersections:{2})".format(npp+1,npa+1,ic))

__ Résultat de l'exécution__ J'ai pu juger comme suit. Il semble que ce module self-made soit beaucoup plus rapide que la méthode utilisant sympy pour le jugement interne / externe.

Le point 1 est en dehors de l'élément 1! (Nombre d'intersections:14)
Le point 2 est en dehors de l'élément 1! (Nombre d'intersections:18)
Le point 3 est dans la fonction 1! (Nombre d'intersections:3)
Le point 4 est à l'extérieur de l'élément 1! (Nombre d'intersections:4)
Le point 5 est dans la fonction 1! (Nombre d'intersections:1)