Merge pull request #8 from JonasCinquini/main

Merging main into dev

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:

iride_csk_frame_grid_utils.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)

mapitaly_at_grid.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}")

reproject_geometry.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)

rm_z_coord.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

test/reproj_test.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/test_generate_grid.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_rotate_nu.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)