Simple script to create an offset polygon around any shapefile.

shapeoffset.py

"""
Simple script to create an offset polygon around any shapefile.
More information:
    https://forum.step.esa.int/t/define-a-buffer-around-shapefile/20966

author: martino.ferrari@c-s.fr
license: GPLv3
"""
import shapefile as shp  # Requires the pyshp package
import matplotlib.pyplot as plt
import numpy as np
import argparse


def __parse_args__():
    """
    Parses command line arguments using argparse utility.
    """

    parser = argparse.ArgumentParser(description='Create an offset polygon around any shapefile.')
    parser.add_argument('path', type=str,
                        help='shapefile path')
    parser.add_argument('offset', type=float,
                        help='angular offset')
    parser.add_argument('--showplots', type=bool, default=False,
                        help='show the computed offset')

    return parser.parse_args()


def compute_offset(points, offset):
    """
    Computes the offset from a list of points using polar coordinate.

    Parameters:
    -----------
     - points: numpy.ndarray containing euclidean coordinates of the geomtery to offset
     - offset: angular offset

    Returns:
    --------
    numpy.ndarray containing the offseted geomtery
    """
    # compute center as mean point of all points
    center = np.mean(points, 1)

    # translate coordinate system to the center
    points_orig = points.T - center
    # convert cartesian coordinate in polar coordinate
    rho = np.sqrt(np.sum(points_orig** 2, 1))
    theta = np.arctan2(points_orig.T[1], points_orig.T[0])

    # setup variable for angular iteration
    N_STEPS = 100
    ANGULAR_RESOLUTION = 2 * np.pi / N_STEPS

    offseted = [[], []]

    # explore polar space
    for alpha in np.linspace(-np.pi, np.pi, N_STEPS):
        # find index of all points around the current angle alpha
        filtered = np.abs(theta - alpha) < ANGULAR_RESOLUTION
        # get filtered distances
        filt_rho = rho[filtered]
        # check if at least one point is in the current angluar step
        if len(filt_rho) > 0:
            filt_theta = theta[filtered]
            # find farther point on the list
            index = np.argmax(filt_rho)
            # compute offseted point
            xoff = (filt_rho[index] + offset) * np.cos(filt_theta[index])
            yoff = (filt_rho[index] + offset) * np.sin(filt_theta[index])
            # added to the offseted point with correct reference system
            offseted[0].append(xoff + center[0])
            offseted[1].append(yoff + center[1])
    # close the polygon
    offseted[0].append(offseted[0][0])
    offseted[1].append(offseted[1][0])
    return np.array(offseted)


def load_points(shapefile_path):
    """
    Loads all the points of a shapefile in memory as a numpy array.

    Parameters:
    -----------
     - shapefile_path: path of the shapefile

    Returns:
    --------
    numpy.ndarray that contains all points
    """
    # open shape files
    sf = shp.Reader(shapefile_path)
    # load the full list of points of all shapes in a singe object
    points = [[], []]
    for shape in sf.shapeRecords():
        for p in shape.shape.points:
            points[0].append(p[0])
            points[1].append(p[1])
    # convert to numpy array for faster math
    return np.array(points)


def save_offset(path, points):
    """
    Saves offseted file to a new shape file named after the original file.

    Parameters:
    -----------
     - path: source file path
     - points: offseted points
    """
    with shp.Writer(path[:-4]+'_offset.shp') as writer:
        writer.field('name', 'C')
        writer.poly([points.T])
        writer.record('polygon1')


def main():
    """
    ShapeOffset main script, take as input the shape file path and the size of
    the offset in meters.
    """
    # parse arguments
    args = __parse_args__()

    # load points of the file
    points = load_points(args.path)
    # compute offset
    offseted = compute_offset(points, args.offset)

    # display offset and show plot if flag is enabled
    if args.showplots:
        plt.scatter(points[0], points[1], marker='.')
        plt.plot(offseted[0], offseted[1])
        plt.show()

    # save file
    save_offset(args.path, offseted)


if __name__ == "__main__":
    main()