centerline.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from shapely.geometry import LineString
from shapely.geometry import MultiLineString
from scipy.spatial import Voronoi
import numpy as np

# from python geospatial analysis cookbook
# https://books.google.com.br/books?id=rPCoCwAAQBAJ&lpg=PA184&ots=1rzuvq36UZ&dq=class%20Centerline(object)%3A%20%20%20%20%20def%20__init__(self%2C%20inputGEOM%2C%20dist%3D0.5)%3A%20%20%20%20%20%20%20%20%20self.inputGEOM%20%3D%20inputGEOM%20%20%20%20%20%20%20%20%20self.dist%20%3D%20abs(dist)&hl=pt-BR&pg=PA184#v=onepage&q&f=false

class Centerline(object):
    def __init__(self, inputGEOM, dist=0.5):
        self.inputGEOM = inputGEOM
        self.dist = abs(dist)

    def create_centerline(self):
        """
        Calculates the centerline of a polygon.
        Densifies the border of a polygon which is then represented
        by a Numpy array of points necessary for creating the
        Voronoi diagram. Once the diagram is created, the ridges
        located within the polygon are joined and returned.
        Returns:
            a MultiLinestring located within the polygon.
        """

        minx = int(min(self.inputGEOM.envelope.exterior.xy[0]))
        miny = int(min(self.inputGEOM.envelope.exterior.xy[1]))

        border = np.array(self.densify_border(self.inputGEOM, minx, miny))

        vor = Voronoi(border)
        vertex = vor.vertices

        lst_lines = []
        for j, ridge in enumerate(vor.ridge_vertices):
            if -1 not in ridge:
                line = LineString([
                    (vertex[ridge[0]][0] + minx, vertex[ridge[0]][1] + miny),
                    (vertex[ridge[1]][0] + minx, vertex[ridge[1]][1] + miny)])

                if line.within(self.inputGEOM) and len(line.coords[0]) > 1:
                    lst_lines.append(line)

        return MultiLineString(lst_lines)

    def densify_border(self, polygon, minx, miny):
        """
        Densifies the border of a polygon by a given factor
        (by default: 0.5).
        The function tests the complexity of the polygons
        geometry, i.e. does the polygon have holes or not.
        If the polygon doesn't have any holes, its exterior
        is extracted and densified by a given factor. If the
        polygon has holes, the boundary of each hole as
        well as its exterior is extracted and densified
        by a given factor.
        Returns:
            a list of points where each point is represented
            by a list of its
            reduced coordinates.
        Example:
            [[X1, Y1], [X2, Y2], ..., [Xn, Yn]
        """

        if len(polygon.interiors) == 0:
            exterior_line = LineString(polygon.exterior)
            points = self.fixed_interpolation(exterior_line, minx, miny)

        else:
            exterior_line = LineString(polygon.exterior)
            points = self.fixed_interpolation(exterior_line, minx, miny)

            for j in range(len(polygon.interiors)):
                interior_line = LineString(polygon.interiors[j])
                points += self.fixed_interpolation(interior_line, minx, miny)

        return points

    def fixed_interpolation(self, line, minx, miny):
        """
        A helping function which is used in densifying
        the border of a polygon.
        It places points on the border at the specified distance.
        By default the distance is 0.5 (meters) which means
        that the first point will be placed 0.5 m from the
        starting point, the second point will be placed at the
        distance of 1.0 m from the first point, etc. Naturally,
        the loop breaks when the summarized distance exceeds
        the length of the line.
        Returns:
            a list of points where each point is represented by
            a list of its reduced coordinates.
        Example:
            [[X1, Y1], [X2, Y2], ..., [Xn, Yn]
        """

        count = self.dist
        newline = []

        startpoint = [line.xy[0][0] - minx, line.xy[1][0] - miny]
        endpoint = [line.xy[0][-1] - minx, line.xy[1][-1] - miny]
        newline.append(startpoint)

        while count < line.length:
            point = line.interpolate(count)
            newline.append([point.x - minx, point.y - miny])
            count += self.dist

        newline.append(endpoint)

        return newline

centerline_test1.py

from osgeo import ogr,osr
from shapely.geometry import shape
import centerline
import json

hidro = r"/home/c/Documents/data/massa_dagua.shp"
out_shp = r"/home/c/Documents/data/centerlinev1.shp"
centerline_distance = 1.0

EPSG_CODE = 4326
SIMPLIFY_TOLERANCE = 1.0

ds = ogr.Open(hidro)
layer = ds.GetLayer()
driver = ogr.GetDriverByName('Esri Shapefile')
out_ds = driver.CreateDataSource(out_shp)
srs = osr.SpatialReference()
srs.ImportFromEPSG(EPSG_CODE)
out_layer = out_ds.CreateLayer('', srs, ogr.wkbMultiLineString)
# Add one attribute
out_layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
## If there are multiple geometries, put the "for" loop here
defn = out_layer.GetLayerDefn()
n = 0
for feature in layer:
    print('starting feature {}'.format(n))
    feat_geom = feature.geometry()
    simple = feat_geom.Simplify(SIMPLIFY_TOLERANCE)
    geojson = simple.ExportToJson()
    geojson = json.loads(geojson)
    if len(geojson['coordinates']) == 0:
        continue
    shp_geom = shape(geojson)

    
    try:
        cent = centerline.Centerline(shp_geom, dist=centerline_distance)
        line = cent.create_centerline()

        # Create a new feature (attribute and geometry)
        out_feat = ogr.Feature(defn)
        out_feat.SetField('id', n)
        # Make a geometry, from Shapely object
        geom = ogr.CreateGeometryFromWkb(line.wkb)
        out_feat.SetGeometry(geom)
        out_layer.CreateFeature(out_feat)
        out_feat = None
        print('finished!')
        n+=1
        
    except Exception as e:
        print(e)
    
out_ds = None