Geometry 常用工具

Geometry 常用工具

• 常用库

  • shapely

  • GeoPandas

  • scikit-geometry

• Voronoi Diagram

• Delaunay Triangulation

1. Polygon Skeleton

原理：Li et al. (2014) Section 2.1

• Constrained Delaunay triangulation 将 Polygon 转为三角面

2. 生成图形

2.1. 创建 Grid

2.2. Voronoi

可实现的方法：

• scipy.spatial.Voronoi

• shapely.ops.voronoi_diagram

实例： GeoDataFrame 生成 Voronoi Diagram

点击查看代码

def generate_voronoi(node_gpd, max_dist=None, boundary=None):
    '''
    creat mesh grids
    
    parameters:
    -----------
    node_gpd (geopandas.GeoDataFrame) : 
        nodes, with CRS information
    max_dist (float, default is None): 
        maximum distance from nodes
    boundary (shapely.geometry.Polygon or shapely.geometry.MultiPolygon, default is None) : 
        a boundary, must share the same CRS with "node_gpd"

    Return:
    -------
    voro (geopandas.GeoDataFrame) :
        
    '''
    from shapely.geometry import MultiPoint, Polygon, MultiPolygon
    from shapely.ops import voronoi_diagram
    
    # voronoi
    points = MultiPoint(node_gpd['geometry'].to_list())
    voro = voronoi_diagram(points)
    # 
    voro = gpd.GeoSeries(list(voro.geoms))
    voro = gpd.GeoDataFrame(list(range(voro.shape[0])), 
                            columns  = ['voro_id'], 
                            geometry = voro.geometry,
                            crs      = node_gpd.crs)
    
    # join voronoi and node
    voro = voro.sjoin(df = node_gpd, 
                      predicate = 'intersects', 
                      how       = 'left',
                      rsuffix   = 'node',
                      lsuffix   = 'voro')
    
    # use original node index
    voro = voro.set_index('index_node').sort_index()
    
    if not (max_dist is None):
        buffer = node_gpd.buffer(max_dist)
        # intersection
        voro['geometry'] = voro.intersection(buffer, align=True)
    
    # clip
    if not (boundary is None):
        assert isinstance(boundary, Polygon) or isinstance(boundary, MultiPolygon)
        voro = voro.clip(boundary)
        
    return voro

2.3. Delaunay triangulation

可实现的方法：

• scipy.spatial.Delaunay

• shapely.ops.triangulate

实例： GeoDataFrame 生成 Voronoi Diagram

点击查看代码

def delaunay_triangulation(node_gpd, node_id, return_type='polygon', boundary=None):
    '''
    creat delaunay triangulation
    
    parameters:
    -----------
    node_gpd (geopandas.GeoDataFrame) : 
        nodes (shapely.geometry.Point), with CRS information
    node_id (str): 
        column name of node id
    return_type (str, default=polygon, {'polygon', 'line_string', 'line_segement'}):
        
    boundary (shapely.geometry.Polygon or shapely.geometry.MultiPolygon, default is None) : 
        a boundary, must share the same CRS with "node_gpd"

    Return:
    -------
    triangle (geopandas.GeoDataFrame) :
        
    '''
    from shapely.geometry import LineString, Point, Polygon, MultiPolygon
    # from shapely.ops import Delaunay
    from scipy.spatial import Delaunay 
    
    # node index and id
    node_id_mapper = {k : v for k, v in enumerate(node_gpd[node_id].to_list())}
    # 
    node_coors = {k : v.coords[0] for k, v in enumerate(node_gpd['geometry'].to_list())}
    # print(node_coors[0])
    
    # coordinates of points
    points = np.concatenate((node_gpd['geometry'].x.values.reshape(-1, 1), 
                             node_gpd['geometry'].y.values.reshape(-1, 1)), axis=1)

    triangle = pd.DataFrame(Delaunay(points).simplices, 
                            columns=['point_1', 'point_2', 'point_3'])
    
    if return_type == 'polygon':

        triangle['geometry'] = triangle[['point_1', 'point_2', 'point_3']].apply(
            lambda x: Polygon([node_coors[i] for i in x.to_list()]),
            axis='columns')
        
        triangle[['point_1', 'point_2', 'point_3']] = triangle[['point_1', 'point_2', 'point_3']].applymap(lambda x : node_id_mapper[x])

    elif return_type == 'line_string':

        triangle['geometry'] = triangle[['point_1', 'point_2', 'point_3', 'point_1']].apply(
            lambda x: LineString([node_coors[i] for i in x.to_list()]),
            axis='columns')

        triangle[['point_1', 'point_2', 'point_3']] = triangle[['point_1', 'point_2', 'point_3']].applymap(lambda x : node_id_mapper[x])

    elif return_type == 'line_segement':

        triangle = pd.DataFrame(
            np.concatenate((triangle[['point_1', 'point_2']].values, 
                            triangle[['point_2', 'point_3']].values, 
                            triangle[['point_1', 'point_3']].values), axis=0),
            columns = ['point_1', 'point_2']
        )
        
        triangle['geometry'] = triangle[['point_1', 'point_2']].apply(
            lambda x: LineString([node_coors[i] for i in x.to_list()]),
            axis='columns')
        
        triangle[['point_1', 'point_2']] = triangle[['point_1', 'point_2']].applymap(lambda x : node_id_mapper[x])
    
    # create gpd.GeoDataFrame
    triangle = gpd.GeoDataFrame(triangle.drop(['geometry'], axis=1), 
                                geometry = triangle['geometry'],
                                crs = node_gpd.crs)
    
    # clip
    if not (boundary is None):
        assert isinstance(boundary, Polygon) or isinstance(boundary, MultiPolygon)
        
        # triangle = triangle.clip(boundary) # quite time-consuming if "boundary" is complex
        
        ix = triangle.centroid.sindex.query(landuse_road, predicate='intersects')
        triangle = triangle.reindex(ix).sort_index().reset_index()
    
    return triangle

2.4. Constrained Delaunay triangulation

3. 处理图形

3.1 Simplify

可实现的方法：

• shapely: object.simplify

  • Ramer–Douglas–Peucker algorithm

  • Visualization

• geopandas.GeoSeries.simplify

• scikit-geometry.simplify

  • CGAL 2D Polyline Simplification

References

Li, Q., Fan, H., Luan, X., Yang, B., & Liu, L. (2014). Polygon-based approach for extracting multilane roads from OpenStreetMap urban road networks. International Journal of Geographical Information Science, 28(11), 2200-2219.

Samuel Peterson, Computing Constrained Delaunay Triangulation, website

Javascript implementation of the Ramer Douglas Peucker Algorithm, website