Merge pull request #8 from JonasCinquini/main

Merging main into dev

Author: eciraci

PtsLine.py

@@ -12,6 +12,9 @@
 class PtsLine:
     """
     Class to compute the equation of a line passing through two points.
+    Note: this class does not check if:
+    - the two points are the same
+    - the two points are aligned [Vertical or Horizontal Line]
     """
     def __init__(self, x_pt1: float, y_pt1: float,
                  x_pt2: float, y_pt2: float) -> None:

generate_grid.py

@@ -21,6 +21,7 @@
 from shapely.geometry import Polygon
 from shapely.affinity import rotate
 from math import ceil
+from typing import Tuple
 
 
 def add_frame_code_field(grid_gdf):
@@ -123,10 +124,13 @@ def get_fishnet_grid(xmin: float, ymin: float, xmax: float, ymax: float,
     return gdf
 
 
-def rotate_polygon_to_north_up(geometry):
+def rotate_polygon_to_north_up(geometry: Polygon) -> Tuple[Polygon, float]:
     """
     Rotate a polygon to align its orientation with the north.
 
+    NOTE: The function assumes that the polygon has approximately
+        a rectangular shape. The rotation angle is calculated
+        based on the orientation of the shorted side of the polygon.
     Parameters:
         geometry (Polygon): Input polygon.
 
@@ -135,9 +139,27 @@ def rotate_polygon_to_north_up(geometry):
         angle (float): Angle of rotation.
     """
     exterior_coords_list = geometry.exterior.coords[:-1]
-    p1 = min(exterior_coords_list, key=lambda t: t[0])
-    ind_p1 = exterior_coords_list.index(p1)
-    p2 = exterior_coords_list[ind_p1+1]
+    # - Find the lowest corner in the north-south direction
+    p_south = min(exterior_coords_list, key=lambda t: t[1])
+
+    # - Find the furthest corners in the x-direction
+    # - to approximate the orientation of the polygon
+    p_east = max(exterior_coords_list, key=lambda t: t[0])
+    ind_p_east = exterior_coords_list.index(p_east)
+    p_west = min(exterior_coords_list, key=lambda t: t[0])
+    ind_p_west = exterior_coords_list.index(p_west)
+
+    # - Find the corner closest to the south corner
+    # - to approximate the orientation of the polygon
+    dist_east = math.sqrt((p_south[0] - p_east[0])**2
+                          + (p_south[1] - p_east[1])**2)
+    dist_west = math.sqrt((p_south[0] - p_west[0])**2
+                          + (p_south[1] - p_west[1])**2)
+    p1 = p_south
+    if dist_east < dist_west:
+        p2 = exterior_coords_list[ind_p_east]
+    else:
+        p2 = exterior_coords_list[ind_p_west]
 
     angle = math.degrees(math.atan2(p2[1] - p1[1], p2[0] - p1[0]))
     centroid = (geometry.centroid.x, geometry.centroid.y)

iride_csk_frame_grid_utils.py

@@ -0,0 +1,37 @@
+# -  Python Dependencies
+from __future__ import print_function
+import os
+import argparse
+from datetime import datetime
+import numpy as np
+import geopandas as gpd
+from rm_z_coord import rm_z_coord
+
+
+def main() -> None:
+    # - Path to data
+    dat_path = os.path.join(os.path.expanduser("~"), 'Desktop', 'MapItaly',
+                            'MAPITALY_CSG1_2_CSK1_2.shp')
+
+    # - Read data
+    gdf = gpd.read_file(dat_path)
+    print(f"# - Data loaded: {dat_path}")
+    print(f"# - Data shape: {gdf.shape}")
+    print(f"# - Data columns: {gdf.columns}")
+
+    # - Remove Z-Coordinate from geometry
+    gdf = rm_z_coord(gdf)
+
+    # - Loop through the GeoDataFrame lines and extract a reference grid
+    # - for each sub-track.
+    for _, row in gdf.iterrows():
+        # - Check Polygon Validity
+        print(row['geometry'].is_valid)
+
+
+# - run main program
+if __name__ == '__main__':
+    start_time = datetime.now()
+    main()
+    end_time = datetime.now()
+    print(f"# - Computation Time: {end_time - start_time}")

mapitaly_at_grid.py

@@ -0,0 +1,34 @@
+#!/usr/bin/env python
+"""
+Written by Enrico Ciraci'
+January 2024
+
+Reproject a Shapely geometry.
+"""
+from shapely.geometry import Polygon, Point
+from pyproj import CRS, Transformer
+from shapely.ops import transform
+
+
+def reproject_geometry(geometry: Polygon | Point,
+                       source_epsg: int, target_epsg: int) -> Polygon | Point:
+    """
+    Reproject a Shapely geometry.
+    Args:
+        geometry: shapely.geometry.Polygon | shapely.geometry.Point
+        source_epsg: source EPSG code
+        target_epsg: target EPSG code
+
+    Returns: shapely.geometry.Polygon | shapely.geometry.Point
+    """
+    # Set up coordinate reference systems
+    target_crs = CRS.from_epsg(target_epsg)
+    # Create a transformer for the coordinate transformation
+    transformer = Transformer.from_crs(source_epsg, target_crs, always_xy=True)
+    try:
+        # - Polygon
+        return geometry.apply(lambda geom:
+                              transform(transformer.transform, geom))
+    except AttributeError:
+        # - Point
+        return transform(transformer.transform, geometry)

reproject_geometry.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python
+import geopandas as gpd
+from shapely.geometry import Polygon
+"""
+Written by Enrico Ciraci'
+January 2024
+
+Remove z coordinate from a GeoDataFrame.
+Convert a GeoDataFrame with z coordinate to a GeoDataFrame without z coordinate.
+"""
+
+
+def rm_z_coord(gdf) -> gpd.GeoDataFrame:
+    """
+    Remove z coordinate from a GeoDataFrame
+    Args:
+        gdf: GeoDataFrame
+
+    Returns:
+
+    """
+    # - remove z coordinate
+    no_z_coord = []
+    for _, row in gdf.iterrows():
+        # - remove z coordinate
+        z_coords = row['geometry'].exterior.coords[:-1]
+        no_z_coord.append(Polygon([(crd[0], crd[1]) for crd in z_coords]))
+    gdf['geometry'] = no_z_coord
+
+    return gdf

rm_z_coord.py

@@ -2,7 +2,7 @@
 """ Unit tests for the reproject_geometry function. """
 import pytest
 from shapely.geometry import Point, Polygon
-from generate_grid import reproject_geometry
+from reproject_geometry import reproject_geometry
 
 @pytest.fixture
 def example_polygon():

test/reproj_test.py

@@ -0,0 +1,22 @@
+import geopandas as gpd
+from shapely.geometry import Polygon
+from generate_grid import generate_grid
+
+
+def test_generate_grid():
+    # Create a sample GeoDataFrame for testing
+    d = {'index': 1,
+         'geometry': Polygon([(12, 36), (11.5, 38), (14, 38.5),
+                              (14, 36.5), (12, 36)])}
+    gdf_t = gpd.GeoDataFrame(d, crs='EPSG:4326', index=[0])
+
+    n_c = 3
+    az_res = 5000
+    buffer_dist = 1000
+
+    # Ensure the function runs without errors
+    result_gdf = generate_grid(gdf_t, n_c, az_res, buffer_dist)
+
+    # Perform assertions based on your expectations for the result
+    assert isinstance(result_gdf, gpd.GeoDataFrame)
+    assert len(result_gdf) > 0

test/test_generate_grid.py

@@ -0,0 +1,20 @@
+import geopandas as gpd
+from shapely.geometry import Polygon
+from iride_csk_frame_grid_utils import rotate_polygon_to_north_up
+from math import isclose
+
+
+def test_rotate_polygon_to_north_up():
+    """
+    Test the rotate_polygon_to_north_up function.
+    """
+    # Create a sample GeoDataFrame for testing
+    polygon = Polygon([(12, 36), (11.5, 38), (14, 38.5),
+                       (14, 36.5), (12, 36)])
+
+    # Ensure the function runs without errors
+    rotated_polygon, angle = rotate_polygon_to_north_up(polygon)
+
+    # Perform assertions based on your expectations for the result
+    assert isinstance(rotated_polygon, Polygon)
+    assert isinstance(angle, float)