Inscribed_rectangle_2

import geopandas as gpd
import numpy as np
import cv2 as cv
from shapely.geometry import Polygon

def largest_interior_rectangle(polygon):
    max_area = 0
    max_rect = None

    for i in range(len(polygon)):
        for j in range(i, len(polygon)):
            for k in range(j, len(polygon)):
                for l in range(k, len(polygon)):
                    if is_convex_quad(polygon[i], polygon[j], polygon[k], polygon[l]):
                        area = (polygon[k][0] - polygon[i][0]) * (polygon[l][1] - polygon[j][1])
                        if area > max_area:
                            max_area = area
                            max_rect = [polygon[i], polygon[j], polygon[k], polygon[l]]

    return max_rect


def is_convex_quad(p1, p2, p3, p4):
    # Compute the convex hull of the quadrilateral
    hull = cv.convexHull(np.array([p1, p2, p3, p4], dtype=np.float32))

    # Check if the midpoint of the diagonal lies inside the convex hull
    return cv.pointPolygonTest(hull, ((p1[0] + p3[0]) / 2, (p1[1] + p3[1]) / 2), False) >= 0


# Read the GeoPackage file
data = gpd.read_file(
     r"C:\Working_desktop\GATE_stuff\Location_allocation\New_data_07022024\Test_processes\test_set_singlepart.gpkg")

# Initialize an empty list to store the largest interior rectangles for each polygon
largest_rectangles = []

# Iterate over each polygon in the GeoDataFrame
for index, row in data.iterrows():
    # Extract the polygon coordinates
    polygon = np.array(row.geometry.exterior.coords)

    # Calculate the largest interior rectangle for the polygon
    rectangle = largest_interior_rectangle(polygon)

    # Append the largest rectangle to the list
    largest_rectangles.append(rectangle)

    # Visualize the polygon and its largest interior rectangle
    img = np.zeros((160, 240, 3), dtype="uint8")
    cv.polylines(img, [polygon.astype(int)], True, (0, 0, 255), 1)
    cv.rectangle(img, (int(rectangle[0][0]), int(rectangle[0][1])), (int(rectangle[2][0]), int(rectangle[2][1])),
                 (255, 0, 0), 1)

    cv.imshow('Largest Interior Rectangle', img)
    cv.waitKey(0)
    cv.destroyAllWindows()

# Print the largest rectangles for each polygon
print("Largest Rectangles:")
for idx, rectangle in enumerate(largest_rectangles):
    print(f"Polygon {idx + 1}: {rectangle}")

# Export the largest inscribed rectangles as a GeoPackage
# Create a list to store the largest rectangles as Polygon objects
rect_polygons = []

# Iterate over the largest rectangles and convert them to Shapely Polygon objects
for idx, rectangle in enumerate(largest_rectangles):
    rect_coords = [(point[0], point[1]) for point in rectangle]
    rect_poly = Polygon(rect_coords)
    rect_polygons.append(rect_poly)

# Create a GeoDataFrame from the list of Polygon objects
gdf_rectangles = gpd.GeoDataFrame(geometry=rect_polygons, crs=data.crs)

# Define the output GeoPackage file path
output_gpkg = (r"C:\Working_desktop\GATE_stuff\Location_allocation\New_data_07022024\Test_processes\rect.gpkg")
# Save the GeoDataFrame to a GeoPackage file
gdf_rectangles.to_file(output_gpkg, driver="GPKG")