From 3fb032388e2490711772bc51984258da7d50e922 Mon Sep 17 00:00:00 2001 From: Caitlin Adams <25995927+caitlinadams@users.noreply.github.com> Date: Fri, 17 Jan 2025 10:27:17 +1100 Subject: [PATCH] Initial reformat with black --- check_gamma_modules.py | 34 +- .../nci/preparation/create_config.py | 62 +-- .../nci/preparation/create_dem_vrt.py | 27 +- sar_antarctica/nci/preparation/dem.py | 298 ++++++++------ sar_antarctica/nci/preparation/geoid.py | 58 +-- sar_antarctica/nci/preparation/orbits.py | 38 +- sar_antarctica/nci/preparation/scenes.py | 35 +- .../nci/processing/GAMMA/GAMMA_utils.py | 3 +- .../nci/processing/pyroSAR/pyrosar_geocode.py | 67 ++-- .../submission/pyrosar_gamma/pyrosar_gamma.py | 46 ++- sar_antarctica/utils/raster.py | 211 +++++----- sar_antarctica/utils/rio_tools.py | 113 +++--- sar_antarctica/utils/spatial.py | 34 +- tests/filesystem/test_filesystem.py | 54 ++- tests/sar_antarctica/test_dem.py | 378 +++++++++++------- tests/sar_antarctica/test_dem_vrt.py | 20 +- tests/sar_antarctica/test_orbits.py | 15 +- tests/sar_antarctica/test_scenes.py | 21 +- utils.py | 22 +- 19 files changed, 945 insertions(+), 591 deletions(-) diff --git a/check_gamma_modules.py b/check_gamma_modules.py index ab2b3a2..361290e 100644 --- a/check_gamma_modules.py +++ b/check_gamma_modules.py @@ -7,7 +7,7 @@ REQUIRED_LIBS_PATH = "/g/data/yp75/projects/pyrosar_processing/sar-pyrosar-nci" if os.environ.get("GAMMA_HOME", None) is None: - + print(f"Setting GAMMA to {GAMMA_HOME_PATH}") os.environ["GAMMA_HOME"] = GAMMA_HOME_PATH @@ -21,23 +21,39 @@ else: print(os.environ.get("LD_LIBRARY_PATH")) -commands_to_skip = ["coord_trans", "phase_sum", "dishgt", "SLC_cat", "ras_ratio_dB", "ptarg"] +commands_to_skip = [ + "coord_trans", + "phase_sum", + "dishgt", + "SLC_cat", + "ras_ratio_dB", + "ptarg", +] -for module in finder(GAMMA_HOME_PATH, ['[A-Z]*'], foldermode=2): +for module in finder(GAMMA_HOME_PATH, ["[A-Z]*"], foldermode=2): print(module) - for submodule in ['bin', 'scripts']: + for submodule in ["bin", "scripts"]: print(f"{module}/{submodule}") - for cmd in sorted(finder(module+"/"+submodule, [r'^\w+$'], regex=True), key=lambda s: s.lower()): + for cmd in sorted( + finder(module + "/" + submodule, [r"^\w+$"], regex=True), + key=lambda s: s.lower(), + ): command_base = os.path.basename(cmd) if command_base in commands_to_skip: print(" skipping " + command_base) else: - + print(" " + command_base) - proc = sp.Popen(cmd, stdin=sp.PIPE, stdout=sp.PIPE, stderr=sp.PIPE, universal_newlines=True) + proc = sp.Popen( + cmd, + stdin=sp.PIPE, + stdout=sp.PIPE, + stderr=sp.PIPE, + universal_newlines=True, + ) out, err = proc.communicate() out += err - usage = re.search('usage:.*(?=\n)', out) + usage = re.search("usage:.*(?=\n)", out) if usage is None: print(" " + err) - print("\n\n") \ No newline at end of file + print("\n\n") diff --git a/sar_antarctica/nci/preparation/create_config.py b/sar_antarctica/nci/preparation/create_config.py index 8164f53..3628abb 100644 --- a/sar_antarctica/nci/preparation/create_config.py +++ b/sar_antarctica/nci/preparation/create_config.py @@ -7,7 +7,16 @@ from sar_antarctica.nci.preparation.orbits import find_latest_orbit_for_scene from sar_antarctica.nci.preparation.dem import get_cop30_dem_for_bounds -def write_file_paths(config_file: Path, scene_file: Path, orbit_file: Path, dem_file: Path, data_dir: Path, ancillary_dir="ancillary", processed_dir="processed_scene"): + +def write_file_paths( + config_file: Path, + scene_file: Path, + orbit_file: Path, + dem_file: Path, + data_dir: Path, + ancillary_dir="ancillary", + processed_dir="processed_scene", +): inputs_header = "[inputs]\n" scene_setting = f"scene = '{str(scene_file)}'\n" orbit_setting = f"orbit = '{str(orbit_file)}'\n" @@ -19,21 +28,24 @@ def write_file_paths(config_file: Path, scene_file: Path, orbit_file: Path, dem_ processed_setting = f"processed = '{processed_dir}'\n" with open(config_file, "w") as cf: - cf.writelines([ - inputs_header, - scene_setting, - orbit_setting, - dem_setting, - "\n", - outputs_header, - data_path_setting, - ancillary_setting, - processed_setting - ]) - + cf.writelines( + [ + inputs_header, + scene_setting, + orbit_setting, + dem_setting, + "\n", + outputs_header, + data_path_setting, + ancillary_setting, + processed_setting, + ] + ) + + @click.command() @click.argument("scene_id", nargs=1) -@click.argument('scene_config', nargs=1) +@click.argument("scene_config", nargs=1) def main(scene_id: str, scene_config: str): print(f"Processing scene: {scene_id} \n") @@ -52,22 +64,24 @@ def main(scene_id: str, scene_config: str): # Identify bounds of scene and use bounding box to build DEM scene = identify(str(scene_file)) scene_bbox = scene.bbox().extent - scene_bounds = (scene_bbox["xmin"], scene_bbox["ymin"], scene_bbox["xmax"], scene_bbox["ymax"]) + scene_bounds = ( + scene_bbox["xmin"], + scene_bbox["ymin"], + scene_bbox["xmax"], + scene_bbox["ymax"], + ) # Set path for dem and create dem_dir = data_dir / "dem" dem_file = dem_dir / f"{scene_id}_dem.tif" - _, _ = get_cop30_dem_for_bounds(bounds=scene_bounds, save_path=dem_file, ellipsoid_heights=True) + _, _ = get_cop30_dem_for_bounds( + bounds=scene_bounds, save_path=dem_file, ellipsoid_heights=True + ) # Write to config file - write_file_paths( - config_file, - scene_file, - latest_poe_file, - dem_file, - data_dir - ) + write_file_paths(config_file, scene_file, latest_poe_file, dem_file, data_dir) + if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/sar_antarctica/nci/preparation/create_dem_vrt.py b/sar_antarctica/nci/preparation/create_dem_vrt.py index babfb15..38af10e 100644 --- a/sar_antarctica/nci/preparation/create_dem_vrt.py +++ b/sar_antarctica/nci/preparation/create_dem_vrt.py @@ -2,6 +2,7 @@ from pathlib import Path from typing import Generator + def find_tiles(source_dir: Path, pattern: str) -> Generator[Path, None, None]: """_summary_ @@ -24,7 +25,12 @@ def find_tiles(source_dir: Path, pattern: str) -> Generator[Path, None, None]: return tiles -def build_vrt(tiles: Generator[Path, None, None] | list[Path], vrt_path: str | os.PathLike, run: bool = True): + +def build_vrt( + tiles: Generator[Path, None, None] | list[Path], + vrt_path: str | os.PathLike, + run: bool = True, +): """Generic function for building a VRT from a generator of tile paths Parameters @@ -42,11 +48,16 @@ def build_vrt(tiles: Generator[Path, None, None] | list[Path], vrt_path: str | o f.writelines(f"{tile}\n" for tile in tiles) if run: - os.system(f'gdalbuildvrt -input_file_list temp.txt {vrt_path}') + os.system(f"gdalbuildvrt -input_file_list temp.txt {vrt_path}") os.remove("temp.txt") -def build_tileindex(tiles: Generator[Path, None, None] | list[Path], tindex_path: str | os.PathLike, run: bool = True): + +def build_tileindex( + tiles: Generator[Path, None, None] | list[Path], + tindex_path: str | os.PathLike, + run: bool = True, +): """Generic function for building a tile index from a generator of tile paths Parameters @@ -64,18 +75,18 @@ def build_tileindex(tiles: Generator[Path, None, None] | list[Path], tindex_path f.writelines(f"{tile}\n" for tile in tiles) if run: - os.system(f'gdaltindex {tindex_path} --optfile temp.txt') + os.system(f"gdaltindex {tindex_path} --optfile temp.txt") os.remove("temp.txt") + def create_glo30_dem_south_vrt(): - """Create a VRT for the Copernicus Global 30m DEM on NCI - """ - + """Create a VRT for the Copernicus Global 30m DEM on NCI""" + SOURCE_DIR = Path("/g/data/v10/eoancillarydata-2/elevation/copernicus_30m_world") PATTERN = "Copernicus_DSM_COG_10_S??_00_????_00_DEM/*.tif" VRT_PATH = Path("/g/data/yp75/projects/ancillary/dem/copdem_south.vrt") tiles = find_tiles(SOURCE_DIR, PATTERN) - build_vrt(tiles, VRT_PATH) \ No newline at end of file + build_vrt(tiles, VRT_PATH) diff --git a/sar_antarctica/nci/preparation/dem.py b/sar_antarctica/nci/preparation/dem.py index 23fd7ec..4647ca8 100644 --- a/sar_antarctica/nci/preparation/dem.py +++ b/sar_antarctica/nci/preparation/dem.py @@ -10,31 +10,36 @@ from .geoid import remove_geoid from ...utils.raster import ( - expand_raster_to_bounds, + expand_raster_to_bounds, reproject_raster, merge_raster_files, bounds_from_profile, - read_vrt_in_bounds + read_vrt_in_bounds, ) from ...utils.spatial import ( - adjust_bounds, + adjust_bounds, get_local_utm, ) -COP30_FOLDER_PATH = Path('/g/data/v10/eoancillarydata-2/elevation/copernicus_30m_world/') -GEOID_TIF_PATH = Path('/g/data/yp75/projects/ancillary/geoid/us_nga_egm2008_1_4326__agisoft.tif') -COP30_GPKG_PATH = Path('/g/data/yp75/projects/ancillary/dem/copdem_tindex.gpkg') +COP30_FOLDER_PATH = Path( + "/g/data/v10/eoancillarydata-2/elevation/copernicus_30m_world/" +) +GEOID_TIF_PATH = Path( + "/g/data/yp75/projects/ancillary/geoid/us_nga_egm2008_1_4326__agisoft.tif" +) +COP30_GPKG_PATH = Path("/g/data/yp75/projects/ancillary/dem/copdem_tindex.gpkg") + def get_cop30_dem_for_bounds( - bounds: tuple, - save_path: Path, - ellipsoid_heights: bool = True, - buffer_pixels : int = 1, - adjust_for_high_lat_and_buffer = True, - cop30_index_path : Path = COP30_GPKG_PATH, - cop30_folder_path : Path = COP30_FOLDER_PATH, - geoid_tif_path : Path = GEOID_TIF_PATH, - ) -> tuple[np.ndarray, dict]: + bounds: tuple, + save_path: Path, + ellipsoid_heights: bool = True, + buffer_pixels: int = 1, + adjust_for_high_lat_and_buffer=True, + cop30_index_path: Path = COP30_GPKG_PATH, + cop30_folder_path: Path = COP30_FOLDER_PATH, + geoid_tif_path: Path = GEOID_TIF_PATH, +) -> tuple[np.ndarray, dict]: """Logic for acquiting the cop30m DEM for a given set of bounds on the NCI. The returned dem will fully encompass the specified bounds. There may be additional data outside of the bounds as all data from the merged tiles is returned. @@ -63,114 +68,145 @@ def get_cop30_dem_for_bounds( dem array and dem rasterio profile """ - assert cop30_index_path or cop30_folder_path, \ - 'Either a `cop30_index_path` or `cop30_folder_path` must be provided' + assert ( + cop30_index_path or cop30_folder_path + ), "Either a `cop30_index_path` or `cop30_folder_path` must be provided" if cop30_index_path and cop30_folder_path: - logging.info(f'both `cop30_index_path` and `cop30_folder_path` provided. `cop30_index_path` used by default') + logging.info( + f"both `cop30_index_path` and `cop30_folder_path` provided. `cop30_index_path` used by default" + ) - logging.info(f'Getting cop30m dem that covers bounds: {bounds}') + logging.info(f"Getting cop30m dem that covers bounds: {bounds}") # check if scene crosses the AM - antimeridian_crossing = check_s1_bounds_cross_antimeridian(bounds, max_scene_width=8) + antimeridian_crossing = check_s1_bounds_cross_antimeridian( + bounds, max_scene_width=8 + ) if antimeridian_crossing: - logging.warning('DEM crosses the dateline/antimeridian') - logging.info('Finding best crs for area') + logging.warning("DEM crosses the dateline/antimeridian") + logging.info("Finding best crs for area") target_crs = get_target_antimeridian_projection(bounds) - logging.warning(f'Data will be returned in EPSG:{target_crs} projection') + logging.warning(f"Data will be returned in EPSG:{target_crs} projection") # split the scene into left and right - logging.info(f'Splitting bounds into left and right side of antimeridian') + logging.info(f"Splitting bounds into left and right side of antimeridian") bounds_left, bounds_right = split_s1_bounds_at_am_crossing(bounds, lat_buff=0) - logging.info(f'Bounds left: {bounds_left}') - logging.info(f'Bounds right: {bounds_right}') + logging.info(f"Bounds left: {bounds_left}") + logging.info(f"Bounds right: {bounds_right}") # use recursion to create dems for the left and right side of AM # when passed back into the top function, this section will be skipped, creating # A valid dem for each side which we can then merge at the desired CRS # Add an additional buffer to ensure full coverage over dateline - left_save_path = '.'.join(str(save_path).split('.')[0:-1]) + "_left." + str(save_path).split('.')[-1] - logging.info(f'Getting tiles for left bounds') + left_save_path = ( + ".".join(str(save_path).split(".")[0:-1]) + + "_left." + + str(save_path).split(".")[-1] + ) + logging.info(f"Getting tiles for left bounds") get_cop30_dem_for_bounds( - bounds_left, - left_save_path, - ellipsoid_heights, + bounds_left, + left_save_path, + ellipsoid_heights, buffer_pixels=10, cop30_index_path=cop30_index_path, cop30_folder_path=cop30_folder_path, - geoid_tif_path=geoid_tif_path) - right_save_path = '.'.join(str(save_path).split('.')[0:-1]) + "_right." + str(save_path).split('.')[-1] - logging.info(f'Getting tiles for right bounds') + geoid_tif_path=geoid_tif_path, + ) + right_save_path = ( + ".".join(str(save_path).split(".")[0:-1]) + + "_right." + + str(save_path).split(".")[-1] + ) + logging.info(f"Getting tiles for right bounds") get_cop30_dem_for_bounds( - bounds_right, - right_save_path, - ellipsoid_heights, + bounds_right, + right_save_path, + ellipsoid_heights, buffer_pixels=10, cop30_index_path=cop30_index_path, cop30_folder_path=cop30_folder_path, - geoid_tif_path=geoid_tif_path) + geoid_tif_path=geoid_tif_path, + ) # reproject to 3031 and merge - logging.info(f'Reprojecting left and right side of antimeridian to EPGS:{target_crs}') + logging.info( + f"Reprojecting left and right side of antimeridian to EPGS:{target_crs}" + ) reproject_raster(left_save_path, left_save_path, target_crs) reproject_raster(right_save_path, right_save_path, target_crs) - logging.info(f'Merging across antimeridian') - dem_arr, dem_profile = merge_raster_files([left_save_path, right_save_path], output_path=save_path) + logging.info(f"Merging across antimeridian") + dem_arr, dem_profile = merge_raster_files( + [left_save_path, right_save_path], output_path=save_path + ) os.remove(left_save_path) os.remove(right_save_path) return dem_arr, dem_profile else: - logging.info(f'Getting cop30m dem for bounds: {bounds}') + logging.info(f"Getting cop30m dem for bounds: {bounds}") if adjust_for_high_lat_and_buffer: - logging.info(f'Expanding bounds by buffer and for high latitude warping') + logging.info(f"Expanding bounds by buffer and for high latitude warping") bounds = expand_bounds(bounds, buffer=0.1) - logging.info(f'Getting cop30m dem for expanded bounds: {bounds}') + logging.info(f"Getting cop30m dem for expanded bounds: {bounds}") if cop30_index_path: - logging.info(f'Finding intersecting DEM files from: {cop30_index_path}') - dem_paths = find_required_dem_paths_from_index(bounds, cop30_index_path=cop30_index_path) + logging.info(f"Finding intersecting DEM files from: {cop30_index_path}") + dem_paths = find_required_dem_paths_from_index( + bounds, cop30_index_path=cop30_index_path + ) else: - logging.info(f'Searching for DEM in folder: {cop30_folder_path}') - dem_paths = find_required_dem_tile_paths_by_filename(bounds,cop30_folder_path=cop30_folder_path) - logging.info(f'{len(dem_paths)} tiles found in bounds') + logging.info(f"Searching for DEM in folder: {cop30_folder_path}") + dem_paths = find_required_dem_tile_paths_by_filename( + bounds, cop30_folder_path=cop30_folder_path + ) + logging.info(f"{len(dem_paths)} tiles found in bounds") if len(dem_paths) == 0: - logging.warning('No DEM tiles found, assuming over water and creating zero dem for bounds') + logging.warning( + "No DEM tiles found, assuming over water and creating zero dem for bounds" + ) fill_value = 0 dem_profile = make_empty_cop30m_profile(bounds) dem_arr, dem_profile = expand_raster_to_bounds( - bounds, - src_profile=dem_profile, + bounds, + src_profile=dem_profile, save_path=save_path, fill_value=fill_value, - buffer_pixels=buffer_pixels) + buffer_pixels=buffer_pixels, + ) else: - logging.info(f'Merging tiles and reading data') + logging.info(f"Merging tiles and reading data") dem_arr, dem_profile = merge_raster_files( - dem_paths, - output_path=save_path, - bounds=bounds, - buffer_pixels=buffer_pixels + dem_paths, + output_path=save_path, + bounds=bounds, + buffer_pixels=buffer_pixels, ) - logging.info(f'Check the dem covers the required bounds') + logging.info(f"Check the dem covers the required bounds") dem_bounds = bounds_from_profile(dem_profile) - logging.info(f'Dem bounds: {dem_bounds}') - logging.info(f'Target bounds: {bounds}') + logging.info(f"Dem bounds: {dem_bounds}") + logging.info(f"Target bounds: {bounds}") bounds_filled_by_dem = box(*bounds).within(box(*dem_bounds)) - logging.info(f'Dem covers target: {bounds_filled_by_dem}') + logging.info(f"Dem covers target: {bounds_filled_by_dem}") if not bounds_filled_by_dem: - warn_msg = "The Cop30 DEM bounds do not fully cover the requested bounds. " \ - "Try increasing the 'buffer_pixels' value. Note at the antimeridian " \ - "This is expected, with bounds being slighly smaller on +ve side. " \ + warn_msg = ( + "The Cop30 DEM bounds do not fully cover the requested bounds. " + "Try increasing the 'buffer_pixels' value. Note at the antimeridian " + "This is expected, with bounds being slighly smaller on +ve side. " "e.g. max_lon is 179.9999 < 180." + ) logging.warning(warn_msg) if ellipsoid_heights: - logging.info(f'Subtracting the geoid from the DEM to return ellipsoid heights') - logging.info(f'Using geoid file: {geoid_tif_path}') + logging.info( + f"Subtracting the geoid from the DEM to return ellipsoid heights" + ) + logging.info(f"Using geoid file: {geoid_tif_path}") dem_arr = remove_geoid( - dem_arr = dem_arr, - dem_profile = dem_profile, - geoid_path = geoid_tif_path, - dem_area_or_point = 'Point', - buffer_pixels = 2, + dem_arr=dem_arr, + dem_profile=dem_profile, + geoid_path=geoid_tif_path, + dem_area_or_point="Point", + buffer_pixels=2, save_path=save_path, ) return dem_arr, dem_profile + def expand_bounds(bounds: tuple, buffer: float) -> tuple: """Expand the bounds for high lattitudes, and add a buffer. The provided bounds sometimes do not contain the full scene due to @@ -200,18 +236,21 @@ def expand_bounds(bounds: tuple, buffer: float) -> tuple: # adjust the bounds at high norther latitudes bounds = adjust_bounds(bounds, src_crs=4326, ref_crs=3995) exp_bounds = list(box(*bounds).buffer(buffer).bounds) - exp_bounds[0] = bounds[0] if exp_bounds[0] < -180 else exp_bounds[0] # keep original - exp_bounds[2] = bounds[2] if exp_bounds[2] > 180 else exp_bounds[2] # keep original + exp_bounds[0] = ( + bounds[0] if exp_bounds[0] < -180 else exp_bounds[0] + ) # keep original + exp_bounds[2] = bounds[2] if exp_bounds[2] > 180 else exp_bounds[2] # keep original return tuple(exp_bounds) + def find_required_dem_tile_paths_by_filename( - bounds: tuple, - check_exists : bool = True, - cop30_folder_path: Path = COP30_FOLDER_PATH, - search_buffer = 0.5, - tifs_in_subfolder = True, - )->list[str]: - """generate a list of the required dem paths based on the bounding coords. The + bounds: tuple, + check_exists: bool = True, + cop30_folder_path: Path = COP30_FOLDER_PATH, + search_buffer=0.5, + tifs_in_subfolder=True, +) -> list[str]: + """generate a list of the required dem paths based on the bounding coords. The function searches the specified folder. Parameters @@ -234,13 +273,13 @@ def find_required_dem_tile_paths_by_filename( # logic to find the correct files based on data being stored in each tile folder min_lat = np.floor(bounds[1]) if bounds[1] < 0 else np.ceil(bounds[1]) - max_lat = np.ceil(bounds[3]) if bounds[3] < 0 else np.floor(bounds[3])+1 + max_lat = np.ceil(bounds[3]) if bounds[3] < 0 else np.floor(bounds[3]) + 1 min_lon = np.floor(bounds[0]) if bounds[0] < 0 else np.floor(bounds[0]) max_lon = np.ceil(bounds[2]) if bounds[2] < 0 else np.ceil(bounds[2]) lat_range = list(range(int(min_lat), int(max_lat))) lon_range = list(range(int(min_lon), int(max_lon))) - logging.info(f'lat tile range: {lat_range}') - logging.info(f'lon tile range: {lon_range}') + logging.info(f"lat tile range: {lat_range}") + logging.info(f"lon tile range: {lon_range}") dem_paths = [] dem_folders = [] @@ -266,18 +305,21 @@ def find_required_dem_tile_paths_by_filename( logging.info(p) return sorted(list(set(dem_paths))) + def find_required_dem_paths_from_index( - bounds: tuple, - cop30_index_path = COP30_GPKG_PATH, - search_buffer=0.5, - ) -> list[str]: - + bounds: tuple, + cop30_index_path=COP30_GPKG_PATH, + search_buffer=0.5, +) -> list[str]: + gdf = gpd.read_file(cop30_index_path) bounding_box = box(*bounds).buffer(search_buffer) - + if gdf.crs is not None: # ensure same crs - bounding_box = gpd.GeoSeries([bounding_box], crs="EPSG:4326").to_crs(gdf.crs).iloc[0] + bounding_box = ( + gpd.GeoSeries([bounding_box], crs="EPSG:4326").to_crs(gdf.crs).iloc[0] + ) # Find rows that intersect with the bounding box intersecting_tiles = gdf[gdf.intersects(bounding_box)] if len(intersecting_tiles) > 0: @@ -285,10 +327,13 @@ def find_required_dem_paths_from_index( else: return [] -def check_s1_bounds_cross_antimeridian(bounds : tuple, max_scene_width : int =20) -> bool: - """Check if the s1 scene bounds cross the antimeridian. The bounds of a sentinel-1 + +def check_s1_bounds_cross_antimeridian( + bounds: tuple, max_scene_width: int = 20 +) -> bool: + """Check if the s1 scene bounds cross the antimeridian. The bounds of a sentinel-1 are valid at the antimeridian, just very large. By setting a max scene width, we - can determine if the antimeridian is crossed. Alternate scenario is a bounds + can determine if the antimeridian is crossed. Alternate scenario is a bounds with a very large width (i.e. close to the width of the earth). Parameters @@ -304,14 +349,15 @@ def check_s1_bounds_cross_antimeridian(bounds : tuple, max_scene_width : int =20 if the bounds cross the antimeridian """ - min_x = -180 + max_scene_width # -160 - max_x = 180 - max_scene_width # 160 + min_x = -180 + max_scene_width # -160 + max_x = 180 - max_scene_width # 160 if (bounds[0] < min_x) and (bounds[0] > -180): if bounds[2] > max_x and bounds[2] < 180: return True return False - -def split_s1_bounds_at_am_crossing(bounds: tuple, lat_buff : float = 0) -> list[tuple]: + + +def split_s1_bounds_at_am_crossing(bounds: tuple, lat_buff: float = 0) -> list[tuple]: """split the s1 bounds into bounds on the left and right of the antimeridian. @@ -326,17 +372,18 @@ def split_s1_bounds_at_am_crossing(bounds: tuple, lat_buff : float = 0) -> list[ ------- list[tuple] a list containing two sets of bounds for the left and right of the antimeridian - """ - max_negative_x = min([x for x in [bounds[0],bounds[2]] if x < 0]) - min_positive_x = min([x for x in [bounds[0],bounds[2]] if x > 0]) - min_y = min([bounds[1],bounds[3]]) - lat_buff - max_y = max([bounds[1],bounds[3]]) + lat_buff + """ + max_negative_x = min([x for x in [bounds[0], bounds[2]] if x < 0]) + min_positive_x = min([x for x in [bounds[0], bounds[2]] if x > 0]) + min_y = min([bounds[1], bounds[3]]) - lat_buff + max_y = max([bounds[1], bounds[3]]) + lat_buff min_y = -90 if min_y < -90 else min_y max_y = 90 if max_y > 90 else max_y bounds_left = (-180, min_y, max_negative_x, max_y) bounds_right = (min_positive_x, min_y, 180, max_y) return [tuple(bounds_left), tuple(bounds_right)] + def get_target_antimeridian_projection(bounds: tuple) -> int: """depending where were are on the earth, the desired crs at the antimeridian will change. e.g. polar stereo @@ -354,9 +401,14 @@ def get_target_antimeridian_projection(bounds: tuple) -> int: The CRS in integer form: e.g. 3031 """ min_lat = min(bounds[1], bounds[3]) - target_crs = 3031 if min_lat < -50 else 3995 if min_lat > 50 else get_local_utm(bounds, antimeridian=True) + target_crs = ( + 3031 + if min_lat < -50 + else 3995 if min_lat > 50 else get_local_utm(bounds, antimeridian=True) + ) return target_crs + def make_empty_cop30m_profile(bounds: tuple) -> dict: """make an empty cop30m dem rasterio profile based on a set of bounds. The desired pixel spacing changes based on lattitude @@ -380,35 +432,35 @@ def make_empty_cop30m_profile(bounds: tuple) -> dict: """ lat_res = 0.0002777777777777778 - mean_lat = abs((bounds[1] + bounds[3])/2) + mean_lat = abs((bounds[1] + bounds[3]) / 2) if mean_lat < 50: lon_res = lat_res elif mean_lat < 60: - lon_res = lat_res*1.5 + lon_res = lat_res * 1.5 elif mean_lat < 70: - lon_res = lat_res*2 + lon_res = lat_res * 2 elif mean_lat < 80: - lon_res = lat_res*3 + lon_res = lat_res * 3 elif mean_lat < 85: - lon_res = lat_res*5 + lon_res = lat_res * 5 elif mean_lat < 90: - lon_res = lat_res*10 + lon_res = lat_res * 10 else: - raise ValueError('cannot resolve cop30m lattitude') + raise ValueError("cannot resolve cop30m lattitude") min_x, min_y, max_x, max_y = bounds transform = Affine.translation(min_x, max_y) * Affine.scale(lon_res, -lat_res) return { - 'driver': 'GTiff', - 'dtype': 'float32', - 'nodata': np.nan, - 'width': abs(int((bounds[2] - bounds[0]) / lon_res)), - 'height': abs(int((bounds[3] - bounds[1]) / lat_res)), - 'count': 1, - 'crs': CRS.from_epsg(4326), - 'transform': transform, - 'blockysize': 1, - 'tiled': False, - 'interleave': 'band' - } \ No newline at end of file + "driver": "GTiff", + "dtype": "float32", + "nodata": np.nan, + "width": abs(int((bounds[2] - bounds[0]) / lon_res)), + "height": abs(int((bounds[3] - bounds[1]) / lat_res)), + "count": 1, + "crs": CRS.from_epsg(4326), + "transform": transform, + "blockysize": 1, + "tiled": False, + "interleave": "band", + } diff --git a/sar_antarctica/nci/preparation/geoid.py b/sar_antarctica/nci/preparation/geoid.py index 2954a8c..cca3c01 100644 --- a/sar_antarctica/nci/preparation/geoid.py +++ b/sar_antarctica/nci/preparation/geoid.py @@ -1,6 +1,7 @@ """ inspired by https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher/blob/dev/src/dem_stitcher/geoid.py """ + import os from typing import Union from pathlib import Path @@ -13,7 +14,10 @@ from ...utils.raster import read_raster_with_bounds from ...utils.rio_tools import translate_profile, reproject_arr_to_match_profile -def read_geoid(geoid_path: Union[str, Path], bounds: tuple, buffer_pixels: int = 0) -> tuple[np.ndarray, dict]: + +def read_geoid( + geoid_path: Union[str, Path], bounds: tuple, buffer_pixels: int = 0 +) -> tuple[np.ndarray, dict]: """Read in the geoid for the bounds provided with a specified buffer. Parameters @@ -35,14 +39,16 @@ def read_geoid(geoid_path: Union[str, Path], bounds: tuple, buffer_pixels: int = FileNotFoundError If ther GEOID file cannot be found """ - + if not os.path.exists(geoid_path): - raise FileNotFoundError(f'Geoid file does not exist at path: {geoid_path}') + raise FileNotFoundError(f"Geoid file does not exist at path: {geoid_path}") - geoid_arr, geoid_profile = read_raster_with_bounds(geoid_path, bounds, buffer_pixels=buffer_pixels) - geoid_arr = geoid_arr.astype('float32') - geoid_arr[geoid_profile['nodata'] == geoid_arr] = np.nan - geoid_profile['nodata'] = np.nan + geoid_arr, geoid_profile = read_raster_with_bounds( + geoid_path, bounds, buffer_pixels=buffer_pixels + ) + geoid_arr = geoid_arr.astype("float32") + geoid_arr[geoid_profile["nodata"] == geoid_arr] = np.nan + geoid_profile["nodata"] = np.nan return geoid_arr, geoid_profile @@ -51,11 +57,11 @@ def remove_geoid( dem_arr: np.ndarray, dem_profile: dict, geoid_path: Union[str, Path], - dem_area_or_point: str = 'Point', + dem_area_or_point: str = "Point", buffer_pixels: int = 2, - save_path: Union[str, Path] = '', + save_path: Union[str, Path] = "", ) -> np.ndarray: - """Subtract the Geoid from a dem file. Result will be + """Subtract the Geoid from a dem file. Result will be ellipsoid referenced heights for the cop30m dem. Parameters @@ -79,38 +85,44 @@ def remove_geoid( np.ndarray original array with the geoid subtracted """ - - assert dem_area_or_point in ['Point', 'Area'] - bounds = array_bounds(dem_profile['height'], dem_profile['width'], dem_profile['transform']) + assert dem_area_or_point in ["Point", "Area"] + + bounds = array_bounds( + dem_profile["height"], dem_profile["width"], dem_profile["transform"] + ) - geoid_arr, geoid_profile = read_geoid(geoid_path, bounds=tuple(bounds), buffer_pixels=buffer_pixels) + geoid_arr, geoid_profile = read_geoid( + geoid_path, bounds=tuple(bounds), buffer_pixels=buffer_pixels + ) - t_dem = dem_profile['transform'] - t_geoid = geoid_profile['transform'] + t_dem = dem_profile["transform"] + t_geoid = geoid_profile["transform"] res_dem = max(t_dem.a, abs(t_dem.e)) res_geoid = max(t_geoid.a, abs(t_geoid.e)) if res_geoid * buffer_pixels <= res_dem: buffer_recommendation = int(np.ceil(res_dem / res_geoid)) warning = ( - 'The dem resolution is larger than the geoid resolution and its buffer; ' - 'Edges resampled with bilinear interpolation will be inconsistent so select larger buffer.' - f'Select a `buffer_pixels = {buffer_recommendation}`' + "The dem resolution is larger than the geoid resolution and its buffer; " + "Edges resampled with bilinear interpolation will be inconsistent so select larger buffer." + f"Select a `buffer_pixels = {buffer_recommendation}`" ) logging.warning(warning) # Translate geoid if necessary as all geoids have Area tag - if dem_area_or_point == 'Point': + if dem_area_or_point == "Point": shift = -0.5 geoid_profile = translate_profile(geoid_profile, shift, shift) - geoid_offset, _ = reproject_arr_to_match_profile(geoid_arr, geoid_profile, dem_profile, resampling='bilinear') + geoid_offset, _ = reproject_arr_to_match_profile( + geoid_arr, geoid_profile, dem_profile, resampling="bilinear" + ) dem_arr_offset = dem_arr + geoid_offset if save_path: - with rasterio.open(save_path, 'w', **dem_profile) as dst: + with rasterio.open(save_path, "w", **dem_profile) as dst: dst.write(dem_arr_offset) - return dem_arr_offset \ No newline at end of file + return dem_arr_offset diff --git a/sar_antarctica/nci/preparation/orbits.py b/sar_antarctica/nci/preparation/orbits.py index 69f7c43..afee3a9 100644 --- a/sar_antarctica/nci/preparation/orbits.py +++ b/sar_antarctica/nci/preparation/orbits.py @@ -3,13 +3,17 @@ import re from typing import Optional -from sar_antarctica.nci.preparation.scenes import parse_scene_file_dates, parse_scene_file_sensor +from sar_antarctica.nci.preparation.scenes import ( + parse_scene_file_dates, + parse_scene_file_sensor, +) # Constants for NCI S1_DIR = Path("/g/data/fj7/Copernicus/Sentinel-1/") POE_DIR = "POEORB" RES_DIR = "RESORB" + def parse_orbit_file_dates(orbit_file_name: str) -> tuple[datetime, datetime, datetime]: """Extracts published_date, start_date, and end_date from the given orbit file. Filename example: S1A_OPER_AUX_POEORB_OPOD_20141207T123431_V20141115T225944_20141117T005944.EOF @@ -33,9 +37,11 @@ def parse_orbit_file_dates(orbit_file_name: str) -> tuple[datetime, datetime, da Did not find a match to the expected date pattern of published_date followed by start_date and end_date """ # Regex pattern to match the dates - pattern = (r"(?P\d{8}T\d{6})_V" - r"(?P\d{8}T\d{6})_" - r"(?P\d{8}T\d{6})\.EOF") + pattern = ( + r"(?P\d{8}T\d{6})_V" + r"(?P\d{8}T\d{6})_" + r"(?P\d{8}T\d{6})\.EOF" + ) # Search for matches in the file name match = re.search(pattern, str(orbit_file_name)) @@ -44,14 +50,17 @@ def parse_orbit_file_dates(orbit_file_name: str) -> tuple[datetime, datetime, da raise ValueError("The input string does not match the expected format.") # Extract and parse the dates into datetime objects - published_date = datetime.strptime(match.group('published_date'), "%Y%m%dT%H%M%S") - start_date = datetime.strptime(match.group('start_date'), "%Y%m%dT%H%M%S") - stop_date = datetime.strptime(match.group('stop_date'), "%Y%m%dT%H%M%S") + published_date = datetime.strptime(match.group("published_date"), "%Y%m%dT%H%M%S") + start_date = datetime.strptime(match.group("start_date"), "%Y%m%dT%H%M%S") + stop_date = datetime.strptime(match.group("stop_date"), "%Y%m%dT%H%M%S") return (published_date, start_date, stop_date) -def find_latest_orbit_for_scene(scene_id: str, orbit_type: Optional[str] = None) -> Path: - """Identifies the most recent orbit file available for a given scene, based + +def find_latest_orbit_for_scene( + scene_id: str, orbit_type: Optional[str] = None +) -> Path: + """Identifies the most recent orbit file available for a given scene, based on the scene's start and end date. Parameters @@ -97,9 +106,11 @@ def find_latest_orbit_for_scene(scene_id: str, orbit_type: Optional[str] = None) for orbit_file in orbit_files: - orbit_published, orbit_start, orbit_stop = parse_orbit_file_dates(orbit_file) - - # Check if scene falls within orbit + orbit_published, orbit_start, orbit_stop = parse_orbit_file_dates( + orbit_file + ) + + # Check if scene falls within orbit if scene_start >= orbit_start and scene_stop <= orbit_stop: orbit_metadata = (orbit_file, orbit_dir, orbit_published) relevant_orbits.append(orbit_metadata) @@ -111,6 +122,5 @@ def find_latest_orbit_for_scene(scene_id: str, orbit_type: Optional[str] = None) raise ValueError("No valid orbit was found.") else: latest_orbit_file = latest_orbit[0] - - return latest_orbit_file + return latest_orbit_file diff --git a/sar_antarctica/nci/preparation/scenes.py b/sar_antarctica/nci/preparation/scenes.py index b629232..f7f4b35 100644 --- a/sar_antarctica/nci/preparation/scenes.py +++ b/sar_antarctica/nci/preparation/scenes.py @@ -4,6 +4,7 @@ SCENE_DIR = Path("/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/") + def parse_scene_file_sensor(scene_id: str) -> str: """Extract Sentinel-1 sensor string (SA1,S1B,S1C,S1D) from scene ID @@ -18,18 +19,20 @@ def parse_scene_file_sensor(scene_id: str) -> str: str Sensor string. Should be one of S1A, S1B, S1C, or S1D - Raises + Raises ------ ValueError Did not find any of S1A, S1B, S1C, or S1D in the scene ID """ # Expect files to be prefaced with any of S1A, S1B, S1C, or S1D, followed by underscore - pattern=r"^(S1[A|B|C|D])_" + pattern = r"^(S1[A|B|C|D])_" match = re.match(pattern, scene_id) if not match: - raise ValueError("No valid sensor was found in the scene ID. Valid sensors are S1A, S1B, S1C, or S1D") + raise ValueError( + "No valid sensor was found in the scene ID. Valid sensors are S1A, S1B, S1C, or S1D" + ) return match.group(1) @@ -55,19 +58,19 @@ def parse_scene_file_dates(scene_id: str) -> tuple[datetime, datetime]: Did not find a match to the expected date pattern of start_date followed by end_date in the scene ID """ # Regex pattern to match the dates - pattern = (r"(?P\d{8}T\d{6})_" - r"(?P\d{8}T\d{6})_") - + pattern = r"(?P\d{8}T\d{6})_" r"(?P\d{8}T\d{6})_" + match = re.search(pattern, scene_id) if not match: raise ValueError("The input string does not match the expected format.") - - start_date = datetime.strptime(match.group('start_date'), "%Y%m%dT%H%M%S") - stop_date = datetime.strptime(match.group('stop_date'), "%Y%m%dT%H%M%S") + + start_date = datetime.strptime(match.group("start_date"), "%Y%m%dT%H%M%S") + stop_date = datetime.strptime(match.group("stop_date"), "%Y%m%dT%H%M%S") return (start_date, stop_date) + def find_scene_file_from_id(scene_id: str) -> Path: """Finds the path to the scene on GADI based on the scene ID @@ -91,11 +94,11 @@ def find_scene_file_from_id(scene_id: str) -> Path: """ # Parse the scene dates -- only start date is needed for search - scene_start, _ = parse_scene_file_dates(scene_id) + scene_start, _ = parse_scene_file_dates(scene_id) # Extract year and month of first path to provide for file search - year = scene_start.strftime('%Y') - month = scene_start.strftime('%m') + year = scene_start.strftime("%Y") + month = scene_start.strftime("%m") # Set path on GADI and search search_path = SCENE_DIR.joinpath(f"{year}/{year}-{month}/") @@ -107,6 +110,8 @@ def find_scene_file_from_id(scene_id: str) -> Path: elif len(file_path) > 1: raise RuntimeError("More than one file found. Review before proceeding") else: - raise RuntimeError("No files found or some other error. Review before proceeding") - - return scene_path \ No newline at end of file + raise RuntimeError( + "No files found or some other error. Review before proceeding" + ) + + return scene_path diff --git a/sar_antarctica/nci/processing/GAMMA/GAMMA_utils.py b/sar_antarctica/nci/processing/GAMMA/GAMMA_utils.py index 6ac3191..06e0687 100644 --- a/sar_antarctica/nci/processing/GAMMA/GAMMA_utils.py +++ b/sar_antarctica/nci/processing/GAMMA/GAMMA_utils.py @@ -3,6 +3,7 @@ GAMMA_HOME_ENV = "GAMMA_HOME" LD_LIBRARY_ENV = "LD_LIBRARY_PATH" + def set_gamma_env_variables(gamma_home_path: str, ld_libraries_path: str): gamma_env = os.environ.get(GAMMA_HOME_ENV, None) ld_lib_env = os.environ.get(LD_LIBRARY_ENV, None) @@ -13,4 +14,4 @@ def set_gamma_env_variables(gamma_home_path: str, ld_libraries_path: str): if ld_lib_env is None: os.environ[LD_LIBRARY_ENV] = ld_libraries_path else: - os.environ[LD_LIBRARY_ENV] = os.path.join(ld_libraries_path, ":", ld_lib_env) \ No newline at end of file + os.environ[LD_LIBRARY_ENV] = os.path.join(ld_libraries_path, ":", ld_lib_env) diff --git a/sar_antarctica/nci/processing/pyroSAR/pyrosar_geocode.py b/sar_antarctica/nci/processing/pyroSAR/pyrosar_geocode.py index b3d1036..31052ec 100644 --- a/sar_antarctica/nci/processing/pyroSAR/pyrosar_geocode.py +++ b/sar_antarctica/nci/processing/pyroSAR/pyrosar_geocode.py @@ -18,17 +18,18 @@ log = logging.getLogger("gammapy") log.setLevel(logging.INFO) + @click.command() @click.argument("workflow_config", nargs=1) @click.argument("scene_config", nargs=1) def cli(workflow_config: str, scene_config: str): - + # Read in config file with open(workflow_config, "rb") as f: - workflow_config_dict = tomli.load(f) + workflow_config_dict = tomli.load(f) with open(scene_config, "rb") as f: - scene_config_dict = tomli.load(f) + scene_config_dict = tomli.load(f) # Split config dicts up to ease readability config_inputs = scene_config_dict["inputs"] @@ -38,8 +39,7 @@ def cli(workflow_config: str, scene_config: str): # Environment variables for GAMMA must be set set_gamma_env_variables( - config_gamma["software_env_var"], - config_gamma["libs_env_var"] + config_gamma["software_env_var"], config_gamma["libs_env_var"] ) # Identify scene @@ -51,8 +51,14 @@ def cli(workflow_config: str, scene_config: str): # Construct output scenes data_dir = Path(config_outputs["data"]) - processed_scene_dir = data_dir / config_outputs["processed"] / pyrosar_scene_id.outname_base(extensions=None) - pyrosar_temp_dir = data_dir / "temp" / pyrosar_scene_id.outname_base(extensions=None) + processed_scene_dir = ( + data_dir + / config_outputs["processed"] + / pyrosar_scene_id.outname_base(extensions=None) + ) + pyrosar_temp_dir = ( + data_dir / "temp" / pyrosar_scene_id.outname_base(extensions=None) + ) pyrosar_temp_log_dir = pyrosar_temp_dir / "logfiles" log.info("creating directories:") @@ -76,11 +82,11 @@ def cli(workflow_config: str, scene_config: str): log.info("running DEM") dem_import( - src=str(dem_tif), - dst=str(dem_gamma), - geoid=None, - logpath=str(pyrosar_temp_log_dir), - outdir=str(dem_tif.parent) + src=str(dem_tif), + dst=str(dem_gamma), + geoid=None, + logpath=str(pyrosar_temp_log_dir), + outdir=str(dem_tif.parent), ) log.info("finished DEM") @@ -90,29 +96,36 @@ def cli(workflow_config: str, scene_config: str): log.info("running geocode") geocode( - scene=pyrosar_scene_id, - dem=str(dem_gamma), + scene=pyrosar_scene_id, + dem=str(dem_gamma), tmpdir=str(pyrosar_temp_dir), - outdir=str(processed_scene_dir), - spacing=config_geocode["spacing"], - scaling=config_geocode["scaling"], + outdir=str(processed_scene_dir), + spacing=config_geocode["spacing"], + scaling=config_geocode["scaling"], func_geoback=1, - nodata=(0, -99), - update_osv=False, - osvdir=str(pyrosar_temp_dir), + nodata=(0, -99), + update_osv=False, + osvdir=str(pyrosar_temp_dir), allow_RES_OSV=False, - cleanup=False, - export_extra=['inc_geo','dem_seg_geo','ls_map_geo','pix_area_gamma0_geo','pix_ratio_geo'], + cleanup=False, + export_extra=[ + "inc_geo", + "dem_seg_geo", + "ls_map_geo", + "pix_area_gamma0_geo", + "pix_ratio_geo", + ], basename_extensions=None, - removeS1BorderNoiseMethod='pyroSAR', - refine_lut=False, - rlks=None, + removeS1BorderNoiseMethod="pyroSAR", + refine_lut=False, + rlks=None, azlks=None, - s1_osv_url_option=1 + s1_osv_url_option=1, ) log.info("finished geocode") + if __name__ == "__main__": - cli() \ No newline at end of file + cli() diff --git a/sar_antarctica/nci/submission/pyrosar_gamma/pyrosar_gamma.py b/sar_antarctica/nci/submission/pyrosar_gamma/pyrosar_gamma.py index 6ea9272..015761a 100644 --- a/sar_antarctica/nci/submission/pyrosar_gamma/pyrosar_gamma.py +++ b/sar_antarctica/nci/submission/pyrosar_gamma/pyrosar_gamma.py @@ -7,6 +7,7 @@ WORKFLOW = "pyrosar_gamma" PROCESSING_DIR = "/g/data/yp75/projects/sar-antractica-processing" + def get_list_of_scenes(scene_source: str) -> list[str]: """Convert script input to list. If a .zip file, produce a list with that. @@ -28,17 +29,22 @@ def get_list_of_scenes(scene_source: str) -> list[str]: scene_list = [scene_source] # Process a .txt file containing .zip files elif scene_source.endswith(".txt"): - with open(scene_source, 'r') as f: - scene_list = [line.strip() for line in f if line.strip().endswith('.zip')] + with open(scene_source, "r") as f: + scene_list = [line.strip() for line in f if line.strip().endswith(".zip")] else: scene_list = [] if scene_list is not None: return scene_list else: - raise RuntimeError("No valid scenes were found for processing. Expected single .zip file or .txt file containing at least one .zip file.") + raise RuntimeError( + "No valid scenes were found for processing. Expected single .zip file or .txt file containing at least one .zip file." + ) + -def update_pbs_template(pbs_template: str, scene_id: str, job_config: dict[str, str | dict[str, Any]]) -> str: +def update_pbs_template( + pbs_template: str, scene_id: str, job_config: dict[str, str | dict[str, Any]] +) -> str: """_summary_ Parameters @@ -73,7 +79,11 @@ def update_pbs_template(pbs_template: str, scene_id: str, job_config: dict[str, """ processing_path = Path(job_config["root"]) - log_path = processing_path / job_config["submission"]["root"] / job_config["submission"]["logs"] + log_path = ( + processing_path + / job_config["submission"]["root"] + / job_config["submission"]["logs"] + ) config_path = processing_path / job_config["configuration"]["root"] job_configuration = job_config["configuration"] @@ -90,12 +100,16 @@ def update_pbs_template(pbs_template: str, scene_id: str, job_config: dict[str, "": job_settings["walltime"], "": job_settings["storage"], "": log_path, - "": config_path / job_configuration["workflow"] / f"{workflow_config}.toml", - "": config_path / job_configuration["scene"] / f"{scene_id}.toml" + "": config_path + / job_configuration["workflow"] + / f"{workflow_config}.toml", + "": config_path / job_configuration["scene"] / f"{scene_id}.toml", } for key, value in replace_dict.items(): - pbs_template = pbs_template.replace(key, value if isinstance(value, str) else str(value)) + pbs_template = pbs_template.replace( + key, value if isinstance(value, str) else str(value) + ) return pbs_template @@ -103,7 +117,9 @@ def update_pbs_template(pbs_template: str, scene_id: str, job_config: dict[str, @click.command() @click.argument("config_file", nargs=1) @click.argument("scene_source", nargs=1) -def pyrosar_gamma_workflow(config_file: str | os.PathLike, scene_source: str | os.PathLike) -> None: +def pyrosar_gamma_workflow( + config_file: str | os.PathLike, scene_source: str | os.PathLike +) -> None: """Take an input of a single scene or file with multiple scenes and submit pyroSAR+GAMMA jobs Parameters @@ -117,7 +133,7 @@ def pyrosar_gamma_workflow(config_file: str | os.PathLike, scene_source: str | o current_file_directory = Path(__file__).resolve().parent with open(config_file, "rb") as f: - config = tomli.load(f) + config = tomli.load(f) # Extract specific configuration dictionaries job_config = config["job"] @@ -139,15 +155,11 @@ def pyrosar_gamma_workflow(config_file: str | os.PathLike, scene_source: str | o scene_script.parent.mkdir(exist_ok=True, parents=True) # Read the workflow template and replace values - workflow_name = settings_config['workflow_config'] + workflow_name = settings_config["workflow_config"] template_file = current_file_directory / f"{workflow_name}.txt" print(template_file) pbs_template = template_file.read_text() - pbs_template = update_pbs_template( - pbs_template, - scene_id, - job_config - ) + pbs_template = update_pbs_template(pbs_template, scene_id, job_config) # Write updated text to pbs script scene_script.write_text(pbs_template) @@ -158,4 +170,4 @@ def pyrosar_gamma_workflow(config_file: str | os.PathLike, scene_source: str | o if __name__ == "__main__": - pyrosar_gamma_workflow() \ No newline at end of file + pyrosar_gamma_workflow() diff --git a/sar_antarctica/utils/raster.py b/sar_antarctica/utils/raster.py index 4684efe..d7e5054 100644 --- a/sar_antarctica/utils/raster.py +++ b/sar_antarctica/utils/raster.py @@ -16,9 +16,11 @@ from rasterio.crs import CRS from rasterio.windows import from_bounds + def bounds_from_profile(profile): # returns the bounds from a rasterio profile dict - return array_bounds(profile['height'], profile['width'], profile['transform']) + return array_bounds(profile["height"], profile["width"], profile["transform"]) + def reproject_raster(src_path: str, out_path: str, crs: int): """Reproject raster to desired crs @@ -40,19 +42,18 @@ def reproject_raster(src_path: str, out_path: str, crs: int): with rasterio.open(src_path) as src: src_crs = src.crs transform, width, height = calculate_default_transform( - src_crs, crs, src.width, src.height, *src.bounds) + src_crs, crs, src.width, src.height, *src.bounds + ) kwargs = src.meta.copy() - # get crs proj + # get crs proj crs = pyproj.CRS(f"EPSG:{crs}") - kwargs.update({ - 'crs': crs, - 'transform': transform, - 'width': width, - 'height': height}) + kwargs.update( + {"crs": crs, "transform": transform, "width": width, "height": height} + ) - with rasterio.open(out_path, 'w', **kwargs) as dst: + with rasterio.open(out_path, "w", **kwargs) as dst: for i in range(1, src.count + 1): reproject( source=rasterio.band(src, i), @@ -61,17 +62,19 @@ def reproject_raster(src_path: str, out_path: str, crs: int): src_crs=src.crs, dst_transform=transform, dst_crs=crs, - resampling=Resampling.nearest) + resampling=Resampling.nearest, + ) + def expand_raster_to_bounds( - trg_bounds : tuple, - src_path : str = '', - src_profile : dict = None, - src_array : np.ndarray = None, - fill_value : float = 0, - buffer_pixels : int = 0, - save_path : str = '' - ) -> tuple[np.ndarray, dict]: + trg_bounds: tuple, + src_path: str = "", + src_profile: dict = None, + src_array: np.ndarray = None, + fill_value: float = 0, + buffer_pixels: int = 0, + save_path: str = "", +) -> tuple[np.ndarray, dict]: """Expand the extent of the input array to the target bounds specified by the user. Either a src_path to expand, src_profile to construct a new raster, or src_profile and src_array o expand must be provided/ @@ -99,8 +102,9 @@ def expand_raster_to_bounds( new array and rasterio profile """ - assert src_path or (src_profile and src_array is not None) or src_profile, \ - "Either src_path, src_array and src_profile, or src_profile must be provided." + assert ( + src_path or (src_profile and src_array is not None) or src_profile + ), "Either src_path, src_array and src_profile, or src_profile must be provided." if src_path: with rasterio.open(src_path) as src: @@ -108,25 +112,27 @@ def expand_raster_to_bounds( src_profile = src.profile src_left, src_bottom, src_right, src_top = src.bounds else: - src_bounds = array_bounds(src_profile['height'], src_profile['width'], src_profile['transform']) + src_bounds = array_bounds( + src_profile["height"], src_profile["width"], src_profile["transform"] + ) src_left, src_bottom, src_right, src_top = src_bounds - + # Define the new bounds trg_left, trg_bottom, trg_right, trg_top = trg_bounds - lon_res = abs(src_profile['transform'].a) # Pixel width - lat_res = abs(src_profile['transform'].e) # Pixel height + lon_res = abs(src_profile["transform"].a) # Pixel width + lat_res = abs(src_profile["transform"].e) # Pixel height # determine the number of new pixels in each direction - new_left_pixels = int(abs(trg_left-src_left)/lon_res) + buffer_pixels - new_right_pixels = int(abs(trg_right-src_right)/lon_res) + buffer_pixels - new_bottom_pixels = int(abs(trg_bottom-src_bottom)/lat_res) + buffer_pixels - new_top_pixels = int(abs(trg_top-src_top)/lat_res) + buffer_pixels - + new_left_pixels = int(abs(trg_left - src_left) / lon_res) + buffer_pixels + new_right_pixels = int(abs(trg_right - src_right) / lon_res) + buffer_pixels + new_bottom_pixels = int(abs(trg_bottom - src_bottom) / lat_res) + buffer_pixels + new_top_pixels = int(abs(trg_top - src_top) / lat_res) + buffer_pixels + # adjust the new bounds with even pixel multiples of existing - new_trg_left = src_left - new_left_pixels*lon_res - new_trg_right = src_right + new_right_pixels*lon_res - new_trg_bottom = src_bottom - new_bottom_pixels*lat_res - new_trg_top = src_top + new_top_pixels*lat_res + new_trg_left = src_left - new_left_pixels * lon_res + new_trg_right = src_right + new_right_pixels * lon_res + new_trg_bottom = src_bottom - new_bottom_pixels * lat_res + new_trg_top = src_top + new_top_pixels * lat_res # keep source if they are already greater than the desired bounds new_trg_left = src_left if src_left < trg_left else new_trg_left @@ -140,44 +146,45 @@ def expand_raster_to_bounds( # Define the new transformation matrix transform = from_origin(new_trg_left, new_trg_top, lon_res, lat_res) - + # Create a new raster dataset with expanded bounds fill_profile = src_profile.copy() - fill_profile.update({ - 'width': new_width, - 'height': new_height, - 'transform': transform - }) - fill_array = np.full((1, new_height, new_width), fill_value=fill_value, dtype=src_profile['dtype']) - + fill_profile.update( + {"width": new_width, "height": new_height, "transform": transform} + ) + fill_array = np.full( + (1, new_height, new_width), fill_value=fill_value, dtype=src_profile["dtype"] + ) + if src_array is not None: # if an existing src array (e.g. dem) is provided to expand trg_array, trg_profile = merge_arrays_with_geometadata( - arrays = [src_array, fill_array], - profiles = [src_profile, fill_profile], - resampling='bilinear', - nodata = src_profile['nodata'], - dtype = src_profile['dtype'], - method='first', - ) + arrays=[src_array, fill_array], + profiles=[src_profile, fill_profile], + resampling="bilinear", + nodata=src_profile["nodata"], + dtype=src_profile["dtype"], + method="first", + ) else: # we are not expanding an existing array # return the fill array that has been constructed based on the src_profile trg_array, trg_profile = fill_array, fill_profile if save_path: - with rasterio.open(save_path, 'w', **trg_profile) as dst: + with rasterio.open(save_path, "w", **trg_profile) as dst: dst.write(trg_array) return trg_array, trg_profile def read_vrt_in_bounds( - vrt_path: str, - bounds: tuple, - output_path: str = '', - return_data: bool =True, - buffer_pixels: int=0, - set_nodata : float = None): + vrt_path: str, + bounds: tuple, + output_path: str = "", + return_data: bool = True, + buffer_pixels: int = 0, + set_nodata: float = None, +): """Read in data from a vrt file in the specified bounds Parameters @@ -202,10 +209,10 @@ def read_vrt_in_bounds( tuple[np.ndarray, dict] new array and rasterio profile """ - + # make upper end of the requested integer # ensures the bounds are covered with requested pixel buffer - buffer_pixels+=0.9 + buffer_pixels += 0.9 if bounds is None: # get all data in tiles @@ -219,10 +226,12 @@ def read_vrt_in_bounds( with rasterio.open(vrt_path) as src: # Define the profile for the GeoTIFF arr_profile = src.profile - arr_profile.update(driver="GTiff") # Ensure the driver is set to GeoTIFF + arr_profile.update( + driver="GTiff" + ) # Ensure the driver is set to GeoTIFF arr = src.read() if set_nodata is not None: - arr_profile['nodata'] = set_nodata + arr_profile["nodata"] = set_nodata return arr, arr_profile else: @@ -245,70 +254,81 @@ def read_vrt_in_bounds( min_x - buffer_x, min_y - buffer_y, max_x + buffer_x, - max_y + buffer_y + max_y + buffer_y, ) # Create a window for the bounding box xmin, ymin, xmax, ymax = buffered_bounds - window = from_bounds(xmin, ymin, xmax, ymax, transform=src_transform) #.round() + window = from_bounds( + xmin, ymin, xmax, ymax, transform=src_transform + ) # .round() # Read data for the specified window - data = src.read(1, window=window) # Read the first band; adjust if you need multiple bands + data = src.read( + 1, window=window + ) # Read the first band; adjust if you need multiple bands # data[~np.isfinite(data)] = 0 # Adjust the transform for the window window_transform = src.window_transform(window) arr_profile = src.profile.copy() - arr_profile['transform'] = window_transform - arr_profile['driver'] = 'GTiff' - arr_profile['count'] = 1 - arr_profile['height'] = data.shape[0] - arr_profile['width'] = data.shape[1] + arr_profile["transform"] = window_transform + arr_profile["driver"] = "GTiff" + arr_profile["count"] = 1 + arr_profile["height"] = data.shape[0] + arr_profile["width"] = data.shape[1] if set_nodata is not None: - arr_profile['nodata'] = set_nodata + arr_profile["nodata"] = set_nodata - # Save the extracted data to a new GeoTIFF + # Save the extracted data to a new GeoTIFF if output_path: - with rasterio.open(output_path, 'w', **arr_profile) as dst: - dst.write(data,1) + with rasterio.open(output_path, "w", **arr_profile) as dst: + dst.write(data, 1) if return_data: return data[np.newaxis, :, :], arr_profile -def merge_raster_files(paths, output_path, bounds=None, return_data=True, buffer_pixels=0, delete_vrt=True): +def merge_raster_files( + paths, output_path, bounds=None, return_data=True, buffer_pixels=0, delete_vrt=True +): # Create a virtual raster (in-memory description of the merged DEMs) vrt_path = str(output_path).replace(".tif", ".vrt") # Temporary VRT file path - gdal.BuildVRT(vrt_path, paths, resolution='highest') + gdal.BuildVRT(vrt_path, paths, resolution="highest") res = read_vrt_in_bounds( vrt_path=vrt_path, bounds=bounds, output_path=output_path, buffer_pixels=buffer_pixels, - return_data=return_data) + return_data=return_data, + ) if delete_vrt: os.remove(vrt_path) if return_data: arr, arr_profile = res return arr, arr_profile + def merge_arrays_with_geometadata( arrays: list[np.ndarray], profiles: list[dict], - resampling: str = 'bilinear', + resampling: str = "bilinear", nodata: Union[float, int] = np.nan, dtype: str = None, - method: str = 'first', + method: str = "first", ) -> tuple[np.ndarray, dict]: # https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher/blob/dev/src/dem_stitcher/merge.py n_dim = arrays[0].shape if len(n_dim) not in [2, 3]: - raise ValueError('Currently arrays must be in BIP format' 'i.e. channels x height x width or flat array') + raise ValueError( + "Currently arrays must be in BIP format" + "i.e. channels x height x width or flat array" + ) if len(set([len(arr.shape) for arr in arrays])) != 1: - raise ValueError('All arrays must have same number of dimensions i.e. 2 or 3') + raise ValueError("All arrays must have same number of dimensions i.e. 2 or 3") if len(n_dim) == 2: arrays_input = [arr[np.newaxis, ...] for arr in arrays] @@ -316,31 +336,36 @@ def merge_arrays_with_geometadata( arrays_input = arrays if (len(arrays)) != (len(profiles)): - raise ValueError('Length of arrays and profiles needs to be the same') + raise ValueError("Length of arrays and profiles needs to be the same") memfiles = [MemoryFile() for p in profiles] datasets = [mfile.open(**p) for (mfile, p) in zip(memfiles, profiles)] [ds.write(arr) for (ds, arr) in zip(datasets, arrays_input)] merged_arr, merged_trans = merge( - datasets, resampling=Resampling[resampling], method=method, nodata=nodata, dtype=dtype + datasets, + resampling=Resampling[resampling], + method=method, + nodata=nodata, + dtype=dtype, ) prof_merged = profiles[0].copy() - prof_merged['transform'] = merged_trans - prof_merged['count'] = merged_arr.shape[0] - prof_merged['height'] = merged_arr.shape[1] - prof_merged['width'] = merged_arr.shape[2] + prof_merged["transform"] = merged_trans + prof_merged["count"] = merged_arr.shape[0] + prof_merged["height"] = merged_arr.shape[1] + prof_merged["width"] = merged_arr.shape[2] if nodata is not None: - prof_merged['nodata'] = nodata + prof_merged["nodata"] = nodata if dtype is not None: - prof_merged['dtype'] = dtype + prof_merged["dtype"] = dtype [ds.close() for ds in datasets] [mfile.close() for mfile in memfiles] return merged_arr, prof_merged + def read_raster_with_bounds(file_path, bounds, buffer_pixels=0): """ Reads a specific region of a raster file defined by bounds and returns the data array and profile. @@ -368,7 +393,7 @@ def read_raster_with_bounds(file_path, bounds, buffer_pixels=0): min_x - buffer_x, min_y - buffer_y, max_x + buffer_x, - max_y + buffer_y + max_y + buffer_y, ) # Create a window from the buffered bounds @@ -382,10 +407,12 @@ def read_raster_with_bounds(file_path, bounds, buffer_pixels=0): # Adjust the profile for the window profile = src.profile.copy() - profile.update({ - "height": window.height, - "width": window.width, - "transform": src.window_transform(window) - }) + profile.update( + { + "height": window.height, + "width": window.width, + "transform": src.window_transform(window), + } + ) return data, profile diff --git a/sar_antarctica/utils/rio_tools.py b/sar_antarctica/utils/rio_tools.py index d0264f9..832b58f 100644 --- a/sar_antarctica/utils/rio_tools.py +++ b/sar_antarctica/utils/rio_tools.py @@ -9,7 +9,12 @@ from rasterio import DatasetReader from rasterio.crs import CRS from rasterio.io import MemoryFile -from rasterio.warp import Resampling, aligned_target, calculate_default_transform, reproject +from rasterio.warp import ( + Resampling, + aligned_target, + calculate_default_transform, + reproject, +) def translate_profile( @@ -33,19 +38,21 @@ def translate_profile( dict Rasterio profile with transform shifted """ - transform = profile['transform'] + transform = profile["transform"] new_origin = transform * (x_shift, y_shift) new_transform = Affine.translation(*new_origin) new_transform = new_transform * transform.scale(transform.a, transform.e) p_new = profile.copy() - p_new['transform'] = new_transform + p_new["transform"] = new_transform return p_new -def translate_dataset(dataset: DatasetReader, x_shift: float, y_shift: float) -> tuple[MemoryFile, DatasetReader]: +def translate_dataset( + dataset: DatasetReader, x_shift: float, y_shift: float +) -> tuple[MemoryFile, DatasetReader]: """Create a new in-memory dataset and translates this. Closes the input dataset. Parameters @@ -78,7 +85,7 @@ def reproject_arr_to_match_profile( ref_profile: dict, nodata: Union[float, int] = None, num_threads: int = 1, - resampling: str = 'bilinear', + resampling: str = "bilinear", ) -> tuple[np.ndarray, dict]: """ Reproject an array to match a reference profile providing the reprojected array and the new profile. @@ -115,25 +122,25 @@ def reproject_arr_to_match_profile( (vertical dim.) x (horizontal dim), but output will be: 1 x (vertical dim.) x (horizontal dim). """ - dst_crs = ref_profile['crs'] - dst_transform = ref_profile['transform'] + dst_crs = ref_profile["crs"] + dst_transform = ref_profile["transform"] reproject_profile = ref_profile.copy() - nodata = nodata or src_profile['nodata'] - src_dtype = src_profile['dtype'] - count = src_profile['count'] + nodata = nodata or src_profile["nodata"] + src_dtype = src_profile["dtype"] + count = src_profile["count"] - reproject_profile.update({'dtype': src_dtype, 'nodata': nodata, 'count': count}) + reproject_profile.update({"dtype": src_dtype, "nodata": nodata, "count": count}) - height, width = ref_profile['height'], ref_profile['width'] + height, width = ref_profile["height"], ref_profile["width"] dst_array = np.zeros((count, height, width)) reproject( src_array, dst_array, - src_transform=src_profile['transform'], - src_crs=src_profile['crs'], + src_transform=src_profile["transform"], + src_crs=src_profile["crs"], dst_transform=dst_transform, dst_crs=dst_crs, dst_nodata=nodata, @@ -157,18 +164,20 @@ def get_bounds_dict(profile: dict) -> dict: dict: The bounds dictionary. """ - lx, ly = profile['width'], profile['height'] - transform = profile['transform'] + lx, ly = profile["width"], profile["height"] + transform = profile["transform"] bounds_dict = { - 'left': transform.c, - 'right': transform.c + transform.a * lx, - 'top': transform.f, - 'bottom': transform.f + transform.e * ly, + "left": transform.c, + "right": transform.c + transform.a * lx, + "top": transform.f, + "bottom": transform.f + transform.e * ly, } return bounds_dict -def reproject_profile_to_new_crs(src_profile: dict, dst_crs: CRS, target_resolution: Union[float, int] = None) -> dict: +def reproject_profile_to_new_crs( + src_profile: dict, dst_crs: CRS, target_resolution: Union[float, int] = None +) -> dict: """Create a new profile into a new CRS based on a dst_crs. May specify resolution. Parameters @@ -188,19 +197,21 @@ def reproject_profile_to_new_crs(src_profile: dict, dst_crs: CRS, target_resolut reprojected_profile = src_profile.copy() bounds_dict = get_bounds_dict(src_profile) - src_crs = src_profile['crs'] - w, h = src_profile['width'], src_profile['height'] - dst_trans, dst_w, dst_h = calculate_default_transform(src_crs, dst_crs, w, h, **bounds_dict) + src_crs = src_profile["crs"] + w, h = src_profile["width"], src_profile["height"] + dst_trans, dst_w, dst_h = calculate_default_transform( + src_crs, dst_crs, w, h, **bounds_dict + ) if target_resolution is not None: tr = target_resolution dst_trans, dst_w, dst_h = aligned_target(dst_trans, dst_w, dst_h, tr) reprojected_profile.update( { - 'crs': dst_crs, - 'transform': dst_trans, - 'width': dst_w, - 'height': dst_h, + "crs": dst_crs, + "transform": dst_trans, + "width": dst_w, + "height": dst_h, } ) return reprojected_profile @@ -210,7 +221,7 @@ def reproject_arr_to_new_crs( src_array: np.ndarray, src_profile: dict, dst_crs: str, - resampling: str = 'bilinear', + resampling: str = "bilinear", target_resolution: float = None, ) -> tuple[np.ndarray, dict]: """ @@ -236,20 +247,28 @@ def reproject_arr_to_new_crs( (reprojected_array, reprojected_profile) of data. """ tr = target_resolution - reprojected_profile = reproject_profile_to_new_crs(src_profile, dst_crs, target_resolution=tr) + reprojected_profile = reproject_profile_to_new_crs( + src_profile, dst_crs, target_resolution=tr + ) resampling = Resampling[resampling] - dst_array = np.zeros((reprojected_profile['count'], reprojected_profile['height'], reprojected_profile['width'])) + dst_array = np.zeros( + ( + reprojected_profile["count"], + reprojected_profile["height"], + reprojected_profile["width"], + ) + ) reproject( # Source parameters source=src_array, - src_crs=src_profile['crs'], - src_transform=src_profile['transform'], + src_crs=src_profile["crs"], + src_transform=src_profile["transform"], # Destination paramaters destination=dst_array, - dst_transform=reprojected_profile['transform'], - dst_crs=reprojected_profile['crs'], - dst_nodata=src_profile['nodata'], + dst_transform=reprojected_profile["transform"], + dst_crs=reprojected_profile["crs"], + dst_nodata=src_profile["nodata"], # Configuration resampling=resampling, ) @@ -297,16 +316,20 @@ def _aligned_target( return dst_transform, dst_width, dst_height -def update_profile_resolution(src_profile: dict, resolution: Union[float, tuple[float]]) -> dict: - transform = src_profile['transform'] - width = src_profile['width'] - height = src_profile['height'] +def update_profile_resolution( + src_profile: dict, resolution: Union[float, tuple[float]] +) -> dict: + transform = src_profile["transform"] + width = src_profile["width"] + height = src_profile["height"] - dst_transform, dst_width, dst_height = _aligned_target(transform, width, height, resolution) + dst_transform, dst_width, dst_height = _aligned_target( + transform, width, height, resolution + ) dst_profile = src_profile.copy() - dst_profile['width'] = dst_width - dst_profile['height'] = dst_height - dst_profile['transform'] = dst_transform + dst_profile["width"] = dst_width + dst_profile["height"] = dst_height + dst_profile["transform"] = dst_transform - return dst_profile \ No newline at end of file + return dst_profile diff --git a/sar_antarctica/utils/spatial.py b/sar_antarctica/utils/spatial.py index c7bfe81..8dcbea6 100644 --- a/sar_antarctica/utils/spatial.py +++ b/sar_antarctica/utils/spatial.py @@ -6,11 +6,14 @@ from pyproj import CRS import logging -def transform_polygon(geometry: Polygon, src_crs: int, dst_crs: int, always_xy: bool = True): + +def transform_polygon( + geometry: Polygon, src_crs: int, dst_crs: int, always_xy: bool = True +): src_crs = pyproj.CRS(f"EPSG:{src_crs}") - dst_crs = pyproj.CRS(f"EPSG:{dst_crs}") + dst_crs = pyproj.CRS(f"EPSG:{dst_crs}") transformer = pyproj.Transformer.from_crs(src_crs, dst_crs, always_xy=always_xy) - # Transform the polygon's coordinates + # Transform the polygon's coordinates if isinstance(geometry, Polygon): # Transform exterior exterior_coords = [ @@ -27,19 +30,22 @@ def transform_polygon(geometry: Polygon, src_crs: int, dst_crs: int, always_xy: # Handle other geometry types as needed raise ValueError("Only Polygon geometries are supported for transformation.") -def adjust_bounds(bounds: tuple, src_crs: int, ref_crs: int, segment_length: float = 0.1) -> tuple: + +def adjust_bounds( + bounds: tuple, src_crs: int, ref_crs: int, segment_length: float = 0.1 +) -> tuple: """ Adjust the bounding box around a scene in src_crs (4326) due to warping at high Latitudes. For example, the min and max boudning values for an antarctic scene in - 4326 may not actually be the true min and max due to distortions at high latitudes. + 4326 may not actually be the true min and max due to distortions at high latitudes. Parameters: - bounds: bounds to adjust. - src_crs: Source EPSG. e.g. 4326 - - ref_crs: reference crs to create the true bbox. i.e. 3031 in southern + - ref_crs: reference crs to create the true bbox. i.e. 3031 in southern hemisphere and 3995 in northern (polar stereographic) - segment_length: distance between generation points along the bounding box sides in - src_crs. e.g. 0.1 degrees in lat/lon + src_crs. e.g. 0.1 degrees in lat/lon Returns: - A polygon bounding box expanded to the true min max @@ -54,8 +60,8 @@ def adjust_bounds(bounds: tuple, src_crs: int, ref_crs: int, segment_length: flo def get_local_utm(bounds, antimeridian=False): - centre_lat = (bounds[1] + bounds[3])/2 - centre_lon = (bounds[0] + bounds[2])/2 + centre_lat = (bounds[1] + bounds[3]) / 2 + centre_lon = (bounds[0] + bounds[2]) / 2 if antimeridian: # force the lon to be next to antimeridian on the side with the scene centre. # e.g. (-177 + 178)/2 = 1, this is > 0 more data on -'ve side @@ -63,12 +69,12 @@ def get_local_utm(bounds, antimeridian=False): utm_crs_list = query_utm_crs_info( datum_name="WGS 84", area_of_interest=AreaOfInterest( - west_lon_degree=centre_lon-0.01, - south_lat_degree=centre_lat-0.01, - east_lon_degree=centre_lon+0.01, - north_lat_degree=centre_lat+0.01, + west_lon_degree=centre_lon - 0.01, + south_lat_degree=centre_lat - 0.01, + east_lon_degree=centre_lon + 0.01, + north_lat_degree=centre_lat + 0.01, ), ) crs = CRS.from_epsg(utm_crs_list[0].code) - crs = str(crs).split(':')[-1] # get the EPSG integer + crs = str(crs).split(":")[-1] # get the EPSG integer return int(crs) diff --git a/tests/filesystem/test_filesystem.py b/tests/filesystem/test_filesystem.py index 5220d69..80954cc 100644 --- a/tests/filesystem/test_filesystem.py +++ b/tests/filesystem/test_filesystem.py @@ -6,6 +6,7 @@ from pathlib import Path import pytest + @dataclasses.dataclass class Scene: id: str @@ -20,35 +21,58 @@ class Scene: scene_1 = Scene( id="S1A_EW_GRDM_1SDH_20220612T120348_20220612T120452_043629_053582_0F66", - file=Path("/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2022/2022-06/65S115E-70S120E/S1A_EW_GRDM_1SDH_20220612T120348_20220612T120452_043629_053582_0F66.zip"), + file=Path( + "/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2022/2022-06/65S115E-70S120E/S1A_EW_GRDM_1SDH_20220612T120348_20220612T120452_043629_053582_0F66.zip" + ), sensor="S1A", - start_date=datetime(2022,6,12,12,3,48), - stop_date=datetime(2022,6,12,12,4,52), - latest_orbit=Path("/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1A/S1A_OPER_AUX_POEORB_OPOD_20220702T081845_V20220611T225942_20220613T005942.EOF"), - latest_poe_orbit=Path("/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1A/S1A_OPER_AUX_POEORB_OPOD_20220702T081845_V20220611T225942_20220613T005942.EOF"), - latest_res_orbit=Path("/g/data/fj7/Copernicus/Sentinel-1/RESORB/S1A/S1A_OPER_AUX_RESORB_OPOD_20220612T143829_V20220612T104432_20220612T140202.EOF"), + start_date=datetime(2022, 6, 12, 12, 3, 48), + stop_date=datetime(2022, 6, 12, 12, 4, 52), + latest_orbit=Path( + "/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1A/S1A_OPER_AUX_POEORB_OPOD_20220702T081845_V20220611T225942_20220613T005942.EOF" + ), + latest_poe_orbit=Path( + "/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1A/S1A_OPER_AUX_POEORB_OPOD_20220702T081845_V20220611T225942_20220613T005942.EOF" + ), + latest_res_orbit=Path( + "/g/data/fj7/Copernicus/Sentinel-1/RESORB/S1A/S1A_OPER_AUX_RESORB_OPOD_20220612T143829_V20220612T104432_20220612T140202.EOF" + ), ) scene_2 = Scene( id="S1B_EW_GRDM_1SDH_20191130T165626_20191130T165726_019159_0242A2_2F58", - file=Path("/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2019/2019-11/65S160E-70S165E/S1B_EW_GRDM_1SDH_20191130T165626_20191130T165726_019159_0242A2_2F58.zip"), + file=Path( + "/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2019/2019-11/65S160E-70S165E/S1B_EW_GRDM_1SDH_20191130T165626_20191130T165726_019159_0242A2_2F58.zip" + ), sensor="S1B", - start_date=datetime(2019,11,30,16,56,26), - stop_date=datetime(2019,11,30,16,57,26), - latest_orbit=Path("/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1B/S1B_OPER_AUX_POEORB_OPOD_20191220T110516_V20191129T225942_20191201T005942.EOF"), - latest_poe_orbit=Path("/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1B/S1B_OPER_AUX_POEORB_OPOD_20191220T110516_V20191129T225942_20191201T005942.EOF"), - latest_res_orbit=Path("/g/data/fj7/Copernicus/Sentinel-1/RESORB/S1B/S1B_OPER_AUX_RESORB_OPOD_20191130T210136_V20191130T154804_20191130T190534.EOF"), + start_date=datetime(2019, 11, 30, 16, 56, 26), + stop_date=datetime(2019, 11, 30, 16, 57, 26), + latest_orbit=Path( + "/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1B/S1B_OPER_AUX_POEORB_OPOD_20191220T110516_V20191129T225942_20191201T005942.EOF" + ), + latest_poe_orbit=Path( + "/g/data/fj7/Copernicus/Sentinel-1/POEORB/S1B/S1B_OPER_AUX_POEORB_OPOD_20191220T110516_V20191129T225942_20191201T005942.EOF" + ), + latest_res_orbit=Path( + "/g/data/fj7/Copernicus/Sentinel-1/RESORB/S1B/S1B_OPER_AUX_RESORB_OPOD_20191130T210136_V20191130T154804_20191130T190534.EOF" + ), ) scenes = [scene_1, scene_2] + @pytest.mark.parametrize("scene", scenes) def test_find_latest_orbit_for_scene(scene: Scene): assert find_latest_orbit_for_scene(scene.id) == scene.latest_orbit - assert find_latest_orbit_for_scene(scene.id, orbit_type="RES") == scene.latest_res_orbit - assert find_latest_orbit_for_scene(scene.id, orbit_type="POE") == scene.latest_poe_orbit + assert ( + find_latest_orbit_for_scene(scene.id, orbit_type="RES") + == scene.latest_res_orbit + ) + assert ( + find_latest_orbit_for_scene(scene.id, orbit_type="POE") + == scene.latest_poe_orbit + ) @pytest.mark.parametrize("scene", scenes) def test_find_scene_file_from_id(scene: Scene): - assert find_scene_file_from_id(scene.id) == scene.file \ No newline at end of file + assert find_scene_file_from_id(scene.id) == scene.file diff --git a/tests/sar_antarctica/test_dem.py b/tests/sar_antarctica/test_dem.py index 761ae0d..a79e029 100644 --- a/tests/sar_antarctica/test_dem.py +++ b/tests/sar_antarctica/test_dem.py @@ -14,241 +14,351 @@ expand_bounds, get_cop30_dem_for_bounds, find_required_dem_paths_from_index, - find_required_dem_tile_paths_by_filename -) -from sar_antarctica.nci.preparation.create_dem_vrt import ( - find_tiles, - build_tileindex + find_required_dem_tile_paths_by_filename, ) +from sar_antarctica.nci.preparation.create_dem_vrt import find_tiles, build_tileindex from sar_antarctica.nci.preparation.geoid import remove_geoid from sar_antarctica.utils.raster import ( merge_raster_files, bounds_from_profile, - expand_raster_to_bounds + expand_raster_to_bounds, ) from data.cop30m_profile import TEST_COP30_PROFILE_1, TEST_COP30_PROFILE_2 CURRENT_DIR = Path(os.path.abspath(__file__)).parent.resolve() -TEST_COP30_FOLDER_PATH = CURRENT_DIR / Path('data/copernicus_30m_world') -TEST_COP30_INDEX_PATH = CURRENT_DIR / Path('data/copernicus_30m_world/copdem_tindex_test.gpkg') -TEST_GEOID_TIF_PATH = CURRENT_DIR / Path('data/geoid/us_nga_egm2008_1_4326__agisoft_clipped.tif') +TEST_COP30_FOLDER_PATH = CURRENT_DIR / Path("data/copernicus_30m_world") +TEST_COP30_INDEX_PATH = CURRENT_DIR / Path( + "data/copernicus_30m_world/copdem_tindex_test.gpkg" +) +TEST_GEOID_TIF_PATH = CURRENT_DIR / Path( + "data/geoid/us_nga_egm2008_1_4326__agisoft_clipped.tif" +) # TEST_COP30_FOLDER_PATH = Path('/g/data/v10/eoancillarydata-2/elevation/copernicus_30m_world/') -# TEST_COP30_INDEX_PATH = Path('/g/data/yp75/projects/ancillary/dem/copdem_tindex.gpkg') # disable test_make_test_gpkg +# TEST_COP30_INDEX_PATH = Path('/g/data/yp75/projects/ancillary/dem/copdem_tindex.gpkg') # disable test_make_test_gpkg # TEST_GEOID_TIF_PATH = Path('/g/data/yp75/projects/ancillary/geoid/us_nga_egm2008_1_4326__agisoft.tif') + def test_pytest(): assert True + def test_make_test_gpkg(): """required to pass for other tests utilising index file""" - tiles = find_tiles(TEST_COP30_FOLDER_PATH, pattern='*.tif') + tiles = find_tiles(TEST_COP30_FOLDER_PATH, pattern="*.tif") build_tileindex(tiles, TEST_COP30_INDEX_PATH) assert os.path.exists(TEST_COP30_INDEX_PATH) -@pytest.mark.parametrize("bounds, expected", [ - ((163.121597, -78.632782, 172.382263, -76.383263), False), # Bounds that do not cross the antimeridian - ((170, -70, 180, -76), False), # Bounds that do not cross the antimeridian - ((-180, 10, -175, 20), False), # Bounds that do not cross the antimeridian - ((-177.884048, -78.176201, 178.838364, -75.697151), True), # Bounds that cross the antimeridian -]) + +@pytest.mark.parametrize( + "bounds, expected", + [ + ( + (163.121597, -78.632782, 172.382263, -76.383263), + False, + ), # Bounds that do not cross the antimeridian + ((170, -70, 180, -76), False), # Bounds that do not cross the antimeridian + ((-180, 10, -175, 20), False), # Bounds that do not cross the antimeridian + ( + (-177.884048, -78.176201, 178.838364, -75.697151), + True, + ), # Bounds that cross the antimeridian + ], +) def test_check_s1_bounds_cross_antimeridian(bounds, expected): assert check_s1_bounds_cross_antimeridian(bounds) == expected -@pytest.mark.parametrize("bounds, lat_buff, expected_left, expected_right", [ - ((-177.884048, -78.176201, 178.838364, -75.697151), 0, (-180, -78.176201, -177.884048, -75.697151), (178.838364, -78.176201, 180, -75.697151)), - ((-177.884048, -78.176201, 178.838364, -75.697151), 0.1, (-180, -78.276201, -177.884048, -75.597151), (178.838364, -78.276201, 180, -75.597151)), -]) + +@pytest.mark.parametrize( + "bounds, lat_buff, expected_left, expected_right", + [ + ( + (-177.884048, -78.176201, 178.838364, -75.697151), + 0, + (-180, -78.176201, -177.884048, -75.697151), + (178.838364, -78.176201, 180, -75.697151), + ), + ( + (-177.884048, -78.176201, 178.838364, -75.697151), + 0.1, + (-180, -78.276201, -177.884048, -75.597151), + (178.838364, -78.276201, 180, -75.597151), + ), + ], +) def test_split_bounds_at_am_crossing(bounds, lat_buff, expected_left, expected_right): left, right = split_s1_bounds_at_am_crossing(bounds, lat_buff) # use pytest.approx to get around small rounding errors with floats - assert all(a == pytest.approx(b,rel=1e-9) for a, b in zip(left, expected_left)) - assert all(a == pytest.approx(b,rel=1e-9) for a, b in zip(right, expected_right)) - -@pytest.mark.parametrize("bounds, target_crs", [ - ((-177.884048, -78.176201, 178.838364, -75.697151), 3031), # antarctic - ((-177.884048, 78.176201, 178.838364, 75.697151), 3995), # arctic - ((-177.884048, -1, 178.838364, 1), 32601), # centre slight left of equator - ((-178.884048, -1, 177.838364, 1), 32660), # centre slight right of equator -]) + assert all(a == pytest.approx(b, rel=1e-9) for a, b in zip(left, expected_left)) + assert all(a == pytest.approx(b, rel=1e-9) for a, b in zip(right, expected_right)) + + +@pytest.mark.parametrize( + "bounds, target_crs", + [ + ((-177.884048, -78.176201, 178.838364, -75.697151), 3031), # antarctic + ((-177.884048, 78.176201, 178.838364, 75.697151), 3995), # arctic + ((-177.884048, -1, 178.838364, 1), 32601), # centre slight left of equator + ((-178.884048, -1, 177.838364, 1), 32660), # centre slight right of equator + ], +) def test_get_target_antimeridian_projection(bounds, target_crs): assert get_target_antimeridian_projection(bounds) == target_crs -@pytest.mark.parametrize("bounds_to_profile", [ - TEST_COP30_PROFILE_1, - TEST_COP30_PROFILE_2 - ]) + +@pytest.mark.parametrize( + "bounds_to_profile", [TEST_COP30_PROFILE_1, TEST_COP30_PROFILE_2] +) def test_make_empty_cop30m_profile(bounds_to_profile): bounds, target_profile = bounds_to_profile profile = make_empty_cop30m_profile(bounds) for key in profile.keys(): if ( - isinstance(profile[key], float) and - math.isnan(profile[key]) and - isinstance(target_profile[key], float) + isinstance(profile[key], float) + and math.isnan(profile[key]) + and isinstance(target_profile[key], float) and math.isnan(target_profile[key]) - ): + ): # Handle NaN comparison continue assert profile[key] == target_profile[key] -@pytest.mark.parametrize("bounds, trg_shape", [ - ((-179.9, -79.2, -179.1, -79.1), (1,362,962)) -]) + +@pytest.mark.parametrize( + "bounds, trg_shape", [((-179.9, -79.2, -179.1, -79.1), (1, 362, 962))] +) def test_make_empty_cop30_dem(bounds, trg_shape): fill_value = 0 empty_dem_profile = make_empty_cop30m_profile(bounds) dem_arr, dem_profile = expand_raster_to_bounds( - bounds, - src_profile=empty_dem_profile, + bounds, + src_profile=empty_dem_profile, fill_value=fill_value, buffer_pixels=1, - save_path='empty.tif' - ) + save_path="empty.tif", + ) dem_bounds = bounds_from_profile(dem_profile) assert box(*bounds).within(box(*dem_bounds)) - assert dem_arr.shape == trg_shape - -@pytest.mark.parametrize("bounds, buffer, expanded_bounds", [ - ((179.1, -79.2, 179.9, -79.1), 0, (179.09161, -79.201308, 179.900922, -79.09867)), - ((179.1, -79.2, 179.9, -79.1), 0.01, (179.08161, -79.21130, 179.91092, -79.08867)), - ((10, -65.2, 20, -60), 0, (8.23494, -66.27035, 23.91474, -58.78356)), - ((10, -65.2, 20, -60), 0, (8.23494, -66.27035, 23.91474, -58.78356)), - ((-90, -85, -80, -80), 0, (-90, -85.07587, -70.54166, -79.85109)) -]) + assert dem_arr.shape == trg_shape + + +@pytest.mark.parametrize( + "bounds, buffer, expanded_bounds", + [ + ( + (179.1, -79.2, 179.9, -79.1), + 0, + (179.09161, -79.201308, 179.900922, -79.09867), + ), + ( + (179.1, -79.2, 179.9, -79.1), + 0.01, + (179.08161, -79.21130, 179.91092, -79.08867), + ), + ((10, -65.2, 20, -60), 0, (8.23494, -66.27035, 23.91474, -58.78356)), + ((10, -65.2, 20, -60), 0, (8.23494, -66.27035, 23.91474, -58.78356)), + ((-90, -85, -80, -80), 0, (-90, -85.07587, -70.54166, -79.85109)), + ], +) def test_expand_bounds(bounds, buffer, expanded_bounds): new_bounds = expand_bounds(bounds, buffer=buffer) - assert pytest.approx(new_bounds[0],rel=1e-5) == pytest.approx(expanded_bounds[0],rel=1e-5) - assert pytest.approx(new_bounds[1],rel=1e-5) == pytest.approx(expanded_bounds[1],rel=1e-5) - assert pytest.approx(new_bounds[2],rel=1e-5) == pytest.approx(expanded_bounds[2],rel=1e-5) - assert pytest.approx(new_bounds[3],rel=1e-5) == pytest.approx(expanded_bounds[3],rel=1e-5) - -@pytest.mark.parametrize("bounds, search_buffer, expected_tiles", [ - ((176.1, -79.2, 176.9, -79.1), 0, []), - ((179.1, -79.2, 179.9, -79.1), 0, [str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif')]), - ((179.1, -79.2, 179.9, -79.1), 0.3, [ - str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E178_00_DEM.tif'), - str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif')]), - ((178.1, -79.2, 179.9, -79.1), 0, [ - str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E178_00_DEM.tif'), - str(TEST_COP30_FOLDER_PATH / 'Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif')]), -]) + assert pytest.approx(new_bounds[0], rel=1e-5) == pytest.approx( + expanded_bounds[0], rel=1e-5 + ) + assert pytest.approx(new_bounds[1], rel=1e-5) == pytest.approx( + expanded_bounds[1], rel=1e-5 + ) + assert pytest.approx(new_bounds[2], rel=1e-5) == pytest.approx( + expanded_bounds[2], rel=1e-5 + ) + assert pytest.approx(new_bounds[3], rel=1e-5) == pytest.approx( + expanded_bounds[3], rel=1e-5 + ) + + +@pytest.mark.parametrize( + "bounds, search_buffer, expected_tiles", + [ + ((176.1, -79.2, 176.9, -79.1), 0, []), + ( + (179.1, -79.2, 179.9, -79.1), + 0, + [ + str( + TEST_COP30_FOLDER_PATH + / "Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif" + ) + ], + ), + ( + (179.1, -79.2, 179.9, -79.1), + 0.3, + [ + str( + TEST_COP30_FOLDER_PATH + / "Copernicus_DSM_COG_10_S80_00_E178_00_DEM.tif" + ), + str( + TEST_COP30_FOLDER_PATH + / "Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif" + ), + ], + ), + ( + (178.1, -79.2, 179.9, -79.1), + 0, + [ + str( + TEST_COP30_FOLDER_PATH + / "Copernicus_DSM_COG_10_S80_00_E178_00_DEM.tif" + ), + str( + TEST_COP30_FOLDER_PATH + / "Copernicus_DSM_COG_10_S80_00_E179_00_DEM.tif" + ), + ], + ), + ], +) def test_find_required_dem_paths(bounds, search_buffer, expected_tiles): # limited to files in TEST_COP30_FOLDER_PATH tiles_from_folder = find_required_dem_tile_paths_by_filename( bounds, cop30_folder_path=TEST_COP30_FOLDER_PATH, search_buffer=search_buffer, - tifs_in_subfolder=False) + tifs_in_subfolder=False, + ) tiles_from_index = find_required_dem_paths_from_index( - bounds, - cop30_index_path=TEST_COP30_INDEX_PATH, - search_buffer=search_buffer) + bounds, cop30_index_path=TEST_COP30_INDEX_PATH, search_buffer=search_buffer + ) assert tiles_from_folder == expected_tiles assert tiles_from_index == expected_tiles assert tiles_from_folder == tiles_from_index -@pytest.mark.parametrize("bounds, trg_shape, buffer_pixels", [ - ((-179.9, -79.2, -179.1, -79.1), (1,362,962), 0), - ((-179.9, -79.2, -179.1, -79.1), (1,370,970), 4), - ((-179.6, -79.9, -179.4, -79.5), (1, 1442, 242), 0), - ((179.1, -79.2, 179.9, -79.1), (1,362,962), 0), - ((179.5, -79.2, 179.6, -79.01), (1,690,126), 2), - ((179.5, -79.2, 179.6, -79.01), (1,688,124), 1), - ((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 0), - ((178.1, -79.2, 179.95, -79.1), (1, 366, 2226), 2), - ]) +@pytest.mark.parametrize( + "bounds, trg_shape, buffer_pixels", + [ + ((-179.9, -79.2, -179.1, -79.1), (1, 362, 962), 0), + ((-179.9, -79.2, -179.1, -79.1), (1, 370, 970), 4), + ((-179.6, -79.9, -179.4, -79.5), (1, 1442, 242), 0), + ((179.1, -79.2, 179.9, -79.1), (1, 362, 962), 0), + ((179.5, -79.2, 179.6, -79.01), (1, 690, 126), 2), + ((179.5, -79.2, 179.6, -79.01), (1, 688, 124), 1), + ((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 0), + ((178.1, -79.2, 179.95, -79.1), (1, 366, 2226), 2), + ], +) def test_dem_read_for_bounds(bounds, trg_shape, buffer_pixels): - os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True) - dem_tiles = find_required_dem_paths_from_index(bounds, cop30_index_path=TEST_COP30_INDEX_PATH) + os.makedirs(CURRENT_DIR / Path("TMP"), exist_ok=True) + dem_tiles = find_required_dem_paths_from_index( + bounds, cop30_index_path=TEST_COP30_INDEX_PATH + ) dem_arr, dem_profile = merge_raster_files( dem_tiles, - output_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'), + output_path=CURRENT_DIR / Path("TMP") / Path("TMP.tif"), bounds=bounds, - buffer_pixels=buffer_pixels - ) - shutil.rmtree(CURRENT_DIR / Path('TMP')) + buffer_pixels=buffer_pixels, + ) + shutil.rmtree(CURRENT_DIR / Path("TMP")) dem_bounds = bounds_from_profile(dem_profile) assert box(*bounds).within(box(*dem_bounds)) assert dem_arr.shape == trg_shape -@pytest.mark.parametrize("bounds, trg_shape, geoid_ref_mean, ellipsoid_ref_mean", [ - ((-179.9, -79.2, -179.1, -79.1), (1,362,962), 44.088665, -9.830775), - ((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 38.270348, -15.338912), - ]) + +@pytest.mark.parametrize( + "bounds, trg_shape, geoid_ref_mean, ellipsoid_ref_mean", + [ + ((-179.9, -79.2, -179.1, -79.1), (1, 362, 962), 44.088665, -9.830775), + ((178.1, -79.2, 179.95, -79.1), (1, 362, 2222), 38.270348, -15.338912), + ], +) def test_remove_geoid(bounds, trg_shape, geoid_ref_mean, ellipsoid_ref_mean): - os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True) - dem_tiles = find_required_dem_paths_from_index(bounds, cop30_index_path=TEST_COP30_INDEX_PATH) + os.makedirs(CURRENT_DIR / Path("TMP"), exist_ok=True) + dem_tiles = find_required_dem_paths_from_index( + bounds, cop30_index_path=TEST_COP30_INDEX_PATH + ) dem_arr, dem_profile = merge_raster_files( dem_tiles, - output_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'), + output_path=CURRENT_DIR / Path("TMP") / Path("TMP.tif"), bounds=bounds, - buffer_pixels=0 - ) + buffer_pixels=0, + ) dem_arr_ellipsoid = remove_geoid( - dem_arr = dem_arr, - dem_profile = dem_profile, - geoid_path = TEST_GEOID_TIF_PATH, - dem_area_or_point = 'Point', - buffer_pixels = 2, - save_path='', + dem_arr=dem_arr, + dem_profile=dem_profile, + geoid_path=TEST_GEOID_TIF_PATH, + dem_area_or_point="Point", + buffer_pixels=2, + save_path="", ) - shutil.rmtree(CURRENT_DIR / Path('TMP')) + shutil.rmtree(CURRENT_DIR / Path("TMP")) assert dem_arr.shape == dem_arr_ellipsoid.shape assert dem_arr.shape == trg_shape - assert np.mean(dem_arr) == pytest.approx(geoid_ref_mean,rel=1e-5) - assert np.mean(dem_arr_ellipsoid) == pytest.approx(ellipsoid_ref_mean,rel=1e-5) - -@pytest.mark.parametrize("bounds, ellipsoid_heights, trg_shape, trg_crs", [ - ((-179.9, -79.2, -179.1, -79.1), False, (1, 366, 966), 4326), - ((178.1, -79.2, 179.95, -79.1), False, (1, 366, 2226), 4326), - ((-179.2, -79.2, 179.1, -79.1), False, (1, 1396, 1676), 3031), # across antimeridian - ((-179.2, -79.2, 179.1, -79.1), True, (1, 1396, 1676), 3031) # across antimeridian - ]) + assert np.mean(dem_arr) == pytest.approx(geoid_ref_mean, rel=1e-5) + assert np.mean(dem_arr_ellipsoid) == pytest.approx(ellipsoid_ref_mean, rel=1e-5) + + +@pytest.mark.parametrize( + "bounds, ellipsoid_heights, trg_shape, trg_crs", + [ + ((-179.9, -79.2, -179.1, -79.1), False, (1, 366, 966), 4326), + ((178.1, -79.2, 179.95, -79.1), False, (1, 366, 2226), 4326), + ( + (-179.2, -79.2, 179.1, -79.1), + False, + (1, 1396, 1676), + 3031, + ), # across antimeridian + ( + (-179.2, -79.2, 179.1, -79.1), + True, + (1, 1396, 1676), + 3031, + ), # across antimeridian + ], +) def test_get_cop30_dem_for_bounds(bounds, ellipsoid_heights, trg_shape, trg_crs): - os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True) + os.makedirs(CURRENT_DIR / Path("TMP"), exist_ok=True) dem_arr, dem_profile = get_cop30_dem_for_bounds( - bounds, - save_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'), - ellipsoid_heights=ellipsoid_heights, + bounds, + save_path=CURRENT_DIR / Path("TMP") / Path("TMP.tif"), + ellipsoid_heights=ellipsoid_heights, buffer_pixels=2, cop30_index_path=TEST_COP30_INDEX_PATH, cop30_folder_path=TEST_COP30_FOLDER_PATH, geoid_tif_path=TEST_GEOID_TIF_PATH, - adjust_for_high_lat_and_buffer=False + adjust_for_high_lat_and_buffer=False, ) - #shutil.rmtree(CURRENT_DIR / Path('TMP')) + # shutil.rmtree(CURRENT_DIR / Path('TMP')) assert dem_arr.shape == trg_shape - assert dem_profile['crs'].to_epsg() == trg_crs + assert dem_profile["crs"].to_epsg() == trg_crs + if __name__ == "__main__": - os.makedirs(CURRENT_DIR / Path('TMP'), exist_ok=True) - bounds = (-179.9, -79.2, -179.1, -79.1) + os.makedirs(CURRENT_DIR / Path("TMP"), exist_ok=True) + bounds = (-179.9, -79.2, -179.1, -79.1) dem_arr, dem_profile = get_cop30_dem_for_bounds( - bounds, - save_path= CURRENT_DIR / Path('TMP') / Path('TMP.tif'), - ellipsoid_heights=True, + bounds, + save_path=CURRENT_DIR / Path("TMP") / Path("TMP.tif"), + ellipsoid_heights=True, buffer_pixels=0, cop30_index_path=TEST_COP30_INDEX_PATH, cop30_folder_path=TEST_COP30_FOLDER_PATH, geoid_tif_path=TEST_GEOID_TIF_PATH, - adjust_for_high_lat_and_buffer=False + adjust_for_high_lat_and_buffer=False, ) print(dem_profile) empty_dem_profile = make_empty_cop30m_profile((0, -90, 1, -86)) print(empty_dem_profile) dem_arr, dem_profile = expand_raster_to_bounds( - bounds, - src_profile=empty_dem_profile, + bounds, + src_profile=empty_dem_profile, fill_value=0, buffer_pixels=1, - save_path=CURRENT_DIR / Path('TMP') / Path('EMPTY.tif'), - ) + save_path=CURRENT_DIR / Path("TMP") / Path("EMPTY.tif"), + ) dem_bounds = bounds_from_profile(dem_profile) print(dem_profile) - - - diff --git a/tests/sar_antarctica/test_dem_vrt.py b/tests/sar_antarctica/test_dem_vrt.py index 8ba24a5..c232961 100644 --- a/tests/sar_antarctica/test_dem_vrt.py +++ b/tests/sar_antarctica/test_dem_vrt.py @@ -1,10 +1,14 @@ -from sar_antarctica.nci.preparation.create_dem_vrt import find_tiles, build_vrt, build_tileindex +from sar_antarctica.nci.preparation.create_dem_vrt import ( + find_tiles, + build_vrt, + build_tileindex, +) from pathlib import Path import os CURRENT_DIR = Path(__file__).parent.resolve() -TEST_DATA_PATH = CURRENT_DIR / Path('data/copernicus_30m_world') +TEST_DATA_PATH = CURRENT_DIR / Path("data/copernicus_30m_world") EXPECTED_DATA_FILES = [ TEST_DATA_PATH / "Copernicus_DSM_COG_10_S80_00_E178_00_DEM.tif", @@ -14,22 +18,23 @@ EXPECTED_VRT_FILE = TEST_DATA_PATH / "copdem_test.vrt" TEST_FILE_LIST_PATH = Path("temp.txt") -TEST_VRT_PATH = TEST_DATA_PATH / 'temp.vrt' -TEST_TINDEX_PATH = TEST_DATA_PATH / 'temp.gpkg' +TEST_VRT_PATH = TEST_DATA_PATH / "temp.vrt" +TEST_TINDEX_PATH = TEST_DATA_PATH / "temp.gpkg" tiles = find_tiles(TEST_DATA_PATH, "Copernicus_DSM_COG_10_???_00_????_00_DEM.tif") list_tiles = list(tiles) + def test_find_tiles_for_vrt(): assert set(list_tiles) == set(EXPECTED_DATA_FILES) - + def test_build_vrt(): build_vrt(list_tiles, TEST_VRT_PATH, run=False) assert TEST_FILE_LIST_PATH.exists - with open(TEST_FILE_LIST_PATH, 'r') as f: + with open(TEST_FILE_LIST_PATH, "r") as f: lines = [Path(line.rstrip()) for line in f.readlines()] assert lines == list_tiles @@ -43,7 +48,7 @@ def test_build_tindex(): build_tileindex(list_tiles, TEST_TINDEX_PATH, run=False) assert TEST_FILE_LIST_PATH.exists - with open(TEST_FILE_LIST_PATH, 'r') as f: + with open(TEST_FILE_LIST_PATH, "r") as f: lines = [Path(line.rstrip()) for line in f.readlines()] assert lines == list_tiles @@ -51,4 +56,3 @@ def test_build_tindex(): assert TEST_TINDEX_PATH.exists os.remove(TEST_TINDEX_PATH) - diff --git a/tests/sar_antarctica/test_orbits.py b/tests/sar_antarctica/test_orbits.py index 6c62f40..f164bfc 100644 --- a/tests/sar_antarctica/test_orbits.py +++ b/tests/sar_antarctica/test_orbits.py @@ -5,6 +5,7 @@ import dataclasses from datetime import datetime + @dataclasses.dataclass class Orbit: file: str @@ -12,21 +13,23 @@ class Orbit: start_date: datetime stop_date: datetime + orbit_1 = Orbit( file="S1A_OPER_AUX_POEORB_OPOD_20141207T123431_V20141115T225944_20141117T005944.EOF", - published_date=datetime(2014, 12, 7,12,34,31), - start_date=datetime(2014,11,15,22,59,44), - stop_date=datetime(2014,11,17,0,59,44) + published_date=datetime(2014, 12, 7, 12, 34, 31), + start_date=datetime(2014, 11, 15, 22, 59, 44), + stop_date=datetime(2014, 11, 17, 0, 59, 44), ) orbit_2 = Orbit( file="S1A_OPER_AUX_POEORB_OPOD_20191220T120706_V20191129T225942_20191201T005942.EOF", - published_date=datetime(2019,12,20,12,7,6), - start_date=datetime(2019,11,29,22,59,42), - stop_date=datetime(2019,12,1,0,59,42) + published_date=datetime(2019, 12, 20, 12, 7, 6), + start_date=datetime(2019, 11, 29, 22, 59, 42), + stop_date=datetime(2019, 12, 1, 0, 59, 42), ) orbits = [orbit_1, orbit_2] + @pytest.mark.parametrize("orbit", orbits) def test_parse_orbit_file_dates(orbit: Orbit): date_tuple = (orbit.published_date, orbit.start_date, orbit.stop_date) diff --git a/tests/sar_antarctica/test_scenes.py b/tests/sar_antarctica/test_scenes.py index 1f85849..28ca97e 100644 --- a/tests/sar_antarctica/test_scenes.py +++ b/tests/sar_antarctica/test_scenes.py @@ -9,6 +9,7 @@ from pathlib import Path import pytest + @dataclasses.dataclass class Scene: id: str @@ -17,24 +18,30 @@ class Scene: start_date: datetime stop_date: datetime + scene_1 = Scene( id="S1A_EW_GRDM_1SDH_20200330T165825_20200330T165929_031907_03AF02_8570", - file=Path("/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2020/2020-03/70S050E-75S055E/S1A_EW_GRDM_1SDH_20200330T165825_20200330T165929_031907_03AF02_8570.zip"), + file=Path( + "/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2020/2020-03/70S050E-75S055E/S1A_EW_GRDM_1SDH_20200330T165825_20200330T165929_031907_03AF02_8570.zip" + ), sensor="S1A", - start_date=datetime(2020,3,30,16,58,25), - stop_date=datetime(2020,3,30,16,59,29) + start_date=datetime(2020, 3, 30, 16, 58, 25), + stop_date=datetime(2020, 3, 30, 16, 59, 29), ) scene_2 = Scene( id="S1B_EW_GRDM_1SDH_20210914T112333_20210914T112403_028693_036C96_3EA8", - file=Path("/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2021/2021-09/60S120E-65S125E/S1B_EW_GRDM_1SDH_20210914T112333_20210914T112403_028693_036C96_3EA8.zip"), + file=Path( + "/g/data/fj7/Copernicus/Sentinel-1/C-SAR/GRD/2021/2021-09/60S120E-65S125E/S1B_EW_GRDM_1SDH_20210914T112333_20210914T112403_028693_036C96_3EA8.zip" + ), sensor="S1B", - start_date=datetime(2021,9,14,11,23,33), - stop_date=datetime(2021,9,14,11,24,3) + start_date=datetime(2021, 9, 14, 11, 23, 33), + stop_date=datetime(2021, 9, 14, 11, 24, 3), ) scenes = [scene_1, scene_2] + @pytest.mark.parametrize("scene", scenes) def test_parse_scene_file_dates(scene: Scene): date_tuple = (scene.start_date, scene.stop_date) @@ -44,5 +51,3 @@ def test_parse_scene_file_dates(scene: Scene): @pytest.mark.parametrize("scene", scenes) def test_parse_scene_file_sensor(scene: Scene): assert parse_scene_file_sensor(scene.id) == scene.sensor - - diff --git a/utils.py b/utils.py index 8579432..81cdbb6 100644 --- a/utils.py +++ b/utils.py @@ -2,11 +2,14 @@ from shapely import segmentize from shapely.geometry import Polygon, box -def transform_polygon(geometry: Polygon, src_crs: int, dst_crs: int, always_xy: bool = True): + +def transform_polygon( + geometry: Polygon, src_crs: int, dst_crs: int, always_xy: bool = True +): src_crs = pyproj.CRS(f"EPSG:{src_crs}") - dst_crs = pyproj.CRS(f"EPSG:{dst_crs}") + dst_crs = pyproj.CRS(f"EPSG:{dst_crs}") transformer = pyproj.Transformer.from_crs(src_crs, dst_crs, always_xy=always_xy) - # Transform the polygon's coordinates + # Transform the polygon's coordinates transformed_exterior = [ transformer.transform(x, y) for x, y in geometry.exterior.coords ] @@ -14,19 +17,22 @@ def transform_polygon(geometry: Polygon, src_crs: int, dst_crs: int, always_xy: transformed_polygon = Polygon(transformed_exterior) return transformed_polygon -def transform_scene_extent(geometry: Polygon, src_crs: int, ref_crs: int, segment_length: float = 0.1) -> Polygon: + +def transform_scene_extent( + geometry: Polygon, src_crs: int, ref_crs: int, segment_length: float = 0.1 +) -> Polygon: """ Adjust the bounding box around a scene in src_crs (4326) due to warping at high Latitudes. For example, the min and max boudning values for an antarctic scene in - 4326 may not actually be the true min and max due to distortions at high latitudes. + 4326 may not actually be the true min and max due to distortions at high latitudes. Parameters: - geometry: Polygon of the scene geometry. - src_crs: Source EPSG. e.g. 4326 - - ref_crs: reference crs to create the true bbox. i.e. 3031 in southern + - ref_crs: reference crs to create the true bbox. i.e. 3031 in southern hemisphere and 3995 in northern (polar stereographic) - segment_length: distance between generation points along the bounding box sides in - src_crs. e.g. 0.1 degrees in lat/lon + src_crs. e.g. 0.1 degrees in lat/lon Returns: - A polygon bounding box expanded to the true min max @@ -39,4 +45,4 @@ def transform_scene_extent(geometry: Polygon, src_crs: int, ref_crs: int, segmen corrected_geometry = transform_polygon(transformed_box, ref_crs, src_crs) - return corrected_geometry \ No newline at end of file + return corrected_geometry