Inscribed_rectangle_3

import geopandas as gpd
import matplotlib.pyplot as plt
import shapely.geometry as sg
import shapely.affinity as sa
import pandas as pd


def main():
    # Load your polygon data
    data = gpd.read_file(r"C:\Working_desktop\GATE_stuff\Location_allocation\New_data_07022024\Test_processes\test_set_singlepart.gpkg")

    # Create an empty GeoDataFrame to store the inscribed rectangles
    inscribed_rectangles = pd.DataFrame(columns=['geometry'])

    # Iterate over each polygon in the dataset
    for index, row in data.iterrows():
        polygon = row['geometry']
        inscribed_rectangle = find_inscribed_rectangle(polygon)

        # Check if an inscribed rectangle is found for the current polygon
        if inscribed_rectangle:
            # Append the inscribed rectangle to the GeoDataFrame
            inscribed_rectangles = pd.concat([inscribed_rectangles, pd.DataFrame({'geometry': [inscribed_rectangle]})],
                                             ignore_index=True)

    # Convert DataFrame to GeoDataFrame
    inscribed_rectangles = gpd.GeoDataFrame(inscribed_rectangles, geometry='geometry')

    # Check if the GeoDataFrame is not empty before trying to export it
    if not inscribed_rectangles.empty:
        # Export the GeoDataFrame containing inscribed rectangles to a new GeoPackage file
        inscribed_rectangles.to_file("inscribed_rectangles.gpkg", driver="GPKG")

        # Visualize the result
        fig, ax = plt.subplots()
        data.plot(ax=ax, color='gray', edgecolor='black', alpha=0.3)  # Plot all polygons in semi-transparent gray
        inscribed_rectangles.plot(ax=ax, color='blue', edgecolor='black')  # Plot inscribed rectangles in blue
        plt.title("Inscribed Rectangles")
        plt.xlabel("X-Coordinate")
        plt.ylabel("Y-Coordinate")
        plt.show()
    else:
        print("No inscribed rectangles found in the dataset.")


def find_inscribed_rectangle(polygon):
    max_area = 0
    max_rectangle = None

    # Get all the vertices of the polygon
    vertices = list(polygon.exterior.coords)

    # Iterate over all possible pairs of non-adjacent vertices
    num_vertices = len(vertices)
    for i in range(num_vertices):
        for j in range(i + 2, num_vertices):
            # Create a rectangle using the current pair of vertices
            rectangle = sg.Polygon([vertices[i], vertices[j],
                                    vertices[j - 1], vertices[i - 1]])

            # Check if the rectangle is completely inside the polygon
            if polygon.contains(rectangle):
                # Calculate the area of the rectangle
                area = rectangle.area

                # Update the maximum area and corresponding rectangle if necessary
                if area > max_area:
                    max_area = area
                    max_rectangle = rectangle

    return max_rectangle


if __name__ == '__main__':
    main()