From 75dd9df68a2832882ac1dceb1659349780d558e1 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Tue, 3 Jan 2023 10:23:59 +0000 Subject: [PATCH 01/12] Added merging function for all match txt files --- s2p/__init__.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/s2p/__init__.py b/s2p/__init__.py index 3da0d407..633486f1 100644 --- a/s2p/__init__.py +++ b/s2p/__init__.py @@ -90,6 +90,15 @@ def check_missing_sift(tiles_pairs): print(f"WARNING: missing {len(missing_sift)}/{len(tiles_pairs)} " "SIFT matches, this may deteriorate output quality") +def merge_all_match_files(): + read_files = glob.glob("**/*sift_matches.txt", recursive=True) + with open(os.path.join(cfg["out_dir"], "merged_sift_matches.txt"), "w") as outfile: + for f in read_files: + if os.stat(f).st_size == 0: + continue + with open(f, "r") as infile: + outfile.write(infile.read()) + def pointing_correction(tile, i): """ From 7d8ff04bdefddc36a0dc19687ec87f98e30d355d Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Tue, 3 Jan 2023 10:24:40 +0000 Subject: [PATCH 02/12] added new imports for merge function --- s2p/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/s2p/__init__.py b/s2p/__init__.py index 633486f1..e1dd2a67 100644 --- a/s2p/__init__.py +++ b/s2p/__init__.py @@ -24,6 +24,8 @@ import json import datetime import multiprocessing +import glob +import os import numpy as np import rasterio From f01c76d4bc19cdf6b8b9bf02a896fed2a87d0a23 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Tue, 3 Jan 2023 10:39:57 +0000 Subject: [PATCH 03/12] pin numpy version --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index e85173f8..5433da16 100644 --- a/setup.py +++ b/setup.py @@ -42,7 +42,7 @@ def finalize_options(self): BdistWheel = None -requirements = ['numpy', +requirements = ['numpy==1.22.2', 'scipy', 'rasterio[s3] @ https://github.com/rasterio/rasterio/archive/refs/tags/1.3.3.tar.gz', 'utm', From 84a72327b9127b1aa7a9b029da6d0db472ad2c42 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 09:42:21 +0000 Subject: [PATCH 04/12] Added ability to correct point locations with rpc --- utils/image_coordinates_to_coordinates.py | 124 ++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 utils/image_coordinates_to_coordinates.py diff --git a/utils/image_coordinates_to_coordinates.py b/utils/image_coordinates_to_coordinates.py new file mode 100644 index 00000000..d0468128 --- /dev/null +++ b/utils/image_coordinates_to_coordinates.py @@ -0,0 +1,124 @@ +import glob +import os +import rasterio +import numpy as np +import geopandas as gpd +import shapely +from shapely import Point +import srtm4 +from rpcm import RPCModel, rpc_from_rpc_file +import geopy +import geopy.distance + + +def localize_row_col_geometry( + geometry: shapely.geometry.Point, + rpc_model: RPCModel, +): + """ + Project points in the row and column space to the geographic space, with altitude correction. + All the altitudes here are in meters. + Args: + geometry: The geometry to add altitudes to, should be a shapely point. + rpc_model: the RPCModel to base the projection on. + returns: + A tuple containing (new longitude, new latitude). + """ + points = [] + col_idx, row_idx = point.coords[0] + points.append((col_idx, row_idx)) + col_indices, row_indices = np.asarray(points).T + coords = [ + rpc_model.localization(col_idx, row_idx, 0) + for col_idx, row_idx in zip(col_indices, row_indices) + ] + corrected_lons, corrected_lats = [], [] + for (lon, lat), col_idx, row_idx in zip(coords, col_indices, row_indices): + + new_altitude = srtm4.srtm4(lon, lat) + new_lon, new_lat = 0, 0 + max_iterations = 10 + ground_dist_tol = 50 + ground_distance = ground_dist_tol + 1 + iteration_no = 0 + while (ground_distance > ground_dist_tol) and ( + iteration_no < max_iterations + ): + iteration_no += 1 + old_lon, old_lat = new_lon, new_lat + new_lon, new_lat = rpc_model.localization( + col_idx, row_idx, new_altitude + ) + new_altitude = srtm4.srtm4(new_lon, new_lat) + ground_distance = geopy.distance.distance( + (old_lat, old_lon), (new_lat, new_lon) + ).m + corrected_lons.append(new_lon) + corrected_lats.append(new_lat) + points.append( + ((lon, lat) for (lon, lat) in zip(corrected_lons, corrected_lats)) + ) + + corrected_lons.append(new_lon) + corrected_lats.append(new_lat) + points.append( + ((lon, lat) for (lon, lat) in zip(corrected_lons, corrected_lats)) + ) + return (new_lon, new_lat) + + +def read_and_prepare_matches(path_to_all_matches: str, match_columns: list=[0, 1]): + """ + Loads in the merged matches file and returns a geodataframe containing Point-like geometry. + + path_to_all_matches: The path to the text file containing all of the matches. + match_columns: The columns to be used as x, y for the matches geodataframe output. + """ + matches = np.loadtxt(path_to_all_matches) + matches = matches.astype(int) + matches = matches[:, match_columns] + matches_x = matches[:, 0] + matches_y = matches[:, 1] + matches_gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(matches_x, matches_y)) + return matches_gdf + + +def correct_points(matches_gdf: gpd.GeoDataFrame, rpc_model_path: str, every_x_match: int=10): + """ + matches_gdf: A geodataframe containing point geometries of the matches created by s2p in image coordinates. + rpc_model_path: A path to the rpc model to be used to calculate real coordinates from image coordinates. + every_x_match: For decreasing the size of the matches geodataframe, takes every point, i.e 5 is every 5. + + return: two numpy arrays, (lons and lats). + """ + corrected_points = [] + rpc_model = rpc_from_rpc_file(rpc_model_path) + matches_gdf = matches_gdf["geometry"][::every_x_match] + + # Iterate over matches and correct them. + for i, point in enumerate(matches_gdf): + localized_point = localize_row_col_geometry(point, rpc_model) + corr_points.append(list(localized_point)) + + # Prepare the lons and lats for returning + corr_points = np.array(corr_points) + lons = corr_points[:, 0] + lats = corr_points[:, 1] + + return lons, lats + + +def points_to_geojson(out_path: str, lons: np.ndarray, lats: np.ndarray): + """ + Saves the lons and lats from correct_points as a GeoJSON. + + out_path: The path to save the output GeoJSON to, + lons: The longitudes returned by correct_points. + lats: The latitudes returned by correct_points. + """ + gdf_corrected = gpd.GeoDataFrame(geometry=gpd.points_from_xy(lon, lat, crs="EPSG:4326")) + gdf_corrected.to_file(out_path, driver="GeoJSON") + + + + \ No newline at end of file From bd5c40d5787da2a2ee2841f24e22330e2e352d0a Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 09:59:25 +0000 Subject: [PATCH 05/12] Added geojson to matching --- s2p/__init__.py | 17 ++++++++++++++++- utils/image_coordinates_to_coordinates.py | 5 +++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/s2p/__init__.py b/s2p/__init__.py index e1dd2a67..f3d73ebb 100644 --- a/s2p/__init__.py +++ b/s2p/__init__.py @@ -46,6 +46,7 @@ from s2p import triangulation from s2p import fusion from s2p import visualisation +from utils.image_coordinates_to_coordinates import matches_to_geojson def remove_missing_tiles(tiles): @@ -92,7 +93,11 @@ def check_missing_sift(tiles_pairs): print(f"WARNING: missing {len(missing_sift)}/{len(tiles_pairs)} " "SIFT matches, this may deteriorate output quality") + def merge_all_match_files(): + """ + Merges together all non-empty match files into one .txt file and saves them in the output directory. + """ read_files = glob.glob("**/*sift_matches.txt", recursive=True) with open(os.path.join(cfg["out_dir"], "merged_sift_matches.txt"), "w") as outfile: for f in read_files: @@ -640,7 +645,16 @@ def main(user_cfg, start_from=0): print('2) correcting pointing globally...') global_pointing_correction(tiles) common.print_elapsed_time() - + + # Create matches GeoJSON. + print("Creating matches GeoJSON") + merge_all_match_files() + matches_to_geojson(f"{cfg['out_dir']}/merged_sift_matches.txt", + cfg['images'][0]['rpcm'], + 10, + [0, 1], + f"{cfg['out_dir']}/matches.geojson" + ) # rectification step: if start_from <= 3: print('3) rectifying tiles...') @@ -704,6 +718,7 @@ def main(user_cfg, start_from=0): if start_from <= 7: print('7) computing global DSM...') global_dsm(tiles) + common.print_elapsed_time() # cleanup diff --git a/utils/image_coordinates_to_coordinates.py b/utils/image_coordinates_to_coordinates.py index d0468128..91817b40 100644 --- a/utils/image_coordinates_to_coordinates.py +++ b/utils/image_coordinates_to_coordinates.py @@ -120,5 +120,10 @@ def points_to_geojson(out_path: str, lons: np.ndarray, lats: np.ndarray): gdf_corrected.to_file(out_path, driver="GeoJSON") +def matches_to_geojson(path_to_matches: str, path_to_rpc: str, every_x_match: int, match_columns: list, out_path: str): + matches_gdf = read_and_prepare_matches(path_to_matches, match_columns) + lons, lats = correct_points(matches_gdf, path_to_rpc, every_x_match) + points_to_geojson(out_path, lons, lats) + \ No newline at end of file From bcdf01acd40f85f51089ac3000ef2fcae4bd9973 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 10:18:29 +0000 Subject: [PATCH 06/12] Added missing reqs --- setup.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 5433da16..6d2074c9 100644 --- a/setup.py +++ b/setup.py @@ -53,7 +53,9 @@ def finalize_options(self): 'ransac', 'rpcm @ https://github.com/20treeAI/rpcm/archive/refs/tags/v1.4.8.tar.gz', 'srtm4 @ https://github.com/20treeAI/srtm4/archive/refs/tags/1.2.4.tar.gz', - 'requests'] + 'requests', + 'geopandas', + 'geopy'] extras_require = { "test": ["pytest", "pytest-cov", "psutil"], From fa7a0dd8cfa2f6f1684a6ce6628bf43a39874c05 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 10:48:38 +0000 Subject: [PATCH 07/12] bugfix --- utils/image_coordinates_to_coordinates.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/utils/image_coordinates_to_coordinates.py b/utils/image_coordinates_to_coordinates.py index 91817b40..77182f06 100644 --- a/utils/image_coordinates_to_coordinates.py +++ b/utils/image_coordinates_to_coordinates.py @@ -92,7 +92,10 @@ def correct_points(matches_gdf: gpd.GeoDataFrame, rpc_model_path: str, every_x_m return: two numpy arrays, (lons and lats). """ corrected_points = [] - rpc_model = rpc_from_rpc_file(rpc_model_path) + if type(rpc_model_path) == str: + rpc_model = rpc_from_rpc_file(rpc_model_path) + else: + rpc_model = rpc_model_path matches_gdf = matches_gdf["geometry"][::every_x_match] # Iterate over matches and correct them. From d25034aa516415e6c54cf6e9660b9cec425aada5 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 12:40:05 +0000 Subject: [PATCH 08/12] bugfix 2 --- utils/image_coordinates_to_coordinates.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/utils/image_coordinates_to_coordinates.py b/utils/image_coordinates_to_coordinates.py index 77182f06..3650e634 100644 --- a/utils/image_coordinates_to_coordinates.py +++ b/utils/image_coordinates_to_coordinates.py @@ -25,8 +25,6 @@ def localize_row_col_geometry( A tuple containing (new longitude, new latitude). """ points = [] - col_idx, row_idx = point.coords[0] - points.append((col_idx, row_idx)) col_indices, row_indices = np.asarray(points).T coords = [ rpc_model.localization(col_idx, row_idx, 0) From 66c45920fc7ec9c32777212e9f2cf9b7cb4864ee Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 12:56:55 +0000 Subject: [PATCH 09/12] use actual function input instead of point --- utils/image_coordinates_to_coordinates.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/utils/image_coordinates_to_coordinates.py b/utils/image_coordinates_to_coordinates.py index 3650e634..edb2199a 100644 --- a/utils/image_coordinates_to_coordinates.py +++ b/utils/image_coordinates_to_coordinates.py @@ -25,6 +25,8 @@ def localize_row_col_geometry( A tuple containing (new longitude, new latitude). """ points = [] + col_idx, row_idx = geometry.coords[0] + points.append((col_idx, row_idx)) col_indices, row_indices = np.asarray(points).T coords = [ rpc_model.localization(col_idx, row_idx, 0) From d8a0c77fcea672032a66d6374319a26df779fde4 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 13:42:29 +0000 Subject: [PATCH 10/12] fixed updated coords variable name --- utils/image_coordinates_to_coordinates.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/utils/image_coordinates_to_coordinates.py b/utils/image_coordinates_to_coordinates.py index edb2199a..cecc1a73 100644 --- a/utils/image_coordinates_to_coordinates.py +++ b/utils/image_coordinates_to_coordinates.py @@ -101,12 +101,12 @@ def correct_points(matches_gdf: gpd.GeoDataFrame, rpc_model_path: str, every_x_m # Iterate over matches and correct them. for i, point in enumerate(matches_gdf): localized_point = localize_row_col_geometry(point, rpc_model) - corr_points.append(list(localized_point)) + corrected_points.append(list(localized_point)) # Prepare the lons and lats for returning - corr_points = np.array(corr_points) - lons = corr_points[:, 0] - lats = corr_points[:, 1] + corrected_points = np.array(corrected_points) + lons = corrected_points[:, 0] + lats = corrected_points[:, 1] return lons, lats From b3ca250142bc47cb08536fc905c7bb2d3f832f0d Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Wed, 4 Jan 2023 14:17:06 +0000 Subject: [PATCH 11/12] fixed another variable name --- utils/image_coordinates_to_coordinates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/image_coordinates_to_coordinates.py b/utils/image_coordinates_to_coordinates.py index cecc1a73..d1ee3a10 100644 --- a/utils/image_coordinates_to_coordinates.py +++ b/utils/image_coordinates_to_coordinates.py @@ -119,7 +119,7 @@ def points_to_geojson(out_path: str, lons: np.ndarray, lats: np.ndarray): lons: The longitudes returned by correct_points. lats: The latitudes returned by correct_points. """ - gdf_corrected = gpd.GeoDataFrame(geometry=gpd.points_from_xy(lon, lat, crs="EPSG:4326")) + gdf_corrected = gpd.GeoDataFrame(geometry=gpd.points_from_xy(lons, lats, crs="EPSG:4326")) gdf_corrected.to_file(out_path, driver="GeoJSON") From e4621ef4e9c69f1fe491555e9fc52c07ffc026e4 Mon Sep 17 00:00:00 2001 From: Tim Iles Date: Thu, 5 Jan 2023 10:07:41 +0000 Subject: [PATCH 12/12] restructured to fix imports --- s2p/__init__.py | 19 ++++++++++--------- s2p/config.py | 3 +++ .../image_coordinates_to_coordinates.py | 0 3 files changed, 13 insertions(+), 9 deletions(-) rename {utils => s2p}/image_coordinates_to_coordinates.py (100%) diff --git a/s2p/__init__.py b/s2p/__init__.py index f3d73ebb..c025ce67 100644 --- a/s2p/__init__.py +++ b/s2p/__init__.py @@ -46,7 +46,7 @@ from s2p import triangulation from s2p import fusion from s2p import visualisation -from utils.image_coordinates_to_coordinates import matches_to_geojson +from s2p.image_coordinates_to_coordinates import matches_to_geojson def remove_missing_tiles(tiles): @@ -647,14 +647,15 @@ def main(user_cfg, start_from=0): common.print_elapsed_time() # Create matches GeoJSON. - print("Creating matches GeoJSON") - merge_all_match_files() - matches_to_geojson(f"{cfg['out_dir']}/merged_sift_matches.txt", - cfg['images'][0]['rpcm'], - 10, - [0, 1], - f"{cfg['out_dir']}/matches.geojson" - ) + if cfg["return_matches"] is True: + merge_all_match_files() + matches_to_geojson(f"{cfg['out_dir']}/merged_sift_matches.txt", + cfg['images'][0]['rpcm'], + 10, + [0, 1], + f"{cfg['out_dir']}/matches.geojson" + ) + # rectification step: if start_from <= 3: print('3) rectifying tiles...') diff --git a/s2p/config.py b/s2p/config.py index b1551eaf..78cf33d8 100644 --- a/s2p/config.py +++ b/s2p/config.py @@ -177,3 +177,6 @@ # If the out_crs is not set this parameter determines if the output CRS uses the EGM96 geoid vertical datum (if True) # or the WGS84 ellipsoid vertical datum (if False). If out_crs is set, this parameter is ignored. cfg['out_geoid'] = False + +# Create an output GeoJSON containing all of the matches created by sift.py +cfg["return_matches"] = False \ No newline at end of file diff --git a/utils/image_coordinates_to_coordinates.py b/s2p/image_coordinates_to_coordinates.py similarity index 100% rename from utils/image_coordinates_to_coordinates.py rename to s2p/image_coordinates_to_coordinates.py