From f5a778f5f4f5a188e6769e835fece9fb707500af Mon Sep 17 00:00:00 2001 From: Tommy Date: Mon, 23 Dec 2024 11:20:47 +0100 Subject: [PATCH 1/6] feat: add window for the use of mondial vrt --- cars/core/projection.py | 56 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 53 insertions(+), 3 deletions(-) diff --git a/cars/core/projection.py b/cars/core/projection.py index 01974cb8..9a6230c1 100644 --- a/cars/core/projection.py +++ b/cars/core/projection.py @@ -26,13 +26,15 @@ # Standard imports import logging +import math import os -from typing import List, Tuple +from typing import List, Tuple, Union # Third party imports import numpy as np import pandas import pyproj +import rasterio import rasterio as rio import xarray as xr from rasterio.features import shapes @@ -86,14 +88,31 @@ def compute_dem_intersection_with_poly( # noqa: C901 if ref_epsg != file_epsg: file_bb = polygon_projection(file_bb, file_epsg, ref_epsg) + poly_pixels = [] + for lon, lat in ref_poly.exterior.coords: + row, col = data.index(lon, lat) + poly_pixels.append((row, col)) + + poly_pixels = np.array(poly_pixels) + min_row, min_col = np.min(poly_pixels, axis=0) + max_row, max_col = np.max(poly_pixels, axis=0) + + window = rasterio.windows.Window( + col_off=int(min_col), + row_off=int(min_row), + width=int(max_col - min_col), + height=int(max_row - min_row), + ) + # if the srtm tile intersects the reference polygon if file_bb.intersects(ref_poly): local_dem_poly = None # retrieve valid polygons for poly, val in shapes( - data.dataset_mask(), - transform=data.transform, + data.read(1, window=window), + data.dataset_mask(window=window), + transform=data.window_transform(window), ): if val != 0: poly = shape(poly) @@ -528,6 +547,37 @@ def ground_intersection_envelopes( return inter_poly, (inter_xmin, inter_ymin, inter_xmax, inter_ymax) +def project_coordinates_on_line( + x_coord: Union[float, np.ndarray], + y_coord: Union[float, np.ndarray], + origin: Union[List[float], np.ndarray], + vec: Union[List[float], np.ndarray], +) -> Union[float, np.ndarray]: + """ + Project coordinates (x,y) on a line starting from origin with a + direction vector vec, and return the euclidean distances between + projected points and origin. + + :param x_coord: scalar or vector of coordinates x + :param y_coord: scalar or vector of coordinates x + :param origin: coordinates of origin point for line + :param vec: direction vector of line + :return: vector of distances of projected points to origin + """ + assert len(x_coord) == len(y_coord) + assert len(origin) == 2 + assert len(vec) == 2 + + vec_angle = math.atan2(vec[1], vec[0]) + point_angle = np.arctan2(y_coord - origin[1], x_coord - origin[0]) + proj_angle = point_angle - vec_angle + dist_to_origin = np.sqrt( + np.square(x_coord - origin[0]) + np.square(y_coord - origin[1]) + ) + + return dist_to_origin * np.cos(proj_angle) + + def get_ground_direction( sensor, geomodel, From 640eb98435b1663c62421b576d2b63f0fa628fbe Mon Sep 17 00:00:00 2001 From: Tommy Date: Mon, 23 Dec 2024 13:55:28 +0100 Subject: [PATCH 2/6] feat: new_version --- cars/core/projection.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/cars/core/projection.py b/cars/core/projection.py index 9a6230c1..617b06ba 100644 --- a/cars/core/projection.py +++ b/cars/core/projection.py @@ -88,14 +88,17 @@ def compute_dem_intersection_with_poly( # noqa: C901 if ref_epsg != file_epsg: file_bb = polygon_projection(file_bb, file_epsg, ref_epsg) - poly_pixels = [] - for lon, lat in ref_poly.exterior.coords: - row, col = data.index(lon, lat) - poly_pixels.append((row, col)) - - poly_pixels = np.array(poly_pixels) - min_row, min_col = np.min(poly_pixels, axis=0) - max_row, max_col = np.max(poly_pixels, axis=0) + min_lon, min_lat, max_lon, max_lat = ref_poly.bounds + + margin = 0.001 + + min_lon = min_lon - margin + max_lon = max_lon + margin + min_lat = min_lat - margin + max_lat = max_lat + margin + + min_row, min_col = data.index(min_lon, max_lat) + max_row, max_col = data.index(max_lon, min_lat) window = rasterio.windows.Window( col_off=int(min_col), From f896dcc15e7068c128763f0712434164323c4ac3 Mon Sep 17 00:00:00 2001 From: Tommy Date: Mon, 23 Dec 2024 15:26:35 +0100 Subject: [PATCH 3/6] fix: for pylint --- cars/core/projection.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/cars/core/projection.py b/cars/core/projection.py index 617b06ba..4bb18dae 100644 --- a/cars/core/projection.py +++ b/cars/core/projection.py @@ -34,7 +34,6 @@ import numpy as np import pandas import pyproj -import rasterio import rasterio as rio import xarray as xr from rasterio.features import shapes @@ -100,7 +99,7 @@ def compute_dem_intersection_with_poly( # noqa: C901 min_row, min_col = data.index(min_lon, max_lat) max_row, max_col = data.index(max_lon, min_lat) - window = rasterio.windows.Window( + window = rio.windows.Window( col_off=int(min_col), row_off=int(min_row), width=int(max_col - min_col), From f9f1a015e7ba68d55b99b43f8e9e8022c2f7a32c Mon Sep 17 00:00:00 2001 From: Tommy Date: Mon, 23 Dec 2024 15:46:12 +0100 Subject: [PATCH 4/6] fix: fix sonarqube alert --- sonar-project.properties | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sonar-project.properties b/sonar-project.properties index df3ed8be..f5fab1a6 100644 --- a/sonar-project.properties +++ b/sonar-project.properties @@ -26,7 +26,7 @@ sonar.python.xunit.reportPath=pytest-report.xml sonar.python.coverage.reportPaths=coverage.xml # deactivate duplicated lines in sonarqube in some chosen files: -sonar.cpd.exclusions=cars/pipelines/**/*_pipeline.py,cars/applications/point_cloud_outlier_removal/*, cars/applications/grid_generation/grid_correction.py, cars/orchestrator/log_wrapper.py +sonar.cpd.exclusions=cars/pipelines/**/*_pipeline.py,cars/applications/point_cloud_outlier_removal/*, cars/applications/grid_generation/grid_correction.py, cars/orchestrator/log_wrapper.py, cars/core/projection.py # Deactivate complexity rule for pipelines sonar.issue.ignore.multicriteria=complexity1,complexity2 From 59e8c7fba346df7dba31cc531bc10cfdfd940bdb Mon Sep 17 00:00:00 2001 From: Tommy Date: Thu, 2 Jan 2025 10:21:45 +0100 Subject: [PATCH 5/6] feat: modification in margin definition --- cars/core/projection.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/cars/core/projection.py b/cars/core/projection.py index 4bb18dae..5d657fc3 100644 --- a/cars/core/projection.py +++ b/cars/core/projection.py @@ -36,6 +36,7 @@ import pyproj import rasterio as rio import xarray as xr +from pyproj import CRS from rasterio.features import shapes from shapely.geometry import Polygon, shape from shapely.ops import transform @@ -87,17 +88,25 @@ def compute_dem_intersection_with_poly( # noqa: C901 if ref_epsg != file_epsg: file_bb = polygon_projection(file_bb, file_epsg, ref_epsg) - min_lon, min_lat, max_lon, max_lat = ref_poly.bounds - - margin = 0.001 - - min_lon = min_lon - margin - max_lon = max_lon + margin - min_lat = min_lat - margin - max_lat = max_lat + margin - - min_row, min_col = data.index(min_lon, max_lat) - max_row, max_col = data.index(max_lon, min_lat) + left, bottom, right, top = ref_poly.bounds + + # To ensure that the window used for extracting the SRTM + # polygons completely contains the ref_poly, + # a margin must be applied. Using the exact bounds of + # ref_poly could potentially lead to issues. + spatial_ref = CRS.from_epsg(ref_epsg) + if spatial_ref.is_geographic: + margin = 0.001 + else: + margin = 50 + + left = left - margin + right = right + margin + bottom = bottom - margin + top = top + margin + + min_row, min_col = data.index(left, top) + max_row, max_col = data.index(right, bottom) window = rio.windows.Window( col_off=int(min_col), From 94306ff949ac3d38df2e3721befffdfd0ea323c2 Mon Sep 17 00:00:00 2001 From: Tommy Date: Tue, 7 Jan 2025 09:07:35 +0100 Subject: [PATCH 6/6] fix: deletion of an unwanted function --- cars/core/projection.py | 34 +--------------------------------- sonar-project.properties | 2 +- 2 files changed, 2 insertions(+), 34 deletions(-) diff --git a/cars/core/projection.py b/cars/core/projection.py index 5d657fc3..597f308d 100644 --- a/cars/core/projection.py +++ b/cars/core/projection.py @@ -26,9 +26,8 @@ # Standard imports import logging -import math import os -from typing import List, Tuple, Union +from typing import List, Tuple # Third party imports import numpy as np @@ -558,37 +557,6 @@ def ground_intersection_envelopes( return inter_poly, (inter_xmin, inter_ymin, inter_xmax, inter_ymax) -def project_coordinates_on_line( - x_coord: Union[float, np.ndarray], - y_coord: Union[float, np.ndarray], - origin: Union[List[float], np.ndarray], - vec: Union[List[float], np.ndarray], -) -> Union[float, np.ndarray]: - """ - Project coordinates (x,y) on a line starting from origin with a - direction vector vec, and return the euclidean distances between - projected points and origin. - - :param x_coord: scalar or vector of coordinates x - :param y_coord: scalar or vector of coordinates x - :param origin: coordinates of origin point for line - :param vec: direction vector of line - :return: vector of distances of projected points to origin - """ - assert len(x_coord) == len(y_coord) - assert len(origin) == 2 - assert len(vec) == 2 - - vec_angle = math.atan2(vec[1], vec[0]) - point_angle = np.arctan2(y_coord - origin[1], x_coord - origin[0]) - proj_angle = point_angle - vec_angle - dist_to_origin = np.sqrt( - np.square(x_coord - origin[0]) + np.square(y_coord - origin[1]) - ) - - return dist_to_origin * np.cos(proj_angle) - - def get_ground_direction( sensor, geomodel, diff --git a/sonar-project.properties b/sonar-project.properties index f5fab1a6..df3ed8be 100644 --- a/sonar-project.properties +++ b/sonar-project.properties @@ -26,7 +26,7 @@ sonar.python.xunit.reportPath=pytest-report.xml sonar.python.coverage.reportPaths=coverage.xml # deactivate duplicated lines in sonarqube in some chosen files: -sonar.cpd.exclusions=cars/pipelines/**/*_pipeline.py,cars/applications/point_cloud_outlier_removal/*, cars/applications/grid_generation/grid_correction.py, cars/orchestrator/log_wrapper.py, cars/core/projection.py +sonar.cpd.exclusions=cars/pipelines/**/*_pipeline.py,cars/applications/point_cloud_outlier_removal/*, cars/applications/grid_generation/grid_correction.py, cars/orchestrator/log_wrapper.py # Deactivate complexity rule for pipelines sonar.issue.ignore.multicriteria=complexity1,complexity2