Complex land parcels often consist of various polygon types, including simple polygons, polygons with interior voids, and composite structures comprising multiple complex polygons. This approach proecsses shapefile data to compute the largest inscribed circle for each polygon using Voronoi diagram methods.

Required Libraries

import numpy as np
from shapely.geometry import Polygon, MultiPolygon, Point
from scipy.spatial import Voronoi
import geopandas as gpd
import matplotlib.pyplot as plt
from pyproj import Transformer

Adaptive Boundary Sampling Strategy

def generate_edge_samples(geometry, sample_count=300):
    if not geometry.is_valid:
        geometry = geometry.buffer(0)
    
    sample_collection = []
    
    if isinstance(geometry, Polygon):
        boundaries = [geometry.exterior] + list(geometry.interiors)
        adjusted_count = sample_count if len(geometry.interiors) > 0 else int(sample_count * 0.75)
        
        total_boundary_length = sum([boundary.length for boundary in boundaries])
        
        for boundary in boundaries:
            segment_length = boundary.length
            segment_samples = int(adjusted_count * (segment_length / total_boundary_length))
            
            sample_distances = np.linspace(0, segment_length, segment_samples)
            sample_collection += [boundary.interpolate(dist) for dist in sample_distances]
    
    return sample_collection

Maximum Inscribed Circle Calcullation

def compute_maximum_inscribed_circle(polygon, sample_density=300):
    boundary_samples = generate_edge_samples(polygon, sample_density)
    unique_points = np.unique(np.array([[pt.x, pt.y] for pt in boundary_samples]), axis=0)
    
    if len(unique_points) < 3:
        raise ValueError("Insufficient unique points for Voronoi diagram construction.")
    
    voronoi_diagram = Voronoi(unique_points)
    
    valid_vertices = [Point(vertex) for vertex in voronoi_diagram.vertices 
                      if polygon.contains(Point(vertex))]
    
    max_radius = 0
    optimal_center = None
    
    for vertex in valid_vertices:
        boundary_distance = polygon.boundary.distance(vertex)
        if boundary_distance > max_radius:
            max_radius = boundary_distance
            optimal_center = vertex
    
    return optimal_center, max_radius

Coordinate Transformation Function

def convert_coordinates(x_coord, y_coord):
    coordinate_converter = Transformer.from_crs("epsg:4545", "epsg:4326", always_xy=True)
    transformed_x, transformed_y = coordinate_converter.transform(x_coord, y_coord)
    return transformed_x, transformed_y

Shapefile Processing Implementation

def process_shapefile_data(input_path, output_path):
    spatial_data = gpd.read_file(input_path)
    
    composite_count = 0
    composite_ids = []
    circle_centers = []
    circle_radii = []
    
    for idx, geometry in enumerate(spatial_data.geometry):
        if isinstance(geometry, MultiPolygon):
            composite_count += 1
            composite_ids.append(idx)
            
            current_max_radius = 0
            current_center = None
            
            for poly in geometry.geoms:
                sample_size = 300 if len(poly.interiors) > 0 else int(300 * 0.75)
                center, radius = compute_maximum_inscribed_circle(poly, sample_size)
                
                if radius > current_max_radius:
                    current_max_radius = radius
                    current_center = center
            
            if current_center:
                transformed_x, transformed_y = convert_coordinates(current_center.x, current_center.y)
                current_center = Point(transformed_x, transformed_y)
            
            circle_centers.append(current_center)
            circle_radii.append(current_max_radius)
            
        elif isinstance(geometry, Polygon):
            sample_size = 300 if len(geometry.interiors) > 0 else int(300 * 0.75)
            center, radius = compute_maximum_inscribed_circle(geometry, sample_size)
            
            if center:
                transformed_x, transformed_y = convert_coordinates(center.x, center.y)
                center = Point(transformed_x, transformed_y)
            
            circle_centers.append(center)
            circle_radii.append(radius)
        else:
            circle_centers.append(None)
            circle_radii.append(0)
    
    spatial_data['center_x'] = [center.x if center else None for center in circle_centers]
    spatial_data['center_y'] = [center.y if center else None for center in circle_centers]
    spatial_data['radius'] = circle_radii
    
    spatial_data.to_file(output_path, encoding='UTF-8')
    return composite_count, composite_ids

Visualization Functions

def visualize_polygon_with_voids(geometry, axis):
    exterior_x, exterior_y = geometry.exterior.xy
    axis.plot(exterior_x, exterior_y, 'b')
    
    for interior in geometry.interiors:
        interior_x, interior_y = interior.xy
        axis.plot(interior_x, interior_y, 'b', linestyle='--')

def plot_inscribed_circle(axis, center_point, circle_radius):
    circle_obj = plt.Circle((center_point.x, center_point.y), circle_radius, 
                           color='r', fill=False)
    axis.add_artist(circle_obj)
    axis.plot(center_point.x, center_point.y, 'ro')

def visualize_processing_results(spatial_data):
    for idx, geometry in enumerate(spatial_data.geometry[29:40]):
        figure, axis = plt.subplots()
        
        if isinstance(geometry, MultiPolygon):
            for poly in geometry.geoms:
                visualize_polygon_with_voids(poly, axis)
                center, radius = compute_maximum_inscribed_circle(poly)
                plot_inscribed_circle(axis, center, radius)
        elif isinstance(geometry, Polygon):
            visualize_polygon_with_voids(geometry, axis)
            center, radius = compute_maximum_inscribed_circle(geometry)
            plot_inscribed_circle(axis, center, radius)
        
        axis.set_aspect('equal')
        plt.title(f'Polygon {idx + 10} with Maximum Inscribed Circle')
        plt.show()

Implemantation Example

input_file = 'shp/cut.shp'
output_file = 'new_shp/cut.shp'

spatial_data = gpd.read_file(input_file)
composite_count, composite_ids = process_shapefile_data(input_file, output_file)
visualize_processing_results(spatial_data)

print("Composite polygon count:", composite_count)
print("Composite polygon IDs:", composite_ids)