import pandas as pd
import geopandas as gpd
import numpy as np
import shapely.geometry as geo

hexagon_radius = 1000
operating_areas_path = "/run/media/sebastian/shoerl_data/astra1802/gis/operating_areas.shp"
output_path = "/run/media/sebastian/shoerl_data/astra1802/gis/waiting_time_areas.shp"

df_operating_areas = gpd.read_file(operating_areas_path)

corners = [
    np.array([np.cos(30.0 * np.pi / 180.0), np.sin(30.0 * np.pi / 180.0)]) * hexagon_radius,
    np.array([np.cos(90.0 * np.pi / 180.0), np.sin(90.0 * np.pi / 180.0)]) * hexagon_radius,
    np.array([np.cos(150.0 * np.pi / 180.0), np.sin(150.0 * np.pi / 180.0)]) * hexagon_radius,
    np.array([np.cos(210.0 * np.pi / 180.0), np.sin(210.0 * np.pi / 180.0)]) * hexagon_radius,
    np.array([np.cos(270.0 * np.pi / 180.0), np.sin(270.0 * np.pi / 180.0)]) * hexagon_radius,
    np.array([np.cos(330.0 * np.pi / 180.0), np.sin(330.0 * np.pi / 180.0)]) * hexagon_radius
]

offset_x = 2.0 * hexagon_radius * np.cos(30.0 * np.pi / 180.0)

t = 2.0 * hexagon_radius * np.cos(30.0 * np.pi / 180.0) / np.sqrt(3.0)
offset_y = hexagon_radius + 0.5 * t

geometry = []

for area_name, area_geometry in zip(df_operating_areas["OAREA"], df_operating_areas["geometry"]):
    min_x, min_y, max_x, max_y = area_geometry.bounds

    max_i = int(np.ceil((max_x - min_x) / offset_x)) + 1
    max_j = int(np.ceil((max_y - min_y) / offset_y)) + 1

    index = 0

    for i in range(max_i):
        for j in range(max_j):
            centroid = np.array([
                min_x + offset_x * i + (0.5 * offset_x if j % 2 == 0 else 0.0),
                min_y + offset_y * j
            ])

            hexagon = geo.Polygon(corners + centroid)

            if hexagon.intersects(area_geometry):
                geometry.append((area_name, index, hexagon))
                index += 1

df = pd.DataFrame.from_records(geometry, columns = ["OAREA", "IDX", "geometry"])
df = gpd.GeoDataFrame(df, crs = {"init": "EPSG:2056"})
df.to_file(output_path)