From 1812fc9864acb2af7f6c0ed0cc738cd17c753ab4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrico=20Cirac=C3=AC?= Date: Sun, 21 Jan 2024 17:14:06 +0100 Subject: [PATCH 1/7] environment.yaml - Updated --- environment.yaml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/environment.yaml b/environment.yaml index e612a90..15a9d28 100644 --- a/environment.yaml +++ b/environment.yaml @@ -19,5 +19,6 @@ dependencies: - xarray - flake8 - pytest - - pydantic - - pytest-cov \ No newline at end of file + - pytest-cov + - codecov + - coverage \ No newline at end of file From 498f07deb695a6eb8e167745efd6b086b31387c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrico=20Cirac=C3=AC?= Date: Sun, 21 Jan 2024 17:14:51 +0100 Subject: [PATCH 2/7] ptsline_test.py - unit test for PtsLine.py added. --- PtsLine.py | 55 ++++++++++++++++++++++++++------------------ test/ptsline_test.py | 40 ++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 22 deletions(-) create mode 100644 test/ptsline_test.py diff --git a/PtsLine.py b/PtsLine.py index fd39a53..ca4ade9 100644 --- a/PtsLine.py +++ b/PtsLine.py @@ -13,20 +13,22 @@ class PtsLine: """ Class to compute the equation of a line passing through two points. """ - def __init__(self, x1: float, y1: float, - x2: float, y2: float) -> None: - self.x1 = x1 - self.y1 = y1 - self.x2 = x2 - self.y2 = y2 - self.m = (y2 - y1) / (x2 - x1) - self.q = y1 - (self.m * x1) + def __init__(self, x_pt1: float, y_pt1: float, + x_pt2: float, y_pt2: float) -> None: + self.x_1 = x_pt1 + self.y_1 = y_pt1 + self.x_2 = x_pt2 + self.y_2 = y_pt2 + self.m_val = (y_pt2 - y_pt1) / (x_pt2 - x_pt1) + self.q_val = y_pt1 - (self.m_val * x_pt1) - def y(self, x: float | np.ndarray) -> float | np.ndarray: - return self.m * x + self.q + def y_val(self, x_pt: float | np.ndarray) -> float | np.ndarray: + """Return the y coordinate of the line at the given x coordinate.""" + return self.m_val * x_pt + self.q_val - def x(self, y: float | np.ndarray) -> float | np.ndarray: - return (y - self.q) / self.m + def x_val(self, y_pt: float | np.ndarray) -> float | np.ndarray: + """Return the x coordinate of the line at the given y coordinate.""" + return (y_pt - self.q_val) / self.m_val class PtsLineIntersect: @@ -35,13 +37,22 @@ class PtsLineIntersect: two points. """ def __init__(self, line_1: PtsLine, line_2: PtsLine) -> None: - self.m1 = line_1.m - self.q1 = line_1.q - self.m2 = line_2.m - self.q2 = line_2.q - - def intersect(self) -> tuple[float, float]: - x = (self.q2 - self.q1) / (self.m1 - self.m2) - y = self.m1 * x + self.q1 - return x, y - + self.m_1 = line_1.m_val + self.q_1 = line_1.q_val + self.m_2 = line_2.m_val + self.q_2 = line_2.q_val + + @property + def intersection(self) -> tuple[float, float]: + """ + Compute the intersection point of two lines. + Returns: Coordinates of the intersection point tuple[float, float]. + """ + x_c = None + y_c = None + try: + x_c = (self.q_2 - self.q_1) / (self.m_1 - self.m_2) + y_c = self.m_1 * x_c + self.q_1 + except ZeroDivisionError: + print("# - Lines are parallel.") + return x_c, y_c diff --git a/test/ptsline_test.py b/test/ptsline_test.py new file mode 100644 index 0000000..ee125da --- /dev/null +++ b/test/ptsline_test.py @@ -0,0 +1,40 @@ +#!/usr/bin/env python +""" +Unit tests for the PtsLine and PtsLineIntersect classes. +""" +import pytest +import numpy as np +from PtsLine import PtsLine, PtsLineIntersect + +@pytest.fixture +def example_lines(): + """Return two lines that intersect at (1, 1).""" + line1 = PtsLine(0, 0, 2, 2) + line2 = PtsLine(0, 2, 2, 0) + return line1, line2 + + +def test_pts_line_y_val(): + """Test the y_val method of the PtsLine class.""" + line = PtsLine(0, 0, 2, 2) + assert line.y_val(1) == 1 + + +def test_pts_line_x_val(): + """Test the x_val method of the PtsLine class.""" + line = PtsLine(0, 0, 2, 2) + assert line.x_val(1) == 1 + + +def test_pts_line_intersect(example_lines): + """Test the intersection property of the PtsLineIntersect class.""" + line1, line2 = example_lines + intersect = PtsLineIntersect(line1, line2) + assert np.allclose(intersect.intersection, (1, 1)) + + +def test_parallel_lines(): + line1 = PtsLine(0, 0, 2, 2) + line2 = PtsLine(1, 1, 3, 3) + intersect = PtsLineIntersect(line1, line2) + assert intersect.intersection == (None, None) \ No newline at end of file From 7c544f0f79b6edb064925b0de81a8272163e9ea5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrico=20Cirac=C3=AC?= Date: Sun, 21 Jan 2024 17:15:37 +0100 Subject: [PATCH 3/7] generate_grid.py - complete description of the processing steps added. --- generate_grid.py | 173 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 119 insertions(+), 54 deletions(-) diff --git a/generate_grid.py b/generate_grid.py index a9634ea..7828df5 100644 --- a/generate_grid.py +++ b/generate_grid.py @@ -3,7 +3,19 @@ Written by Enrico Ciraci' January 2024 -Test along/track grid generation procedure. +Compute a regular grid along the provided satellite track. + + 1. Merge frames polygons into a single track polygon. + 2. Compute the centroid of the track polygon. + 3. Project the track polygon to 3857 Web Mercator Projection. + 4. Rotate the track polygon to align it with the North-South direction. + 5. Compute the trapezoid corners coordinates. + 6. Compute the trapezoid diagonals equations. + 7. Extend diagonals using a user define buffer. + 8. Split the vertical and horizontal dimensions into a number of segments + defined by az_res and n_c parameters. + 9. Compute the grid cells corners coordinates. + 10. Generate output shapefile. positional arguments: input_file Input file. @@ -14,32 +26,66 @@ Output directory. --buffer_dist BUFFER_DIST, -B BUFFER_DIST Buffer distance. + --az_res AZ_RES, -R AZ_RES + Cross track grid resolution (m). + --n_c N_C, -C N_C Number of columns in the grid. --plot, -P Plot intermediate results. + Python Dependencies geopandas: Open source project to make working with geospatial data in python easier: https://geopandas.org pyproj: Python interface to PROJ (cartographic projections and coordinate transformations library): https://pyproj4.github.io/pyproj/stable/index.html -shapely: Python package for manipulation and analysis of planar geometric +shapely: Python package for manipulation and analy_sis of planar geometric objects: https://shapely.readthedocs.io/en/stable/ +matplotlib: Comprehensive library for creating static, animated, and + interactive visualizations in Python: + https://matplotlib.org + """ # - Python Dependencies from __future__ import print_function import os import argparse -import numpy as np from datetime import datetime +import numpy as np import geopandas as gpd -from shapely.geometry import Polygon +from shapely.geometry import Polygon, Point from shapely.affinity import rotate +from pyproj import CRS, Transformer +from shapely.ops import transform import matplotlib.pyplot as plt from iride_csk_frame_grid_utils import (reproject_geodataframe, rotate_polygon_to_north_up) from PtsLine import PtsLine, PtsLineIntersect +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) + + def main() -> None: """ Test along/track grid generation procedure. @@ -70,40 +116,59 @@ def main() -> None: # - Parse arguments args = parser.parse_args() + # - Number of Cells + n_c = args.n_c # - Along Track - Number of Columns + az_res = args.az_res # - Cross Track Resolution (km) - Azimuth Resolution + # - set path to input shapefile input_shapefile = args.input_file - #input_shapefile \ - # = os.path.join('.', 'data', 'csk_frame_map_italy_epsg4326_dissolve', - # 'csk_frame_map_italy_epsg4326_dissolve.shp') output_f_name \ = os.path.basename(input_shapefile).replace('.shp','_grid.shp') # - set path to output shapefile out_dir = args.out_dir - #out_dir = os.path.join(r'C:\Users\e.ciraci\Desktop\test') os.makedirs(out_dir, exist_ok=True) # - import input data gdf = gpd.read_file(input_shapefile) # - get input data crs source_crs = gdf.crs.to_epsg() + # - extract only needed columns + gdf = gdf[['geometry']] - # - Number of Cells - n_c = args.n_c # - Along Track - Number of Columns - az_res = args.az_res # - Cross Track Resolution (km) - Azimuth Resolution + # - 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 + + # - Merge frames polygons into a single track polygon + gdf_t \ + = (gpd.GeoDataFrame(geometry=[gdf.unary_union], crs=source_crs) + .explode(index_parts=False).reset_index(drop=True)) + + # - Compute Track Centroid - Need to reproject the track geometry + # - to minimize distortion in the calculation. + # - 1. Project to WGS 84 Web Mercator Projection EPSG:3857 + # - 2. Compute Centroid + r_geom \ + = reproject_geometry(gdf_t['geometry'], + source_crs, 3857).loc[0] + proj_centroid = Point(r_geom.centroid.x, r_geom.centroid.y) + ll_centroid = reproject_geometry(proj_centroid, 3857, source_crs) - # - Extract Polygon centroid - ll_centroid = (gdf['geometry'].centroid.x, gdf['geometry'].centroid.y) # - Longitude of the centroid - convert to radians - lat_cent = ll_centroid[1][0] * np.pi / 180 + lat_cent = ll_centroid.y * np.pi / 180 # - Estimate an average distortion factor associated to the # - usage of Web Mercator Projection (EPSG:3857) # - Reference: https://en.wikipedia.org/wiki/Mercator_projection d_scale = np.cos(lat_cent) # - reproject to WGS 84 Web Mercator Projection EPSG:3857 - gdf = reproject_geodataframe(gdf, 3857) + gdf = reproject_geodataframe(gdf_t, 3857) - # - remove z coordinate + # - If still present remove z coordinate d3_coord = gdf['geometry'].loc[0].exterior.coords[:-1] d2_coords = Polygon([(coord[0], coord[1]) for coord in d3_coord]) @@ -143,44 +208,44 @@ def main() -> None: # - Compute trapezoid diagonals equations # - Diagonal 1 - ln1 = PtsLine(pt_ul[0], pt_ul[1], pt_lr[0], pt_lr[1]) + ln_1 = PtsLine(pt_ul[0], pt_ul[1], pt_lr[0], pt_lr[1]) # - Diagonal 2 - ln2 = PtsLine(pt_ur[0], pt_ur[1], pt_ll[0], pt_ll[1]) + ln_2 = PtsLine(pt_ur[0], pt_ur[1], pt_ll[0], pt_ll[1]) # - Extend diagonals using a user define buffer - xs = [] - ys = [] + x_s = [] + y_s = [] # - Corner 1 - Upper Left - x1 = pt_ul[0] - args.buffer_dist - y1 = ln1.y(x1) - ul_ext = (x1, y1) - xs.append(x1) - ys.append(y1) + x_1 = pt_ul[0] - args.buffer_dist + y_1 = ln_1.y_val(x_1) + ul_ext = (x_1, y_1) + x_s.append(x_1) + y_s.append(y_1) # - Corner 2 - x2 = pt_ur[0] + args.buffer_dist - y2 = ln2.y(x2) - ur_ext = (x2, y2) - xs.append(x2) - ys.append(y2) + x_2 = pt_ur[0] + args.buffer_dist + y_2 = ln_2.y_val(x_2) + ur_ext = (x_2, y_2) + x_s.append(x_2) + y_s.append(y_2) # - Corner 3 - x3 = pt_lr[0] + args.buffer_dist - y3 = ln1.y(x3) - lr_ext = (x3, y3) - xs.append(x3) - ys.append(y3) + x_3 = pt_lr[0] + args.buffer_dist + y_3 = ln_1.y_val(x_3) + lr_ext = (x_3, y_3) + x_s.append(x_3) + y_s.append(y_3) # - Corner 4 - x4 = pt_ll[0] - args.buffer_dist - y4 = ln2.y(x4) - ll_ext = (x4, y4) - xs.append(x4) - ys.append(y4) - xs.append(x1) - ys.append(y1) + x_4 = pt_ll[0] - args.buffer_dist + y_4 = ln_2.y_val(x_4) + ll_ext = (x_4, y_4) + x_s.append(x_4) + y_s.append(y_4) + x_s.append(x_1) + y_s.append(y_1) # - Trapezoid major axis equation - ln3 = PtsLine(ul_ext[0], ul_ext[1], ur_ext[0], ur_ext[1]) + ln_3 = PtsLine(ul_ext[0], ul_ext[1], ur_ext[0], ur_ext[1]) # - Trapezoid minor axis equation - ln4 = PtsLine(ll_ext[0], ll_ext[1], lr_ext[0], lr_ext[1]) + ln_4 = PtsLine(ll_ext[0], ll_ext[1], lr_ext[0], lr_ext[1]) # - Trapezoid left side equation ln5 = PtsLine(ul_ext[0], ul_ext[1], ll_ext[0], ll_ext[1]) # - Trapezoid right side equation @@ -193,8 +258,8 @@ def main() -> None: x_south = np.linspace(ll_ext[0], lr_ext[0], n_c+1) # - Evaluate the y coordinates of the grid vertical lines - y_north = ln3.y(x_north) - y_south = ln4.y(x_south) + y_north = ln_3.y_val(x_north) + y_south = ln_4.y_val(x_south) # - Compute grid number of rows n_r = int(np.ceil(((max(y_north) - min(y_south)) * d_scale) / az_res)) @@ -207,9 +272,9 @@ def main() -> None: x_vert_l = [] x_vert_r = [] - for y in y_vert: - x_vert_l.append(ln5.x(y)) - x_vert_r.append(ln6.x(y)) + for y_p in y_vert: + x_vert_l.append(ln5.x_val(y_p)) + x_vert_r.append(ln6.x_val(y_p)) # - Generate list of vertical and horizontal lines horiz_lines = [PtsLine(x_vert_l[i], float(y_vert[i]), @@ -221,10 +286,10 @@ def main() -> None: # - Initialize Grid Corner Matrix corners_matrix = [] - for hr in horiz_lines: + for h_r in horiz_lines: corners_horiz = [] - for vr in vert_lines: - corners_horiz.append(PtsLineIntersect(hr, vr).intersect()) + for v_r in vert_lines: + corners_horiz.append(PtsLineIntersect(h_r, v_r).intersection) corners_matrix.append(corners_horiz) # - Convert to numpy array corners_matrix = np.array(corners_matrix) @@ -252,7 +317,7 @@ def main() -> None: if args.plot: # - Plot rotated geometry - fig, ax = plt.subplots(figsize=(5, 7)) + _, ax = plt.subplots(figsize=(5, 7)) ax.set_title('Rotated Geometry') ax.set_xlabel('Easting') ax.set_ylabel('Northing') @@ -265,7 +330,7 @@ def main() -> None: ax.scatter(ur_ext[0], ur_ext[1], color='red') ax.scatter(lr_ext[0], lr_ext[1], color='red') ax.scatter(ll_ext[0], ll_ext[1], color='red') - ax.scatter(xs, ys, color='orange', marker='x') + ax.scatter(x_s, y_s, color='orange', marker='x') ax.plot([pt_ul[0], pt_lr[0]], [pt_ul[1], pt_lr[1]], color='cyan') ax.scatter(x_north, y_north, color='green') ax.scatter(x_south, y_south, color='green') From abd7232af3273cb751295dcc367d9f4fc011f577 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrico=20Cirac=C3=AC?= Date: Sun, 21 Jan 2024 17:16:20 +0100 Subject: [PATCH 4/7] reproj_test.py - unit test for generate_grid.reproject_geometry added. --- test/reproj_test.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 test/reproj_test.py diff --git a/test/reproj_test.py b/test/reproj_test.py new file mode 100644 index 0000000..a1fdba2 --- /dev/null +++ b/test/reproj_test.py @@ -0,0 +1,29 @@ +#!/usr/bin/env python +""" Unit tests for the reproject_geometry function. """ +import pytest +from shapely.geometry import Point, Polygon +from generate_grid import reproject_geometry + +@pytest.fixture +def example_polygon(): + """Return a square polygon.""" + return Polygon([(0, 0), (0, 1), (1, 1), (1, 0)]) + +@pytest.fixture +def example_point(): + """Return a point in the center of the square polygon.""" + return Point(0.5, 0.5) + + +def test_reproject_geometry_polygon(example_polygon): + """Test the reproject_geometry function with a polygon.""" + result \ + = reproject_geometry(example_polygon, 4326, 3857) + assert isinstance(result, Polygon) + + +def test_reproject_geometry_point(example_point): + """Test the reproject_geometry function with a point.""" + result \ + = reproject_geometry(example_point,4326, 3857) + assert isinstance(result, Point) \ No newline at end of file From 2a0bd884c272762a36a69deb17a05925eb4ecdd9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrico=20Cirac=C3=AC?= Date: Sun, 21 Jan 2024 17:16:42 +0100 Subject: [PATCH 5/7] unit test added to GitHib workflow. --- .github/workflows/preliminary_test_conda.yaml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/preliminary_test_conda.yaml b/.github/workflows/preliminary_test_conda.yaml index ac5331f..ae0ed14 100644 --- a/.github/workflows/preliminary_test_conda.yaml +++ b/.github/workflows/preliminary_test_conda.yaml @@ -58,3 +58,7 @@ jobs: flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + - name: Run Unit test with pytest + shell: bash -l {0} + run: | + python -m pytest --import-mode=append test/ \ No newline at end of file From de6158e749178cf1ee6e07a095ddd443ab182971 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrico=20Cirac=C3=AC?= Date: Sun, 21 Jan 2024 17:40:48 +0100 Subject: [PATCH 6/7] Minors change. --- generate_grid.py | 1 + 1 file changed, 1 insertion(+) diff --git a/generate_grid.py b/generate_grid.py index 7828df5..cb68979 100644 --- a/generate_grid.py +++ b/generate_grid.py @@ -17,6 +17,7 @@ 9. Compute the grid cells corners coordinates. 10. Generate output shapefile. + positional arguments: input_file Input file. From 07d0e7ee1d22380d5dc326311fac88aaaf21b137 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Enrico=20Cirac=C3=AC?= Date: Sun, 21 Jan 2024 17:55:19 +0100 Subject: [PATCH 7/7] Running test with GitHub workflow with python 3.10 and 3.11. --- .github/workflows/preliminary_test_conda.yaml | 2 +- environment.yaml | 3 --- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/.github/workflows/preliminary_test_conda.yaml b/.github/workflows/preliminary_test_conda.yaml index ae0ed14..f120c59 100644 --- a/.github/workflows/preliminary_test_conda.yaml +++ b/.github/workflows/preliminary_test_conda.yaml @@ -15,7 +15,7 @@ jobs: - os: ubuntu-latest label: linux-64 prefix: /usr/share/miniconda3/envs/my-env - python-version: 3.9 + python-version: '3.11' - os: ubuntu-latest label: linux-64b diff --git a/environment.yaml b/environment.yaml index 15a9d28..b885041 100644 --- a/environment.yaml +++ b/environment.yaml @@ -5,7 +5,6 @@ dependencies: - cartopy - gdal - matplotlib - - netCDF4 - numpy - pip - pyproj @@ -15,8 +14,6 @@ dependencies: - pandas - geopandas - shapely - - rasterio - - xarray - flake8 - pytest - pytest-cov