end0tknr's kipple - web写経開発

太宰府天満宮の狛犬って、妙にカワイイ

hands-on shapely for python

https://shapely.readthedocs.io/en/stable/manual.html

上記の readthedocs の中から、いくつかを実装

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path        import Path
from matplotlib.patches     import PathPatch
from matplotlib.collections import PatchCollection

# cf.https://shapely.readthedocs.io/en/stable/manual.html
from shapely.geometry       import Polygon
from shapely.geometry       import MultiPolygon


def main():
    polygon =  Polygon(
        shell=[ (10,10),(10,40),(70,40),(70,10),(60,10),
                (60,30),(50,30),(50,10),(10,10) ],
        holes=[ [(20,20),(40,20),(40,30),(20,30),(20,20)] ]
    )
    
    get_properties( polygon )   # 属性
    relations( polygon )        # 位置関係
    modifies( polygon )         # 変形?
    

def relations( polygon ):
    polygons = MultiPolygon(
        [Polygon( ( (-10,10),(-10,20),( 0,20),( 0,10),(-10,10) ) ),
         Polygon( ( (  0,10),(  0,20),(10,20),(10,10),(  0,10) ) ),
         Polygon( ( (  0,10),(  0,20),(20,20),(20,10),(  0,10) ) ),
         Polygon( ( ( 30,10),( 30,20),(40,20),(40,10),( 30,10) )) ]
    )

    relation_methos = ["EQUAL","CONTAIN","I-SECT","OVERLAP",
                       "CROSS","TOUCH","WITHIN","DISTANCE" ]
    print( "\t".join(relation_methos) )

    for polygon_tmp in polygons.geoms:
        relations = []
        relations.append( str( polygon.equals( polygon_tmp ) ))
        relations.append( str( polygon.contains( polygon_tmp ) ))
        relations.append( str( polygon.intersects( polygon_tmp ) ))
        relations.append( str( polygon.overlaps( polygon_tmp ) ))
        relations.append( str( polygon.crosses( polygon_tmp ) ))
        relations.append( str( polygon.touches( polygon_tmp ) ))
        relations.append( str( polygon.within( polygon_tmp ) ))
        relations.append( str( polygon.distance( polygon_tmp ) ))
        print( "\t".join(relations) )
    
def modifies( polygon ):
    polygon_2 =  Polygon(
        shell=[ (0,10),(0,20),(20,20),(20,10),(0,10) ]
    )
    polygon_3 =  Polygon(
        shell=[ (30,10),(30,20),(40,20),(40,10),(30,10) ]
    )
    polygon = polygon.difference(polygon_2)             # 図形の差分
    polygon = polygon.difference(polygon_3)             # 図形の差分
    #polygon = polygon.intersection(polygon_2)          #  〃 の積
    #polygon = polygon.union(polygon_2)                 #  〃 の和
    #polygon = polygon.symmetric_difference(polygon_2)  #  〃 のXOR
    #return

    fig, ax = plt.subplots()
    
    if polygon.type == "MultiPolygon":
        for polygon_tmp in polygon.geoms:
            plot_polygon(ax, polygon_tmp, facecolor='lightblue', edgecolor='red')
    else:
        plot_polygon(ax, polygon, facecolor='lightblue', edgecolor='red')
        
    plt.show()

    
def get_properties( polygon ):
    print( polygon )
    print( polygon.type )
    print( "外周座標:",list(polygon.exterior.coords) )
    for interior in polygon.interiors:
        print( "内周座標:", list(interior.coords) )
    
    print( "面積:", polygon.area )
    print( "BOUNDING BOX:", polygon.bounds )
    print( "BOUNDING BOX:", polygon.bounds )
    print( "EDGE LENGTH:", polygon.length )

    
# Plots a Polygon to pyplot `ax`
# cf. https://stackoverflow.com/questions/55522395
def plot_polygon(ax, poly, **kwargs):
    path = Path.make_compound_path(
        Path(np.asarray(poly.exterior.coords)[:, :2]),
        *[Path(np.asarray(ring.coords)[:, :2]) for ring in poly.interiors])

    patch = PathPatch(path, **kwargs)
    collection = PatchCollection([patch], **kwargs)
    
    ax.add_collection(collection, autolim=True)
    ax.autoscale_view()
    return collection

    
if __name__ == '__main__':
    main()