map_extend geopandas

mapextend_test.py

import numpy as np
from shapely.ops import unary_union, voronoi_diagram
from shapely.geometry import Polygon, MultiPolygon
import geopandas as gpd
import pandas as pd

def read_data(file_path):
    return gpd.read_file(file_path)

def create_outline(df):
    df_dis = df.dissolve().boundary
    return gpd.GeoDataFrame(geometry=gpd.GeoSeries(df_dis), crs=df.crs)

def create_union(df_dis, df):
    return gpd.overlay(df_dis, df, how='union').reset_index(drop=True)

def generate_points(outline, delta_d=100):
    distances = np.arange(0, outline.length, delta_d)
    points = [outline.interpolate(distance) for distance in distances]
    unique_points = list(set(points))
    return unary_union(unique_points)

def create_buffer(df, buffer_distance=15000):
    buffer_df = df.dissolve().buffer(buffer_distance)
    return gpd.GeoDataFrame(geometry=gpd.GeoSeries(buffer_df), crs=df.crs)

def create_voronoi_diagram(mp, buffer_df):
    voronoi_polygons = voronoi_diagram(mp, envelope=buffer_df.unary_union)
    vor = gpd.GeoDataFrame(geometry=gpd.GeoSeries(voronoi_polygons), crs=buffer_df.crs).explode(index_parts=False)
    vor = vor.set_geometry("geometry")
    vor = vor[vor.geometry.apply(lambda x: isinstance(x, Polygon) or isinstance(x, MultiPolygon))]
    return gpd.clip(vor, buffer_df.unary_union, keep_geom_type=True)

def perform_spatial_join(vor, result):
    joined = gpd.sjoin(vor, result, how='left')
    col = list(joined.columns)
    col.remove("geometry")
    joined = joined.dissolve(col).reset_index()
    return joined

def remove_overlaps(joined, df):
    joined['geometry'] = joined['geometry'].apply(lambda x: unary_union([x]))
    joined = gpd.overlay(joined, df, how='difference')
    joined['geometry'] = joined['geometry'].apply(lambda x: unary_union([x]))
    return joined

def process_polygons(joined):
    joined = joined.explode(index_parts=False)
    result_gdf = gpd.GeoDataFrame(columns=joined.columns, crs=joined.crs)

    for idx, polygon in joined.iterrows():
        if result_gdf.empty:
            result_gdf = gpd.GeoDataFrame([polygon], columns=joined.columns, crs=joined.crs)
        else:
            difference = gpd.overlay(gpd.GeoDataFrame([polygon], columns=joined.columns, crs=joined.crs), result_gdf, how='difference')
            result_gdf = gpd.GeoDataFrame(pd.concat([result_gdf, difference], ignore_index=True), columns=joined.columns, crs=joined.crs)
    
    return result_gdf

def concatenate_dataframes(result_gdf, df):
    extended_df = pd.concat([result_gdf, df])
    extended_df = extended_df.explode(index_parts=False).reset_index(drop=True)
    extended_df['geometry'] = extended_df['geometry'].apply(lambda x: unary_union([x]) if isinstance(x, (Polygon, MultiPolygon)) else x)
    return extended_df

def save_to_file(df, file_path):
    df.to_file(file_path, driver='GPKG')

def main():
    file_path = 'somedata.gpkg'
    df = read_data(file_path)

    df_dis = create_outline(df)
    result = create_union(df_dis, df)

    outline = df_dis.loc[0, 'geometry']
    mp = generate_points(outline, delta_d=100)

    buffer_df = create_buffer(df, buffer_distance=15000)
    vor = create_voronoi_diagram(mp, buffer_df)

    joined = perform_spatial_join(vor, result)
    joined = remove_overlaps(joined, df)

    result_gdf = process_polygons(joined)
    extended_df = concatenate_dataframes(result_gdf, df)

    save_to_file(extended_df, 'extended_data.gpkg')

if __name__ == "__main__":
    main()

a little bit not working, got some overlapping polygons, need to sort out why it does that.

-edited 2024/07/09

overlapping polygons solved, but is it satisfactory?