diff --git a/.gitignore b/.gitignore index bbca6646..3a894b2c 100644 --- a/.gitignore +++ b/.gitignore @@ -240,3 +240,7 @@ Sessionx.vim tags # Persistent undo [._]*.un~ + +# Test data +tests/data/merge/* +tests/data/water_mask_output.* diff --git a/CHANGELOG.md b/CHANGELOG.md index 624466fb..ec4453d8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,20 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/) and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [1.0.0] +### Added +* `merge_tops_bursts.py` file and workflow for merge burst products created using insar_tops_bursts. +* `merge_tops_bursts` entrypoint +* `merge_tops_bursts` README template and creation functionality +* several classes and functions to `burst.py` and `utils.py` to support `merge_tops_burst`. +* tests for the added functionality. +* `tests/data/merge.zip` example data for testing merge workflow. +* `tests/data/create_merge_test_data.py` for generating merge workflow test data. + +### Changed +* `insar_tops_burst.py` to add four radar coordinate datasets to burst InSAR products (wrapped phase, lat, lon, los). +* README files generated in `insar_tops_burst.py` are now use blocks and extends the `insar_burst_base.md.txt.j2` jinja template. + ## [0.10.0] ### Added * Support for a new water masking dataset based off of OpenStreetMaps and ESA WorldCover data. diff --git a/README.md b/README.md index b4aa256e..c64afdab 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,34 @@ python -m hyp3_isce2 ++process insar_tops_burst \ This command will create a Sentinel-1 interferogram that contains a deformation signal related to a 2020 Iranian earthquake. +### Product Merging Utility Usage +**This feature is under active development and is subject to change!** + +Burst InSAR products created using the `insar_tops_burst` workflow can be merged together using the `merge_tops_burst` workflow. This can be useful when the deformation signal you'd like to observe spans multiple bursts. It can be called using the following syntax: +``` +python -m hyp3_isce2 ++process merge_tops_bursts \ + PATH_TO_UNZIPPED_PRODUCTS \ + --apply-water-mask True +``` +Where `PATH_TO_UNZIPPED_PRODUCTS` is the path to a directory containing unzipped burst InSAR products. For example: +```bash +PATH_TO_UNZIPPED_PRODUCTS/ +├─ S1_136232_IW2_20200604_20200616_VV_INT80_663F/ +├─ S1_136231_IW2_20200604_20200616_VV_INT80_529D/ +``` +In order to be merging eligible, all burst products must: +1. Have the same reference and secondary dates +1. Have the same polarization +1. Have the same multilooking +1. Be from the same relative orbit +1. Be contiguous + +The workflow should through an error if any of these conditions are not met. + +**Merging burst InSAR products requires extra data that is not contained in the production HyP3 Burst InSAR products. For the time being, to be merging eligible burst products must be created locally using your own installation of `hyp3-isce2` from the `merge_bursts` branch of this repository!** + +As mentioned above this feature is under active development, so we welcome any feedback you have! + ### Options To learn about the arguments for each workflow, look at the help documentation (`python -m hyp3_isce2 ++process [WORKFLOW_NAME] --help`). diff --git a/images/coverage.svg b/images/coverage.svg index dd6df127..6c15cace 100644 --- a/images/coverage.svg +++ b/images/coverage.svg @@ -9,13 +9,13 @@ - + coverage coverage - 61% - 61% + 75% + 75% diff --git a/pyproject.toml b/pyproject.toml index 8be3acb1..dbf2fbb4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -55,11 +55,13 @@ Documentation = "https://hyp3-docs.asf.alaska.edu" insar_tops_burst = "hyp3_isce2.insar_tops_burst:main" insar_tops = "hyp3_isce2.insar_tops:main" insar_stripmap = "hyp3_isce2.insar_stripmap:main" +merge_tops_bursts = "hyp3_isce2.merge_tops_bursts:main" [project.entry-points.hyp3] insar_tops_burst = "hyp3_isce2.insar_tops_burst:main" insar_tops = "hyp3_isce2.insar_tops:main" insar_stripmap = "hyp3_isce2.insar_stripmap:main" +merge_tops_bursts = "hyp3_isce2.merge_tops_bursts:main" [tool.pytest.ini_options] testpaths = ["tests"] @@ -77,6 +79,3 @@ readme = {file = ["README.md"], content-type = "text/markdown"} where = ["src"] [tool.setuptools_scm] - -[tool.ruff] -line-length = 120 diff --git a/src/hyp3_isce2/__main__.py b/src/hyp3_isce2/__main__.py index 58ff025f..9d50385e 100644 --- a/src/hyp3_isce2/__main__.py +++ b/src/hyp3_isce2/__main__.py @@ -19,7 +19,7 @@ def main(): parser = argparse.ArgumentParser(prefix_chars='+', formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument( '++process', - choices=['insar_tops_burst', 'insar_tops', 'insar_stripmap'], + choices=['insar_tops_burst', 'insar_tops', 'insar_stripmap', 'merge_tops_bursts'], default='insar_tops_burst', help='Select the HyP3 entrypoint to use', # HyP3 entrypoints are specified in `pyproject.toml` ) diff --git a/src/hyp3_isce2/burst.py b/src/hyp3_isce2/burst.py index fbfeda14..2ab7a3cd 100644 --- a/src/hyp3_isce2/burst.py +++ b/src/hyp3_isce2/burst.py @@ -5,16 +5,21 @@ import time from concurrent.futures import ThreadPoolExecutor from dataclasses import dataclass +from datetime import datetime, timedelta from pathlib import Path from secrets import token_hex from typing import Iterator, List, Optional, Tuple, Union import asf_search +import numpy as np import requests from isceobj.Sensor.TOPS.Sentinel1 import Sentinel1 +from isceobj.TopsProc.runMergeBursts import multilook from lxml import etree from shapely import geometry +from hyp3_isce2.utils import load_isce2_image, load_product, write_isce2_image_from_obj + log = logging.getLogger(__name__) @@ -32,6 +37,21 @@ class BurstParams: burst_number: int +@dataclass +class BurstPosition: + """Parameters describing the position of a burst and its valid data.""" + + n_lines: int + n_samples: int + first_valid_line: int + n_valid_lines: int + first_valid_sample: int + n_valid_samples: int + azimuth_time_interval: float + range_pixel_size: float + sensing_stop: datetime + + class BurstMetadata: """Metadata for a burst.""" @@ -325,11 +345,7 @@ def download_bursts(param_list: Iterator[BurstParams]) -> List[BurstMetadata]: return bursts -def get_product_name( - reference_scene: str, - secondary_scene: str, - pixel_spacing: int -) -> str: +def get_product_name(reference_scene: str, secondary_scene: str, pixel_spacing: int) -> str: """Get the name of the interferogram product. Args: @@ -354,16 +370,18 @@ def get_product_name( pixel_spacing = str(int(pixel_spacing)) product_id = token_hex(2).upper() - return '_'.join([ - platform, - burst_id, - image_plus_swath, - reference_date, - secondary_date, - polarization, - product_type + pixel_spacing, - product_id - ]) + return '_'.join( + [ + platform, + burst_id, + image_plus_swath, + reference_date, + secondary_date, + polarization, + product_type + pixel_spacing, + product_id, + ] + ) def get_burst_params(scene_name: str) -> BurstParams: @@ -402,16 +420,206 @@ def validate_bursts(reference_scene: str, secondary_scene: str) -> None: sec_polarization = sec_split[4] if ref_burst_id != sec_burst_id: - raise ValueError( - f'The reference and secondary burst IDs are not the same: {ref_burst_id} and {sec_burst_id}.' - ) + raise ValueError(f'The reference and secondary burst IDs are not the same: {ref_burst_id} and {sec_burst_id}.') if ref_polarization != sec_polarization: raise ValueError( f'The reference and secondary polarizations are not the same: {ref_polarization} and {sec_polarization}.' ) - if ref_polarization != "VV" and ref_polarization != "HH": - raise ValueError( - f'{ref_polarization} polarization is not currently supported, only VV and HH.' - ) + if ref_polarization != 'VV' and ref_polarization != 'HH': + raise ValueError(f'{ref_polarization} polarization is not currently supported, only VV and HH.') + + +def load_burst_position(swath_xml_path: str, burst_number: int) -> BurstPosition: + """Get the tiff resolution and position parameters for a burst. + + Args: + swath_xml_path: The path to the swath xml file. + burst_number: The burst number. + + Returns: + A BurstPosition object describing the burst. + """ + + product = load_product(swath_xml_path) + burst_props = product.bursts[burst_number] + + pos = BurstPosition( + n_lines=burst_props.numberOfLines, + n_samples=burst_props.numberOfSamples, + first_valid_line=burst_props.firstValidLine, + n_valid_lines=burst_props.numValidLines, + first_valid_sample=burst_props.firstValidSample, + n_valid_samples=burst_props.numValidSamples, + azimuth_time_interval=burst_props.azimuthTimeInterval, + range_pixel_size=burst_props.rangePixelSize, + sensing_stop=burst_props.sensingStop, + ) + return pos + + +def evenize(length: int, first_valid: int, valid_length: int, looks: int) -> Tuple[int]: + """Get dimensions for an image that are integer multiples of looks. + This applies to both the full image and the valid data region. + Works with either the image's lines or samples. + + Args: + length: The length of the image. + first_valid: The first valid pixel of the image. + valid_length: The length of the valid data region. + looks: The number of looks. + """ + # even_length must be a multiple of looks + n_remove = length % looks + even_length = length - n_remove + + # even_first_valid is the first multiple of looks after first_valid + even_first_valid = int(np.ceil(first_valid / looks)) * looks + + # account for the shift introduced by the shift of first_valid + n_first_valid_shift = even_first_valid - first_valid + new_valid_length = valid_length - n_first_valid_shift + + # even_valid_length must be a multiple of looks + n_valid_length_remove = new_valid_length % looks + even_valid_length = new_valid_length - n_valid_length_remove + + if (even_first_valid + even_valid_length) > even_length: + raise ValueError('The computed valid data region extends beyond the image bounds.') + + return even_length, even_first_valid, even_valid_length + + +def evenly_subset_position(position: BurstPosition, rg_looks, az_looks) -> BurstPosition: + """Get the parameters necessary to multilook a burst using even dimensions. + + Multilooking using the generated parameters ensures that there is a clear link + between pixels in the full resolution and multilooked pixel positions. + + Args: + position: The BurstPosition object describing the burst. + rg_looks: The number of range looks. + az_looks: The number of azimuth looks. + + Returns: + A BurstPosition object describing the burst. + """ + even_n_samples, even_first_valid_sample, even_n_valid_samples = evenize( + position.n_samples, position.first_valid_sample, position.n_valid_samples, rg_looks + ) + even_n_lines, even_first_valid_line, even_n_valid_lines = evenize( + position.n_lines, position.first_valid_line, position.n_valid_lines, az_looks + ) + n_lines_remove = position.n_lines - even_n_lines + even_sensing_stop = position.sensing_stop - timedelta(seconds=position.azimuth_time_interval * (n_lines_remove)) + + clip_position = BurstPosition( + n_lines=even_n_lines, + n_samples=even_n_samples, + first_valid_line=even_first_valid_line, + n_valid_lines=even_n_valid_lines, + first_valid_sample=even_first_valid_sample, + n_valid_samples=even_n_valid_samples, + azimuth_time_interval=position.azimuth_time_interval, + range_pixel_size=position.range_pixel_size, + sensing_stop=even_sensing_stop, + ) + return clip_position + + +def multilook_position(position: BurstPosition, rg_looks: int, az_looks: int) -> BurstPosition: + """Multilook a BurstPosition object. + + Args: + position: The BurstPosition object to multilook. + rg_looks: The number of range looks. + az_looks: The number of azimuth looks. + """ + multilook_position = BurstPosition( + n_lines=int(position.n_lines / az_looks), + n_samples=int(position.n_samples / rg_looks), + first_valid_line=int(position.first_valid_line / az_looks), + n_valid_lines=int(position.n_valid_lines / az_looks), + first_valid_sample=int(position.first_valid_sample / rg_looks), + n_valid_samples=int(position.n_valid_samples / rg_looks), + azimuth_time_interval=position.azimuth_time_interval * az_looks, + range_pixel_size=position.range_pixel_size * rg_looks, + sensing_stop=position.sensing_stop, + ) + return multilook_position + + +def safely_multilook( + in_file: str, position: BurstPosition, rg_looks: int, az_looks: int, subset_to_valid: bool = True +) -> None: + """Multilook an image, but only over a subset of the data whose dimensions are + integer divisible by range/azimuth looks. Do the same for the valid data region. + + Args: + in_file: The path to the input ISCE2-formatted image. + position: The BurstPosition object describing the burst. + rg_looks: The number of range looks. + az_looks: The number of azimuth looks. + subset_to_valid: Whether to subset the image to the valid data region specified by position. + """ + image_obj, array = load_isce2_image(in_file) + + dtype = image_obj.toNumpyDataType() + identity_value = np.identity(1, dtype=dtype) + mask = np.zeros((position.n_lines, position.n_samples), dtype=dtype) + + if subset_to_valid: + last_line = position.first_valid_line + position.n_valid_lines + last_sample = position.first_valid_sample + position.n_valid_samples + mask[position.first_valid_line: last_line, position.first_valid_sample: last_sample] = identity_value + else: + mask[:, :] = identity_value + + if len(array.shape) == 2: + array = array[: position.n_lines, : position.n_samples].copy() + array *= mask + else: + array = array[:, : position.n_lines, : position.n_samples].copy() + for band in range(array.shape[0]): + array[band, :, :] *= mask + + original_path = Path(image_obj.filename) + clip_path = str(original_path.parent / f'{original_path.stem}.clip{original_path.suffix}') + multilook_path = str(original_path.parent / f'{original_path.stem}.multilooked{original_path.suffix}') + + image_obj.setFilename(clip_path) + image_obj.setWidth(position.n_samples) + image_obj.setLength(position.n_lines) + write_isce2_image_from_obj(image_obj, array) + + multilook(clip_path, multilook_path, alks=az_looks, rlks=rg_looks) + + +def multilook_radar_merge_inputs( + swath_number: int, rg_looks: int, az_looks: int, base_dir: Optional[Path] = None +) -> BurstPosition: + """Multilook the radar datasets needed for post-generation product merging. + + Args: + swath_number: The swath number. + rg_looks: The number of range looks. + az_looks: The number of azimuth looks. + base_dir: The working directory. If not set, defaults to the current working directory. + """ + if base_dir is None: + base_dir = Path.cwd() + ifg_dir = base_dir / 'fine_interferogram' + geom_dir = base_dir / 'geom_reference' + + swath = f'IW{swath_number}' + position_params = load_burst_position(str(ifg_dir / f'{swath}.xml'), 0) + even_position_params = evenly_subset_position(position_params, rg_looks, az_looks) + safely_multilook(str(ifg_dir / swath / 'burst_01.int'), even_position_params, rg_looks, az_looks) + + for geom in ['lat_01.rdr', 'lon_01.rdr', 'los_01.rdr']: + geom_path = str(geom_dir / swath / geom) + safely_multilook(geom_path, even_position_params, rg_looks, az_looks, subset_to_valid=False) + + multilooked_params = multilook_position(even_position_params, rg_looks, az_looks) + return multilooked_params diff --git a/src/hyp3_isce2/insar_tops.py b/src/hyp3_isce2/insar_tops.py index 85124486..639ffd80 100644 --- a/src/hyp3_isce2/insar_tops.py +++ b/src/hyp3_isce2/insar_tops.py @@ -72,6 +72,7 @@ def insar_tops( orbit_directory=str(orbit_dir), aux_cal_directory=str(aux_cal_dir), dem_filename=str(dem_path), + geocode_dem_filename=str(dem_path), roi=roi, swaths=swaths, azimuth_looks=azimuth_looks, diff --git a/src/hyp3_isce2/insar_tops_burst.py b/src/hyp3_isce2/insar_tops_burst.py index 2eff970a..64d420d1 100644 --- a/src/hyp3_isce2/insar_tops_burst.py +++ b/src/hyp3_isce2/insar_tops_burst.py @@ -5,11 +5,11 @@ import os import subprocess import sys -from collections import namedtuple +from dataclasses import dataclass from datetime import datetime, timezone from pathlib import Path from shutil import copyfile, make_archive -from typing import Optional +from typing import Iterable, Optional import isce from hyp3lib.aws import upload_file_to_s3 @@ -18,24 +18,27 @@ from hyp3lib.util import string_is_true from isceobj.TopsProc.runMergeBursts import multilook from lxml import etree -from osgeo import gdal +from osgeo import gdal, gdalconst from pyproj import CRS import hyp3_isce2 import hyp3_isce2.metadata.util from hyp3_isce2 import topsapp from hyp3_isce2.burst import ( + BurstPosition, download_bursts, get_burst_params, get_isce2_burst_bbox, get_product_name, get_region_of_interest, + multilook_radar_merge_inputs, validate_bursts, ) from hyp3_isce2.dem import download_dem_for_isce2 from hyp3_isce2.logging import configure_root_logger from hyp3_isce2.s1_auxcal import download_aux_cal from hyp3_isce2.utils import ( + ParameterFile, get_esa_credentials, image_math, isce2_copy, @@ -52,6 +55,14 @@ log = logging.getLogger(__name__) +@dataclass +class ISCE2Dataset: + name: str + suffix: str + band: Iterable[int] + dtype: Optional[int] = gdalconst.GDT_Float32 + + def insar_tops_burst( reference_scene: str, secondary_scene: str, @@ -183,7 +194,7 @@ def make_readme( 'dem_pixel_spacing': '30 m', 'apply_water_mask': apply_water_mask, } - content = hyp3_isce2.metadata.util.render_template('insar_burst/readme.md.txt.j2', payload) + content = hyp3_isce2.metadata.util.render_template('insar_burst/insar_burst_readme.md.txt.j2', payload) output_file = product_dir / f'{product_name}_README.md.txt' with open(output_file, 'w') as f: @@ -197,6 +208,7 @@ def make_parameter_file( swath_number: int, azimuth_looks: int, range_looks: int, + multilook_position: BurstPosition, apply_water_mask: bool, dem_name: str = 'GLO_30', dem_resolution: int = 30, @@ -261,78 +273,114 @@ def make_parameter_file( slant_range_center = (slant_range_near + slant_range_far) / 2 s = ref_time.split('T')[1].split(':') - utc_time = (int(s[0]) * 60 + int(s[1]) * 60) + float(s[2]) - - output_strings = [ - f'Reference Granule: {reference_scene}\n', - f'Secondary Granule: {secondary_scene}\n', - f'Reference Pass Direction: {ref_orbit_direction}\n', - f'Reference Orbit Number: {ref_orbit_number}\n', - f'Secondary Pass Direction: {sec_orbit_direction}\n', - f'Secondary Orbit Number: {sec_orbit_number}\n', - f'Baseline: {baseline_perp}\n', - f'UTC time: {utc_time}\n', - f'Heading: {ref_heading}\n', - f'Spacecraft height: {SPACECRAFT_HEIGHT}\n', - f'Earth radius at nadir: {EARTH_RADIUS}\n', - f'Slant range near: {slant_range_near}\n', - f'Slant range center: {slant_range_center}\n', - f'Slant range far: {slant_range_far}\n', - f'Range looks: {range_looks}\n', - f'Azimuth looks: {azimuth_looks}\n', - 'INSAR phase filter: yes\n', - f'Phase filter parameter: {phase_filter_strength}\n', - 'Range bandpass filter: no\n', - 'Azimuth bandpass filter: no\n', - f'DEM source: {dem_name}\n', - f'DEM resolution (m): {dem_resolution}\n', - f'Unwrapping type: {unwrapper_type}\n', - 'Speckle filter: yes\n', - f'Water mask: {apply_water_mask}\n', - ] + utc_time = ((int(s[0]) * 60 + int(s[1])) * 60) + float(s[2]) + + parameter_file = ParameterFile( + reference_granule=reference_scene, + secondary_granule=secondary_scene, + reference_orbit_direction=ref_orbit_direction, + reference_orbit_number=ref_orbit_number, + secondary_orbit_direction=sec_orbit_direction, + secondary_orbit_number=sec_orbit_number, + baseline=float(baseline_perp), + utc_time=utc_time, + heading=ref_heading, + spacecraft_height=SPACECRAFT_HEIGHT, + earth_radius_at_nadir=EARTH_RADIUS, + slant_range_near=slant_range_near, + slant_range_center=slant_range_center, + slant_range_far=slant_range_far, + range_looks=int(range_looks), + azimuth_looks=int(azimuth_looks), + insar_phase_filter=True, + phase_filter_parameter=float(phase_filter_strength), + range_bandpass_filter=False, + azimuth_bandpass_filter=False, + dem_source=dem_name, + dem_resolution=dem_resolution, + unwrapping_type=unwrapper_type, + speckle_filter=True, + water_mask=apply_water_mask, + radar_n_lines=multilook_position.n_lines, + radar_n_samples=multilook_position.n_samples, + radar_first_valid_line=multilook_position.first_valid_line, + radar_n_valid_lines=multilook_position.n_valid_lines, + radar_first_valid_sample=multilook_position.first_valid_sample, + radar_n_valid_samples=multilook_position.n_valid_samples, + multilook_azimuth_time_interval=multilook_position.azimuth_time_interval, + multilook_range_pixel_size=multilook_position.range_pixel_size, + radar_sensing_stop=multilook_position.sensing_stop, + ) + parameter_file.write(out_path) - output_string = ''.join(output_strings) - with open(out_path.__str__(), 'w') as outfile: - outfile.write(output_string) +def find_product(pattern: str) -> str: + """Find a single file within the working directory's structure + + Args: + pattern: Glob pattern for file + + Returns + Path to file + """ + search = Path.cwd().glob(pattern) + product = str(list(search)[0]) + return product -def translate_outputs(isce_output_dir: Path, product_name: str, pixel_size: float) -> None: - """Translate ISCE outputs to a standard GTiff format with a UTM projection +def translate_outputs(product_name: str, pixel_size: float, include_radar: bool = False, use_multilooked=False) -> None: + """Translate ISCE outputs to a standard GTiff format with a UTM projection. + Assume you are in the top level of an ISCE run directory Args: - isce_output_dir: Path to the ISCE output directory product_name: Name of the product pixel_size: Pixel size + include_radar: Flag to include the full resolution radar geometry products in the output """ - src_ds = gdal.Open(str(isce_output_dir / 'filt_topophase.unw.geo')) + src_ds = gdal.Open('merged/filt_topophase.unw.geo') src_geotransform = src_ds.GetGeoTransform() src_projection = src_ds.GetProjection() - target_ds = gdal.Open(str(isce_output_dir / 'dem.crop'), gdal.GA_Update) + target_ds = gdal.Open('merged/dem.crop', gdal.GA_Update) target_ds.SetGeoTransform(src_geotransform) target_ds.SetProjection(src_projection) del src_ds, target_ds - ISCE2Dataset = namedtuple('ISCE2Dataset', ['name', 'suffix', 'band']) datasets = [ - ISCE2Dataset('filt_topophase.unw.geo', 'unw_phase', 2), - ISCE2Dataset('phsig.cor.geo', 'corr', 1), - ISCE2Dataset('dem.crop', 'dem', 1), - ISCE2Dataset('filt_topophase.unw.conncomp.geo', 'conncomp', 1), + ISCE2Dataset('merged/filt_topophase.unw.geo', 'unw_phase', [2]), + ISCE2Dataset('merged/phsig.cor.geo', 'corr', [1]), + ISCE2Dataset('merged/dem.crop', 'dem', [1]), + ISCE2Dataset('merged/filt_topophase.unw.conncomp.geo', 'conncomp', [1]), ] + suffix = '01' + if use_multilooked: + suffix += '.multilooked' + + rdr_datasets = [ + ISCE2Dataset( + find_product(f'fine_interferogram/IW*/burst_{suffix}.int.vrt'), + 'wrapped_phase_rdr', + [1], + gdalconst.GDT_CFloat32, + ), + ISCE2Dataset(find_product(f'geom_reference/IW*/lat_{suffix}.rdr.vrt'), 'lat_rdr', [1]), + ISCE2Dataset(find_product(f'geom_reference/IW*/lon_{suffix}.rdr.vrt'), 'lon_rdr', [1]), + ISCE2Dataset(find_product(f'geom_reference/IW*/los_{suffix}.rdr.vrt'), 'los_rdr', [1, 2]), + ] + if include_radar: + datasets += rdr_datasets + for dataset in datasets: out_file = str(Path(product_name) / f'{product_name}_{dataset.suffix}.tif') - in_file = str(isce_output_dir / dataset.name) - gdal.Translate( destName=out_file, - srcDS=in_file, - bandList=[dataset.band], + srcDS=dataset.name, + bandList=dataset.band, format='GTiff', + outputType=dataset.dtype, noData=0, creationOptions=['TILED=YES', 'COMPRESS=LZW', 'NUM_THREADS=ALL_CPUS'], ) @@ -342,13 +390,13 @@ def translate_outputs(isce_output_dir: Path, product_name: str, pixel_size: floa cmd = ( 'gdal_calc.py ' f'--outfile {product_name}/{product_name}_{wrapped_phase.suffix}.tif ' - f'-A {isce_output_dir / wrapped_phase.name} --A_band={wrapped_phase.band} ' + f'-A merged/{wrapped_phase.name} --A_band={wrapped_phase.band} ' '--calc angle(A) --type Float32 --format GTiff --NoDataValue=0 ' '--creation-option TILED=YES --creation-option COMPRESS=LZW --creation-option NUM_THREADS=ALL_CPUS' ) subprocess.run(cmd.split(' '), check=True) - ds = gdal.Open(str(isce_output_dir / 'los.rdr.geo'), gdal.GA_Update) + ds = gdal.Open('merged/los.rdr.geo', gdal.GA_Update) ds.GetRasterBand(1).SetNoDataValue(0) ds.GetRasterBand(2).SetNoDataValue(0) del ds @@ -361,7 +409,7 @@ def translate_outputs(isce_output_dir: Path, product_name: str, pixel_size: floa cmd = ( 'gdal_calc.py ' f'--outfile {product_name}/{product_name}_{incidence_angle.suffix}.tif ' - f'-A {isce_output_dir / incidence_angle.name} --A_band={incidence_angle.band} ' + f'-A merged/{incidence_angle.name} --A_band={incidence_angle.band} ' '--calc (90-A)*pi/180 --type Float32 --format GTiff --NoDataValue=0 ' '--creation-option TILED=YES --creation-option COMPRESS=LZW --creation-option NUM_THREADS=ALL_CPUS' ) @@ -375,18 +423,18 @@ def translate_outputs(isce_output_dir: Path, product_name: str, pixel_size: floa cmd = ( 'gdal_calc.py ' f'--outfile {product_name}/{product_name}_{azimuth_angle.suffix}.tif ' - f'-A {isce_output_dir / azimuth_angle.name} --A_band={azimuth_angle.band} ' + f'-A merged/{azimuth_angle.name} --A_band={azimuth_angle.band} ' '--calc (90+A)*pi/180 --type Float32 --format GTiff --NoDataValue=0 ' '--creation-option TILED=YES --creation-option COMPRESS=LZW --creation-option NUM_THREADS=ALL_CPUS' ) subprocess.run(cmd.split(' '), check=True) - ds = gdal.Open(str(isce_output_dir / 'filt_topophase.unw.geo')) + ds = gdal.Open('merged/filt_topophase.unw.geo') geotransform = ds.GetGeoTransform() del ds epsg = utm_from_lon_lat(geotransform[0], geotransform[3]) - files = [str(path) for path in Path(product_name).glob('*.tif')] + files = [str(path) for path in Path(product_name).glob('*.tif') if not path.name.endswith('rdr.tif')] for file in files: gdal.Warp( file, @@ -436,8 +484,8 @@ def convert_raster_from_isce2_gdal(input_image, ref_image, output_image): outputBounds=[minx, miny, maxx, maxy], xRes=pixel_size, yRes=pixel_size, - targetAlignedPixels=True - ) + targetAlignedPixels=True, + ) def main(): @@ -446,8 +494,8 @@ def main(): parser.add_argument('--bucket', help='AWS S3 bucket HyP3 for upload the final product(s)') parser.add_argument('--bucket-prefix', default='', help='Add a bucket prefix to product(s)') - parser.add_argument('--esa-username', default=None, help="Username for ESA\'s Copernicus Data Space Ecosystem") - parser.add_argument('--esa-password', default=None, help="Password for ESA\'s Copernicus Data Space Ecosystem") + parser.add_argument('--esa-username', default=None, help="Username for ESA's Copernicus Data Space Ecosystem") + parser.add_argument('--esa-password', default=None, help="Password for ESA's Copernicus Data Space Ecosystem") parser.add_argument( '--looks', choices=['20x4', '10x2', '5x1'], default='20x4', help='Number of looks to take in range and azimuth' ) @@ -479,7 +527,7 @@ def main(): range_looks, azimuth_looks = [int(looks) for looks in args.looks.split('x')] apply_water_mask = args.apply_water_mask - isce_output_dir = insar_tops_burst( + insar_tops_burst( reference_scene=reference_scene, secondary_scene=secondary_scene, azimuth_looks=azimuth_looks, @@ -491,13 +539,16 @@ def main(): ) log.info('ISCE2 TopsApp run completed successfully') + + multilook_position = multilook_radar_merge_inputs(swath_number, rg_looks=range_looks, az_looks=azimuth_looks) + pixel_size = get_pixel_size(args.looks) product_name = get_product_name(reference_scene, secondary_scene, pixel_spacing=int(pixel_size)) product_dir = Path(product_name) product_dir.mkdir(parents=True, exist_ok=True) - translate_outputs(isce_output_dir, product_name, pixel_size=pixel_size) + translate_outputs(product_name, pixel_size=pixel_size, include_radar=True, use_multilooked=True) unwrapped_phase = f'{product_name}/{product_name}_unw_phase.tif' water_mask = f'{product_name}/{product_name}_water_mask.tif' @@ -533,6 +584,7 @@ def main(): azimuth_looks=azimuth_looks, range_looks=range_looks, swath_number=swath_number, + multilook_position=multilook_position, apply_water_mask=apply_water_mask, ) output_zip = make_archive(base_name=product_name, format='zip', base_dir=product_name) diff --git a/src/hyp3_isce2/merge_tops_bursts.py b/src/hyp3_isce2/merge_tops_bursts.py new file mode 100644 index 00000000..ae71429c --- /dev/null +++ b/src/hyp3_isce2/merge_tops_bursts.py @@ -0,0 +1,1253 @@ +"""A workflow for merging standard burst InSAR products.""" +import argparse +import copy +import datetime +import logging +import os +import shutil +import sys +from concurrent.futures import ThreadPoolExecutor +from dataclasses import dataclass, field +from itertools import combinations +from pathlib import Path +from secrets import token_hex +from shutil import make_archive +from tempfile import TemporaryDirectory +from typing import Iterable, List, Optional, Tuple + +import asf_search +import isce +import isceobj +import lxml.etree as ET +import numpy as np +from contrib.Snaphu.Snaphu import Snaphu +from hyp3lib.util import string_is_true +from isceobj.Orbit.Orbit import Orbit +from isceobj.Planet.Planet import Planet +from isceobj.Sensor.TOPS.Sentinel1 import Sentinel1 +from isceobj.TopsProc.runIon import maskUnwrap +from isceobj.TopsProc.runMergeBursts import mergeBox, mergeBursts2 +from iscesys.Component import createTraitSeq +from iscesys.Component.ProductManager import ProductManager +from mroipac.filter.Filter import Filter +from mroipac.icu.Icu import Icu +from osgeo import gdal +from shapely import geometry +from stdproc.rectify.geocode.Geocodable import Geocodable +from zerodop.geozero import createGeozero + +import hyp3_isce2 +import hyp3_isce2.burst as burst_utils +from hyp3_isce2.dem import download_dem_for_isce2 +from hyp3_isce2.utils import ( + ParameterFile, + create_image, + image_math, + load_product, + make_browse_image, + read_product_metadata, + resample_to_radar_io, +) +from hyp3_isce2.water_mask import create_water_mask + + +log_format = '%(asctime)s - %(name)s - %(levelname)s - %(message)s' +logging.basicConfig(stream=sys.stdout, format=log_format, level=logging.INFO, force=True) +log = logging.getLogger(__name__) + +from hyp3_isce2.insar_tops_burst import get_pixel_size, translate_outputs # noqa + + +BURST_IFG_DIR = 'fine_interferogram' +BURST_GEOM_DIR = 'geom_reference' +FILT_WRP_IFG_NAME = 'filt_topophase.flat' +UNW_IFG_NAME = 'filt_topophase.unw' +WRP_IFG_NAME = 'topophase.flat' +COH_NAME = 'phsig.cor' +LOS_NAME = 'los.rdr' +LAT_NAME = 'lat.rdr' +LON_NAME = 'lon.rdr' +CCOM_NAME = 'filt_topophase.unw.conncomp' +GEOCODE_LIST = [UNW_IFG_NAME, COH_NAME, LOS_NAME, FILT_WRP_IFG_NAME, CCOM_NAME] + + +@dataclass +class BurstProduct: + """A dataclass to hold burst metadata""" + + granule: str + reference_date: datetime.datetime + secondary_date: datetime.datetime + burst_id: int + swath: str + polarization: str + burst_number: int + product_path: Path + n_lines: int + n_samples: int + range_looks: int + azimuth_looks: int + first_valid_line: int + n_valid_lines: int + first_valid_sample: int + n_valid_samples: int + az_time_interval: float + rg_pixel_size: float + start_utc: datetime.datetime + stop_utc: datetime.datetime + relative_orbit: int + isce2_burst_number: int = field(default=None) + + def to_burst_params(self): + """Convert to a burst_utils.BurstParams object""" + return burst_utils.BurstParams(self.granule, self.swath, self.polarization, self.burst_number) + + +def get_burst_metadata(product_paths: Iterable[Path]) -> Iterable[BurstProduct]: + """Create a list of BurstProduct objects from a set of burst product paths + + Args: + product_paths: A list of paths to unzipped burst product directories + + Returns: + A list of BurstProduct objects representing the burst metadata + """ + meta_file_paths = [path / f'{path.name}.txt' for path in product_paths] + metas = [read_product_metadata(path) for path in meta_file_paths] + + # TODO why does asf_search not return values in order? + results = [asf_search.granule_search(item['ReferenceGranule'])[0] for item in metas] + + relative_orbits = [result.properties['pathNumber'] for result in results] + granules = [Path(result.properties['url']).parts[2] for result in results] + pattern = '%Y%m%dT%H%M%S' + reference_granules = [datetime.datetime.strptime(item['ReferenceGranule'].split('_')[3], pattern) for item in metas] + secondary_granules = [datetime.datetime.strptime(item['SecondaryGranule'].split('_')[3], pattern) for item in metas] + swaths = [result.properties['burst']['subswath'] for result in results] + burst_ids = [result.properties['burst']['relativeBurstID'] for result in results] + burst_indexes = [result.properties['burst']['burstIndex'] for result in results] + polarization = [result.properties['polarization'] for result in results] + start_utc = [ + datetime.datetime.strptime(result.properties['startTime'], '%Y-%m-%dT%H:%M:%S.%fZ') for result in results + ] + relative_orbits = [result.properties['pathNumber'] for result in results] + n_lines = [int(meta['Radarnlines']) for meta in metas] + n_samples = [int(meta['Radarnsamples']) for meta in metas] + range_looks = [int(meta['Rangelooks']) for meta in metas] + azimuth_looks = [int(meta['Azimuthlooks']) for meta in metas] + first_valid_line = [int(meta['Radarfirstvalidline']) for meta in metas] + n_valid_lines = [int(meta['Radarnvalidlines']) for meta in metas] + first_valid_sample = [int(meta['Radarfirstvalidsample']) for meta in metas] + n_valid_samples = [int(meta['Radarnvalidsamples']) for meta in metas] + az_time_interval = [float(meta['Multilookazimuthtimeinterval']) for meta in metas] + rg_pixel_size = [float(meta['Multilookrangepixelsize']) for meta in metas] + stop_utc = [datetime.datetime.strptime(meta['Radarsensingstop'], '%Y-%m-%dT%H:%M:%S.%f') for meta in metas] + products = [] + for i in range(len(granules)): + product = BurstProduct( + granules[i], + reference_granules[i], + secondary_granules[i], + burst_ids[i], + swaths[i], + polarization[i], + burst_indexes[i], + product_paths[i], + n_lines[i], + n_samples[i], + range_looks[i], + azimuth_looks[i], + first_valid_line[i], + n_valid_lines[i], + first_valid_sample[i], + n_valid_samples[i], + az_time_interval[i], + rg_pixel_size[i], + start_utc[i], + stop_utc[i], + relative_orbits[i], + ) + products.append(product) + return products + + +def prep_metadata_dirs(base_path: Optional[Path] = None) -> Tuple[Path, Path]: + """Create the annotation and manifest directories + + Args: + base_path: The base path to create the directories in. Defaults to the current working directory. + + Returns: + A tuple of the annotation and manifest directories + """ + if base_path is None: + base_path = Path.cwd() + annotation_dir = base_path / 'annotation' + annotation_dir.mkdir(exist_ok=True) + manifest_dir = base_path / 'manifest' + manifest_dir.mkdir(exist_ok=True) + return annotation_dir, manifest_dir + + +def download_metadata_xmls(params: Iterable[burst_utils.BurstParams], base_dir: Optional[Path] = None) -> None: + """Download annotation xmls for a set of burst parameters to a directory + named 'annotation' in the current working directory + + Args: + params: A list of burst_utils.BurstParams objects + base_dir: The base directory to download the metadata to. Defaults to the current working directory. + """ + if base_dir is None: + base_dir = Path.cwd() + annotation_dir, manifest_dir = prep_metadata_dirs(base_dir) + + granules = list(set([param.granule for param in params])) + download_params = [burst_utils.BurstParams(granule, 'IW1', 'VV', 0) for granule in granules] + + with burst_utils.get_asf_session() as asf_session: + with ThreadPoolExecutor(max_workers=10) as executor: + xml_futures = [ + executor.submit(burst_utils.download_metadata, asf_session, param) for param in download_params + ] + metadata_xmls = {param.granule: future.result() for param, future in zip(download_params, xml_futures)} + + et_args = {'encoding': 'UTF-8', 'xml_declaration': True} + for param in params: + metadata_xml = metadata_xmls[param.granule] + burst_metadata = burst_utils.BurstMetadata(metadata_xml, param) + ET.ElementTree(burst_metadata.annotation).write(annotation_dir / burst_metadata.annotation_name, **et_args) + ET.ElementTree(burst_metadata.manifest).write(manifest_dir / f'{burst_metadata.safe_name}.xml', **et_args) + + +def get_scene_roi(s1_obj_bursts: Iterable[isceobj.Sensor.TOPS.BurstSLC.BurstSLC]) -> Tuple: + """Get the bounding box of a set of ISCE2 Sentinel1 burst objects + + Args: + s1_obj_bursts: A list of ISCE2 Sentinel1 burst objects + + Returns: + A tuple of (west, south, east, north) bounding box coordinates + """ + for i, isce_burst in enumerate(s1_obj_bursts): + snwe = isce_burst.getBbox() + bbox = geometry.box(snwe[2], snwe[0], snwe[3], snwe[1]) + if i == 0: + overall_bbox = bbox + else: + overall_bbox = overall_bbox.union(bbox) + return overall_bbox.bounds + + +class Sentinel1BurstSelect(Sentinel1): + """A Modified version of the ISCE2 Sentinel1 class that allows for subsetting of bursts""" + + def select_bursts(self, start_utcs: Iterable[datetime.datetime]) -> None: + """Subset the burst list to only include bursts with start times in start_utcs + + Args: + start_utcs: A list of burst start times to subset the burst list to + """ + start_utcs = sorted(start_utcs) + cropList = createTraitSeq('burst') + tiffList = [] + eapList = [] + + print('Number of Bursts before cropping: ', len(self.product.bursts)) + for start_utc in start_utcs: + start_utc = start_utc.replace(microsecond=0) + match = [ + x for x in enumerate(self.product.bursts) if x[1].burstStartUTC.replace(microsecond=0) == start_utc + ] + + if not len(match) > 0: + raise ValueError(f'No match for burst at time {start_utc} found.') + ind, isce2_burst = match[0] + + cropList.append(isce2_burst) + if len(self._tiffSrc): + tiffList.append(self._tiffSrc[ind]) + eapList.append(self._elevationAngleVsTau[ind]) + + # Actual cropping + self.product.bursts = cropList + self.product.numberOfBursts = len(self.product.bursts) + + self._tiffSrc = tiffList + self._elevationAngleVsTau = eapList + print('Number of Bursts after cropping: ', len(self.product.bursts)) + + def update_burst_properties(self, products: Iterable[BurstProduct]) -> None: + """Update burst properties based on the burst metadata and previous subset operations + + Args: + products: A list of BurstProduct objects + """ + width = self._burstWidth + length = self._burstLength + if width is None: + width = self.product.bursts[0].numberOfSamples + if length is None: + length = self.product.bursts[0].numberOfLines + + for index, burst in enumerate(self.product.bursts): + product = products[index] + if product.start_utc.replace(microsecond=0) != burst.burstStartUTC.replace(microsecond=0): + raise ValueError('Burst product and ISCE2 burst do not match (different start times).') + + burst.firstValidLine = product.first_valid_line + burst.numValidLines = product.n_valid_lines + burst.firstValidSample = product.first_valid_sample + burst.numValidSamples = product.n_valid_samples + + outfile = os.path.join(self.output, 'burst_%02d' % (index + 1) + '.slc') + slcImage = isceobj.createSlcImage() + slcImage.setByteOrder('l') + slcImage.setFilename(outfile) + slcImage.setAccessMode('read') + slcImage.setWidth(width) + slcImage.setLength(length) + slcImage.setXmin(0) + slcImage.setXmax(width) + burst.image = slcImage + burst.numberOfSamples = width + burst.numberOfLines = length + + print('Updating burst number from {0} to {1}'.format(burst.burstNumber, index + 1)) + burst.burstNumber = index + 1 + + def write_xml(self) -> None: + """Write the product xml to the directory specified by self.output""" + pm = ProductManager() + pm.configure() + + outxml = self.output + if outxml.endswith('/'): + outxml = outxml[:-1] + pm.dumpProduct(self.product, os.path.join(outxml + '.xml')) + + +def load_isce_s1_obj(swath: int, polarization: str, base_dir: Optional[Path] = None) -> Sentinel1BurstSelect: + """Load a modified ISCE2 Sentinel1 instance for a swath and polarization given annotation and manifest directories + + Args: + swath: The swath number + polarization: The polarization + base_dir: The base directory containing the annotation and manifest directories. Defaults to the woking dir. + + Returns: + A modified ISCE2 Sentinel1 instance + """ + if base_dir is None: + base_dir = Path.cwd() + base_dir = Path(base_dir) + + annotation_dir = base_dir / 'annotation' + manifest_dir = base_dir / 'manifest' + + annotation_xmls = [str(path) for path in annotation_dir.glob(f's1?-??{swath}-slc-{polarization.lower()}*')] + if len(annotation_xmls) == 0: + raise ValueError( + f'No annotation files for swath {swath} and polarization {polarization} found in annotation directory' + ) + manifest_xmls = [str(path) for path in manifest_dir.glob('S1*.xml')] + + obj = Sentinel1BurstSelect() + obj.configure() + obj.xml = annotation_xmls + obj.tiff = ['' for _ in range(len(annotation_xmls))] + obj.manifest = manifest_xmls + obj.swath = swath + obj.polarization = polarization.lower() + obj.parse() + return obj + + +def create_burst_cropped_s1_obj( + swath: str, products: Iterable[BurstProduct], polarization: str = 'VV', base_dir: Optional[Path] = None +) -> Sentinel1BurstSelect: + """Create an ISCE2 Sentinel1BurstSelect instance for a set of burst products, and write the xml. + Also updates the BurstProduct objects with the ISCE2 burst number. + + Args: + swath: The swath id (e.g., IW1) of the burst products + products: A list of BurstProduct objects to create the ISCE2 Sentinel1 instance for + polarization: The polarization of the burst products + base_dir: The base directory to write the xml to. Defaults to the current working directory. + + Returns: + A tuple of the updated BurstProduct objects and the ISCE2 Sentinel1 instance + """ + if base_dir is None: + base_dir = Path.cwd() + + swaths_in_products = list(set([int(product.swath[2:3]) for product in products])) + if len(swaths_in_products) > 1 or swaths_in_products[0] != swath: + raise ValueError(f'Products provided are not all in swath {swath}') + + obj = load_isce_s1_obj(swath, polarization, base_dir=base_dir) + + out_path = base_dir / BURST_IFG_DIR + out_path.mkdir(exist_ok=True) + obj.output = str(out_path / f'IW{swath}') + + products = sorted(products, key=lambda x: x.start_utc) + obj.select_bursts([b.start_utc for b in products]) + obj.update_burst_properties(products) + obj.write_xml() + return obj + + +def modify_for_multilook( + burst_products: Iterable[BurstProduct], swath_obj: Sentinel1BurstSelect, outdir: str = BURST_IFG_DIR +) -> None: + """Modify a Sentinel1 instance so that it is compatible with previously multilooked burst products + + Args: + burst_products: A list of BurstProduct objects containing the needed metadata + swath_obj: A Sentinel1BurstSelect (or Sentinel1) instance representing the parent swath + outdir: The directory to write the xml to + """ + multilook_swath_obj = copy.deepcopy(swath_obj) + multilook_swath_obj.output = os.path.join(outdir, 'IW{0}_multilooked'.format(multilook_swath_obj.swath)) + multilook_swath_obj.output = str(Path(outdir) / f'IW{multilook_swath_obj.swath}_multilooked') + for new_metadata, burst_obj in zip(burst_products, multilook_swath_obj.product.bursts): + if new_metadata.start_utc.replace(microsecond=0) != burst_obj.burstStartUTC.replace(microsecond=0): + raise ValueError('Burst product and ISCE2 burst do not match (different start times).') + burst_obj.numberOfSamples = new_metadata.n_samples + burst_obj.numberOfLines = new_metadata.n_lines + burst_obj.firstValidSample = new_metadata.first_valid_sample + burst_obj.numValidSamples = new_metadata.n_valid_samples + burst_obj.firstValidLine = new_metadata.first_valid_line + burst_obj.numValidLines = new_metadata.n_valid_lines + burst_obj.sensingStop = new_metadata.stop_utc + burst_obj.azimuthTimeInterval = new_metadata.az_time_interval + burst_obj.rangePixelSize = new_metadata.rg_pixel_size + + return multilook_swath_obj + + +def download_dem_for_multiple_bursts(s1_objs: Iterable[Sentinel1BurstSelect], base_dir=None): + """Download the DEM for the region covered in a set of ISCE2 Sentinel1 instances + + Args: + s1_objs: A list of Sentinel1BurstSelect instances + base_dir: The base directory to download the DEM to. Defaults to the current working directory. + """ + if base_dir is None: + base_dir = Path.cwd() + base_dir = Path(base_dir) + + burst_objs = [] + for s1_obj in s1_objs: + burst_objs += s1_obj.product.bursts + dem_roi = get_scene_roi(burst_objs) + download_dem_for_isce2(dem_roi, dem_name='glo_30', dem_dir=base_dir, buffer=0, resample_20m=False) + + +def translate_image(in_path: str, out_path: str, image_type: str) -> None: + """Translate a HyP3 burst product image to an ISCE2 compatible image + + Args: + in_path: The path to the input image + out_path: The path to the output image + image_type: The type of image to translate can be one of 'ifg', 'los', 'lat', or 'lon' + """ + info = gdal.Info(in_path, format='json') + width = info['size'][0] + if image_type in 'ifg': + out_img = isceobj.createIntImage() + n_bands = 1 + out_img.initImage(out_path, 'read', width, bands=n_bands) + elif image_type in ['lat', 'lon']: + n_bands = 1 + out_img = isceobj.createImage() + out_img.initImage(out_path, 'read', width, 'DOUBLE', bands=n_bands) + elif image_type == 'los': + out_img = isceobj.createImage() + n_bands = 2 + out_img.initImage(out_path, 'read', width, 'DOUBLE', bands=n_bands, scheme='bil') + out_img.imageType = 'bil' + else: + raise NotImplementedError(f'{image_type} is not a valid format') + + with TemporaryDirectory() as tmpdir: + out_tmp_path = str(Path(tmpdir) / Path(out_path).name) + gdal.Translate( + out_tmp_path, + in_path, + bandList=[n + 1 for n in range(n_bands)], + format='ENVI', + creationOptions=['INTERLEAVE=BIL'], + ) + shutil.copy(out_tmp_path, out_path) + out_img.renderHdr() + + +def spoof_isce2_setup( + burst_products: Iterable[BurstProduct], s1_obj: Sentinel1BurstSelect, base_dir: Optional[Path] = None +) -> None: + """For a set of ASF burst products, create spoofed geom_reference and fine_interferogram directories + that are in the state they would be in after running topsApp.py from the 'startup' step to the 'burstifg' step. + + Args: + burst_products: A list of BurstProduct objects + s1_obj: An ISCE2 Sentinel1 instance representing the parent swath + base_dir: The base directory to write the spoofed directories to. Defaults to the current working directory. + """ + if base_dir is None: + base_dir = Path.cwd() + base_dir = Path(base_dir) + + ifg_dir = base_dir / 'fine_interferogram' + ifg_dir.mkdir(exist_ok=True) + + geom_dir = base_dir / 'geom_reference' + geom_dir.mkdir(exist_ok=True) + + swaths = list(set([product.swath for product in burst_products])) + for swath in swaths: + ifg_swath_path = ifg_dir / swath + ifg_swath_path.mkdir(exist_ok=True) + + geom_swath_path = geom_dir / swath + geom_swath_path.mkdir(exist_ok=True) + + file_types = { + 'ifg': 'wrapped_phase_rdr', + 'los': 'los_rdr', + 'lat': 'lat_rdr', + 'lon': 'lon_rdr', + } + for product in burst_products: + for image_type in file_types: + if image_type == 'ifg': + img_dir = ifg_dir + name = f'burst_{product.isce2_burst_number:02}.int' + else: + img_dir = geom_dir + name = f'{image_type}_{product.isce2_burst_number:02}.rdr' + in_path = str(product.product_path / f'{product.product_path.stem}_{file_types[image_type]}.tif') + out_path = str(img_dir / product.swath / name) + # translate_image(in_path, out_path, product.n_samples, image_type) + translate_image(in_path, out_path, image_type) + + +def get_swath_list(base_dir: Path) -> list[str]: + """Get the list of swaths from a directory of burst products + + Args: + base_dir: The directory containing the burst products + + Returns: + A list of swaths + """ + swathList = [] + for x in [1, 2, 3]: + swath_path = Path(base_dir) / f'IW{x}' + if swath_path.exists(): + swathList.append(x) + + return swathList + + +def get_merged_orbit(products: Iterable[Sentinel1]) -> Orbit: + """Create a merged orbit from a set of ISCE2 Sentinel1 products + + Args: + product: A list of ISCE2 Sentinel1 products + + Returns: + The merged orbit + """ + # Create merged orbit + orb = Orbit() + orb.configure() + + burst = products[0].bursts[0] + # Add first burst orbit to begin with + for sv in burst.orbit: + orb.addStateVector(sv) + + for pp in products: + # Add all state vectors + for bb in pp.bursts: + for sv in bb.orbit: + if (sv.time < orb.minTime) or (sv.time > orb.maxTime): + orb.addStateVector(sv) + + bb.orbit = orb + + return orb + + +def get_frames_and_indexes(burst_ifg_dir: Path) -> Tuple: + """Get the frames and burst indexes from a directory of burst interferograms. + + Args: + burst_ifg_dir: The directory containing the burst interferograms + + Returns: + A tuple of the frames and burst indexes + """ + frames = [] + burst_index = [] + + swath_list = get_swath_list(burst_ifg_dir) + for swath in swath_list: + ifg = load_product(os.path.join(burst_ifg_dir, 'IW{0}_multilooked.xml'.format(swath))) + min_burst = ifg.bursts[0].burstNumber - 1 + max_burst = ifg.bursts[-1].burstNumber + frames.append(ifg) + burst_index.append([int(swath), min_burst, max_burst]) + + return frames, burst_index + + +def merge_bursts(range_looks: int, azimuth_looks: int, merge_dir: str = 'merged') -> None: + """Merge burst products into a multi-swath product, and multilook. + Note: Can't test this function without polluting home directory with ISCE2 files since + mergeBursts2 sets up pathing assuming it is in the working directory. + + Args: + azimuth_looks: The number of azimuth looks + range_looks: The number of range looks + merge_dir: The directory to write the merged product to + """ + frames, burstIndex = get_frames_and_indexes(BURST_IFG_DIR) + + box = mergeBox(frames) + file_types = { + 'ifg': (BURST_IFG_DIR, 'burst_%02d.int', WRP_IFG_NAME), + 'los': (BURST_GEOM_DIR, 'los_%02d.rdr', LOS_NAME), + 'lat': (BURST_GEOM_DIR, 'lat_%02d.rdr', LAT_NAME), + 'lon': (BURST_GEOM_DIR, 'lon_%02d.rdr', LON_NAME), + } + + for file_type in file_types: + directory, file_pattern, name = file_types[file_type] + burst_paths = os.path.join(directory, 'IW%d', file_pattern) + out_path = os.path.join(merge_dir, name) + mergeBursts2(frames, burst_paths, burstIndex, box, out_path, virtual=True, validOnly=True) + with TemporaryDirectory() as tmpdir: + out_tmp_path = str(Path(tmpdir) / Path(out_path).name) + gdal.Translate(out_tmp_path, out_path + '.vrt', format='ENVI', creationOptions=['INTERLEAVE=BIL']) + shutil.copy(out_tmp_path, out_path) + + +def goldstein_werner_filter(in_path: Path, out_path: Path, coh_path: Path, filter_strength: float = 0.6) -> None: + """Apply the Goldstein-Werner filter to the merged interferogram. + See https://doi.org/10.1029/1998GL900033 for method details. + + Args: + in_path: The path to the input interferogram + out_path: The path to the output filtered interferogram + coh_path: The path to the coherence file + filter_strength: The filter strength (0,1] + """ + ifg_image = create_image(str(in_path), image_subtype='ifg', action='load') + width_ifg = ifg_image.getWidth() + + out_path = str(out_path) + filt_image = create_image(out_path, width_ifg, 'write', image_subtype='ifg') + + filter_obj = Filter() + filter_obj.wireInputPort(name='interferogram', object=ifg_image) + filter_obj.wireOutputPort(name='filtered interferogram', object=filt_image) + filter_obj.goldsteinWerner(alpha=filter_strength) + + ifg_image.finalizeImage() + filt_image.finalizeImage() + del filt_image + + filt_image = create_image(out_path, width_ifg, 'read', image_subtype='ifg') + phsigImage = create_image(str(coh_path), width_ifg, 'write', image_subtype='cor') + + icu_obj = Icu(name='topsapp_filter_icu') + icu_obj.configure() + icu_obj.unwrappingFlag = False + icu_obj.useAmplitudeFlag = False + icu_obj.icu(intImage=filt_image, phsigImage=phsigImage) + + filt_image.finalizeImage() + phsigImage.finalizeImage() + phsigImage.renderHdr() + + +def mask_coherence(out_name: str, merge_dir: Optional[Path] = None) -> None: + """Mask the coherence image with a water mask that has been resampled to radar coordinates + + Args: + out_name: The name of the output masked coherence file + mergedir: The output directory containing the merged interferogram + """ + if merge_dir is None: + merge_dir = Path.cwd() / 'merged' + merge_dir = Path(merge_dir) + + input_files = ('water_mask', LAT_NAME, LON_NAME, 'water_mask.rdr', COH_NAME, out_name) + mask_geo, lat, lon, mask_rdr, coh, masked_coh = [str(Path(merge_dir) / name) for name in input_files] + create_water_mask(input_image='full_res.dem.wgs84', output_image=mask_geo, gdal_format='ISCE') + resample_to_radar_io(mask_geo, lat, lon, mask_rdr) + image_math(coh, mask_rdr, masked_coh, 'a*b') + + +def snaphu_unwrap( + range_looks: int, + azimuth_looks: int, + corrfile: Optional[Path] = None, + base_dir: Optional[Path] = None, + cost_mode='DEFO', + init_method='MST', + defomax=4.0, +) -> None: + """Unwrap the merged interferogram using SNAPHU + + Args: + azimuth_looks: The number of azimuth looks + range_looks: The number of range looks + corrfile: The path to the coherence file, will default to the coherence file in mergedir + mergedir: The output directory containing the merged interferogram + cost_mode: The cost mode to use for SNAPHU + init_method: The initialization method to use for SNAPHU + defomax: The maximum deformation allowed for SNAPHU + """ + if base_dir is None: + base_dir = Path.cwd() / 'merged' + base_dir = Path(base_dir) + + burst_ifg_dir = base_dir.parent / BURST_IFG_DIR + if not corrfile: + corrfile = base_dir / COH_NAME + wrap_name = base_dir / FILT_WRP_IFG_NAME + unwrap_name = base_dir / UNW_IFG_NAME + + img = isceobj.createImage() + img.load(str(wrap_name) + '.xml') + width = img.getWidth() + + swath = get_swath_list(str(burst_ifg_dir))[0] + ifg = load_product(str(burst_ifg_dir / f'IW{swath}_multilooked.xml')) + wavelength = ifg.bursts[0].radarWavelength + + # tmid + tstart = ifg.bursts[0].sensingStart + tend = ifg.bursts[-1].sensingStop + tmid = tstart + 0.5 * (tend - tstart) + + # Sometimes tmid may exceed the time span, so use mid burst instead + burst_index = int(np.around(len(ifg.bursts) / 2)) + orbit = ifg.bursts[burst_index].orbit + peg = orbit.interpolateOrbit(tmid, method='hermite') + + # Calculate parameters for SNAPHU + ref_elp = Planet(pname='Earth').ellipsoid + llh = ref_elp.xyz_to_llh(peg.getPosition()) + hdg = orbit.getENUHeading(tmid) + ref_elp.setSCH(llh[0], llh[1], hdg) + earth_radius = ref_elp.pegRadCur + altitude = llh[2] + azfact = 0.8 + rngfact = 0.8 + corrLooks = range_looks * azimuth_looks / (azfact * rngfact) + maxComponents = 20 + + snp = Snaphu() + snp.setInitOnly(False) + snp.setInput(str(wrap_name)) + snp.setOutput(str(unwrap_name)) + snp.setWidth(width) + snp.setCostMode(cost_mode) + snp.setEarthRadius(earth_radius) + snp.setWavelength(wavelength) + snp.setAltitude(altitude) + snp.setCorrfile(str(corrfile)) + snp.setInitMethod(init_method) + snp.setCorrLooks(corrLooks) + snp.setMaxComponents(maxComponents) + snp.setDefoMaxCycles(defomax) + snp.setRangeLooks(range_looks) + snp.setAzimuthLooks(azimuth_looks) + snp.setCorFileFormat('FLOAT_DATA') + snp.prepare() + snp.unwrap() + + # Render XML + create_image(str(unwrap_name), width, 'read', image_subtype='unw', action='finalize') + + # Check if connected components was created + if not snp.dumpConnectedComponents: + raise RuntimeError('SNAPHU did not create connected components file') + + create_image(str(unwrap_name) + '.conncomp', width, 'read', image_subtype='conncomp', action='finalize') + maskUnwrap(str(unwrap_name), str(wrap_name)) + + +def geocode_products( + range_looks: int, + azimuth_looks: int, + dem_path: Path, + base_dir: Optional[Path] = None, + to_be_geocoded=GEOCODE_LIST, +) -> None: + """Geocode a set of ISCE2 products + + Args: + azimuth_looks: The number of azimuth looks + range_looks: The number of range looks + dem_path: The path to the DEM + base_dir: The output directory containing the merged interferogram + to_be_geocoded: A list of products to geocode + """ + if base_dir is None: + base_dir = Path.cwd() / 'merged' + base_dir = Path(base_dir) + burst_ifg_dir = base_dir.parent / BURST_IFG_DIR + + to_be_geocoded = [str(base_dir / file) for file in to_be_geocoded] + swath_list = get_swath_list(str(burst_ifg_dir)) + + frames = [] + for swath in swath_list: + reference_product = load_product(str(burst_ifg_dir / f'IW{swath}.xml')) + frames.append(reference_product) + + orbit = get_merged_orbit(frames) + + bboxes = [] + for frame in frames: + bboxes.append(frame.getBbox()) + + snwe = [ + min([x[0] for x in bboxes]), + max([x[1] for x in bboxes]), + min([x[2] for x in bboxes]), + max([x[3] for x in bboxes]), + ] + + # Identify the 4 corners and dimensions + top_swath = min(frames, key=lambda x: x.sensingStart) + left_swath = min(frames, key=lambda x: x.startingRange) + + # Get required values from product + burst = frames[0].bursts[0] + t0 = top_swath.sensingStart + dtaz = burst.azimuthTimeInterval + r0 = left_swath.startingRange + dr = burst.rangePixelSize + wvl = burst.radarWavelength + planet = Planet(pname='Earth') + + # Setup DEM + demImage = isceobj.createDemImage() + demImage.load(dem_path + '.xml') + + # Geocode one by one + ge = Geocodable() + for prod in to_be_geocoded: + geo_obj = createGeozero() + geo_obj.configure() + + geo_obj.snwe = snwe + geo_obj.demCropFilename = os.path.join(base_dir, 'dem.crop') + geo_obj.numberRangeLooks = range_looks + geo_obj.numberAzimuthLooks = azimuth_looks + geo_obj.lookSide = -1 + + # create the instance of the input image and the appropriate geocode method + inImage, method = ge.create(prod) + geo_obj.method = method + + geo_obj.slantRangePixelSpacing = dr + geo_obj.prf = 1.0 / dtaz + geo_obj.orbit = orbit + geo_obj.width = inImage.getWidth() + geo_obj.length = inImage.getLength() + geo_obj.dopplerCentroidCoeffs = [0.0] + geo_obj.radarWavelength = wvl + + geo_obj.rangeFirstSample = r0 + ((range_looks - 1) / 2.0) * dr + geo_obj.setSensingStart(t0 + datetime.timedelta(seconds=(((azimuth_looks - 1) / 2.0) * dtaz))) + geo_obj.wireInputPort(name='dem', object=demImage) + geo_obj.wireInputPort(name='planet', object=planet) + geo_obj.wireInputPort(name='tobegeocoded', object=inImage) + + geo_obj.geocode() + + +def get_product_name(product: BurstProduct, pixel_size: int) -> str: + """Create a product name for a merged burst product. Follows the convention used by ASF burst products, + but replaces the burst id with the relative orbit number, and removes the swath compenent with ''. + + Args: + product: The BurstProduct object to create the name for. Can be any product in the merged product. + pixel_size: The pixel size of the product + + Returns: + The merged product name + """ + swath_placeholder = '' + reference_date = product.reference_date.strftime('%Y%m%d') + secondary_date = product.secondary_date.strftime('%Y%m%d') + product_id = token_hex(2).upper() + return '_'.join( + [ + 'S1', + f'{product.relative_orbit:03d}', + swath_placeholder, + reference_date, + secondary_date, + product.polarization.upper(), + f'INT{int(pixel_size)}', + product_id, + ] + ) + + +def get_product_metadata_info(base_dir: Path) -> List: + """Get the metadata for a set of ASF burst products + + Args: + base_dir: The directory containing UNZIPPED ASF burst products + + Returns: + A list of metadata dictionaries + """ + product_paths = list(Path(base_dir).glob('S1_??????_IW?_*')) + meta_file_paths = [path / f'{path.name}.txt' for path in product_paths] + metas = [read_product_metadata(path) for path in meta_file_paths] + return metas + + +def make_parameter_file( + out_path: Path, + metas: List, + range_looks: int, + azimuth_looks: int, + filter_strength: float, + water_mask: bool, + dem_name: str = 'GLO_30', + dem_resolution: int = 30, + base_dir: Optional[Path] = None, +): + """Create a parameter file for the ASF merged burst product and write it to out_path + + Args: + out_path: The path to write the parameter file to + metas: A list of metadata dictionaries for the burst products + range_looks: The number of range looks + azimuth_looks: The number of azimuth looks + filter_strength: The Goldstein-Werner filter strength + water_mask: Whether or not to use a water mask + dem_name: The name of the source DEM + dem_resolution: The resolution of the source DEM + base_dir: The base directory to write the parameter file to. Defaults to the current working directory. + """ + if base_dir is None: + base_dir = Path.cwd() + base_dir = Path(base_dir) + + SPACECRAFT_HEIGHT = 693000.0 + EARTH_RADIUS = 6337286.638938101 + + reference_scenes = [meta['ReferenceGranule'] for meta in metas] + secondary_scenes = [meta['SecondaryGranule'] for meta in metas] + ref_orbit_number = metas[0]['ReferenceOrbitNumber'] + sec_orbit_number = metas[0]['SecondaryOrbitNumber'] + baseline_perp = metas[0]['Baseline'] + + burst_ifg_dir = base_dir / BURST_IFG_DIR + insar_product = load_product(burst_ifg_dir / f'IW{get_swath_list(burst_ifg_dir)[0]}.xml') + + orbit_direction = insar_product.bursts[0].passDirection + ref_heading = insar_product.orbit.getHeading() + ref_time = insar_product.bursts[0].sensingStart + + utc_time = (((ref_time.hour * 60) + ref_time.minute) * 60) + ref_time.second + (ref_time.microsecond * 10e-7) + slant_range_near = insar_product.startingRange + slant_range_center = insar_product.midRange + slant_range_far = insar_product.farRange + + parameter_file = ParameterFile( + reference_granule=', '.join(reference_scenes), + secondary_granule=', '.join(secondary_scenes), + reference_orbit_direction=orbit_direction, + reference_orbit_number=ref_orbit_number, + secondary_orbit_direction=orbit_direction, + secondary_orbit_number=sec_orbit_number, + baseline=baseline_perp, + utc_time=utc_time, + heading=ref_heading, + spacecraft_height=SPACECRAFT_HEIGHT, + earth_radius_at_nadir=EARTH_RADIUS, + slant_range_near=slant_range_near, + slant_range_center=slant_range_center, + slant_range_far=slant_range_far, + range_looks=range_looks, + azimuth_looks=azimuth_looks, + insar_phase_filter=True, + phase_filter_parameter=filter_strength, + range_bandpass_filter=False, + azimuth_bandpass_filter=False, + dem_source=dem_name, + dem_resolution=dem_resolution, + unwrapping_type='snaphu_mcf', + speckle_filter=True, + water_mask=water_mask, + ) + parameter_file.write(out_path) + + +def make_readme( + product_dir: Path, + reference_scenes: Iterable, + secondary_scenes: Iterable, + range_looks: int, + azimuth_looks: int, + apply_water_mask: bool, +) -> None: + """Create a README file for the merged burst product and write it to product_dir + + Args: + product_dir: The path to the directory containing the merged burst product, + the directory name should be the product name + reference_scenes: A list of reference scenes + secondary_scenes: A list of secondary scenes + range_looks: The number of range looks + azimuth_looks: The number of azimuth looks + apply_water_mask: Whether or not a water mask was applied + """ + product_name = product_dir.name + wrapped_phase_path = product_dir / f'{product_name}_wrapped_phase.tif' + info = gdal.Info(str(wrapped_phase_path), format='json') + secondary_granule_datetime_str = secondary_scenes[0].split('_')[3] + + payload = { + 'processing_date': datetime.datetime.now(datetime.timezone.utc), + 'plugin_name': hyp3_isce2.__name__, + 'plugin_version': hyp3_isce2.__version__, + 'processor_name': isce.__name__.upper(), # noqa + 'processor_version': isce.__version__, # noqa + 'projection': hyp3_isce2.metadata.util.get_projection(info['coordinateSystem']['wkt']), + 'pixel_spacing': info['geoTransform'][1], + 'product_name': product_name, + 'reference_burst_name': ', '.join(reference_scenes), + 'secondary_burst_name': ', '.join(secondary_scenes), + 'range_looks': range_looks, + 'azimuth_looks': azimuth_looks, + 'secondary_granule_date': datetime.datetime.strptime(secondary_granule_datetime_str, '%Y%m%dT%H%M%S'), + 'dem_name': 'GLO-30', + 'dem_pixel_spacing': '30 m', + 'apply_water_mask': apply_water_mask, + } + content = hyp3_isce2.metadata.util.render_template('insar_burst/insar_burst_merge_readme.md.txt.j2', payload) + + output_file = product_dir / f'{product_name}_README.md.txt' + with open(output_file, 'w') as f: + f.write(content) + + +def check_burst_group_validity(products) -> None: + """Check that a set of burst products are valid for merging. This includes: + All products have the same: + - date + - relative orbit + - polarization + - multilook + All products must also be contiguous. This means: + - A given swath has a continuous series of bursts + - Neighboring swaths have at at most one burst separation + + This function will raise a ValueError if any of these conditions are not met. + + Args: + products: A list of BurstProduct objects + """ + reference_dates = set([product.reference_date.date() for product in products]) + secondary_dates = set([product.secondary_date.date() for product in products]) + polarizations = set([product.polarization for product in products]) + relative_orbits = set([product.relative_orbit for product in products]) + range_looks = [product.range_looks for product in products] + azimuth_looks = [product.azimuth_looks for product in products] + looks = set([f'{rgl}x{azl}' for rgl, azl in zip(range_looks, azimuth_looks)]) + + sets = { + 'reference_date': reference_dates, + 'secondary_date': secondary_dates, + 'polarization': polarizations, + 'relative_orbit': relative_orbits, + 'looks': looks, + } + for key, value in sets.items(): + if len(value) > 1: + key_name = key.replace('_', ' ') + value_names = ', '.join([str(v) for v in value]) + raise ValueError(f'All products must have the same {key_name}. Found {value_names}.') + + swath_ids = {} + for swath in set([product.swath for product in products]): + swath_products = [product for product in products if product.swath == swath] + swath_products.sort(key=lambda x: x.burst_id) + ids = np.array([p.burst_id for p in swath_products]) + if not np.all(ids - ids.min() == np.arange(len(ids))): + raise ValueError(f'Products for swath {swath} are not contiguous') + swath_ids[swath] = ids + + for swath1, swath2 in combinations(swath_ids.keys(), 2): + separations = np.concatenate([swath_ids[swath1] - elem for elem in swath_ids[swath2]]) + if separations.min() > 1: + raise ValueError(f'Products from swaths {swath1} and {swath2} do not overlap') + + +def get_product_multilook(product_dir: Path) -> Tuple: + """Get the multilook values for a set of ASF burst products. + You should have already checked that all products have the same multilook, + so you can just use the first product's values. + + Args: + product_dir: The path to the directory containing the UNZIPPED ASF burst product directories + + Returns: + The number of azimuth looks and range looks + """ + product_path = list(product_dir.glob('S1_??????_IW?_*'))[0] + metadata_path = product_path / f'{product_path.name}.txt' + meta = read_product_metadata(metadata_path) + return int(meta['Rangelooks']), int(meta['Azimuthlooks']) + + +def prepare_products(directory: Path) -> None: + """Set up a directory for ISCE2-based burst merging using a set of ASF burst products. + This includes: + - Downloading annotation xml files + - Downloading manifest xml files + - Downloading a DEM + - Spoofing the geom_reference and fine_interferogram directories + + Args: + directory: The path to the directory containing the UNZIPPED ASF burst product directories + """ + product_paths = list(directory.glob('S1_??????_IW?_*')) + products = get_burst_metadata(product_paths) + check_burst_group_validity(products) + download_metadata_xmls([product.to_burst_params() for product in products]) + swaths = list(set([int(product.swath[2:3]) for product in products])) + swath_objs = [] + for swath in swaths: + swath_products = [product for product in products if int(product.swath[2:3]) == swath] + swath_products = sorted(swath_products, key=lambda x: x.start_utc) + swath_obj = create_burst_cropped_s1_obj(swath, swath_products) + for product, burst_obj in zip(swath_products, swath_obj.product.bursts): + product.isce2_burst_number = burst_obj.burstNumber + + multilooked_swath_obj = modify_for_multilook(swath_products, swath_obj) + multilooked_swath_obj.write_xml() + + spoof_isce2_setup(swath_products, swath_obj) + swath_objs.append(copy.deepcopy(swath_obj)) + del swath_obj + + download_dem_for_multiple_bursts(swath_objs) + + +def run_isce2_workflow( + range_looks: int, azimuth_looks: int, mergedir='merged', filter_strength=0.5, apply_water_mask=False +) -> None: + """Run the ISCE2 workflow for burst merging, filtering, unwrapping, and geocoding + + Args: + azimuth_looks: The number of azimuth looks + range_looks: The number of range looks + mergedir: The output directory containing the merged interferogram + filter_strength: The Goldstein-Werner filter strength + apply_water_mask: Whether or not to apply a water body mask to the coherence file before unwrapping + """ + mergedir = Path(mergedir) + mergedir.mkdir(exist_ok=True) + merge_bursts(range_looks, azimuth_looks, merge_dir=str(mergedir)) + goldstein_werner_filter( + mergedir / WRP_IFG_NAME, mergedir / FILT_WRP_IFG_NAME, mergedir / COH_NAME, filter_strength=filter_strength + ) + if apply_water_mask: + log.info('Water masking requested, downloading water mask') + mask_coherence(f'masked.{COH_NAME}') + corrfile = os.path.join(str(mergedir), f'masked.{COH_NAME}') + else: + corrfile = os.path.join(str(mergedir), COH_NAME) + snaphu_unwrap(range_looks, azimuth_looks, corrfile=corrfile, base_dir=str(mergedir)) + geocode_products(range_looks, azimuth_looks, dem_path='full_res.dem.wgs84', base_dir=str(mergedir)) + + +def package_output( + product_directory: Path, looks: str, filter_strength: float, water_mask: bool, archive=False +) -> None: + """Package the output of the ISCE2 workflow into a the standard ASF burst product format + + Args: + product_directory: The path to the directory containing the UNZIPPED ASF burst product directories + looks: The number of looks [20x4, 10x2, 5x1] + filter_strength: The Goldstein-Werner filter strength + archive: Whether or not to create a zip archive of the output + """ + pixel_size = get_pixel_size(looks) + range_looks, azimuth_looks = [int(look) for look in looks.split('x')] + + product_path = list(product_directory.glob('S1_??????_IW?_*'))[0] + example_metadata = get_burst_metadata([product_path])[0] + product_name = get_product_name(example_metadata, pixel_size) + out_product_dir = Path(product_name) + out_product_dir.mkdir(parents=True, exist_ok=True) + + metas = get_product_metadata_info(product_directory) + make_parameter_file( + Path(f'{product_name}/{product_name}.txt'), + metas, + range_looks, + azimuth_looks, + filter_strength, + water_mask, + ) + + translate_outputs(product_name, pixel_size=pixel_size, include_radar=False) + reference_scenes = [meta['ReferenceGranule'] for meta in metas] + secondary_scenes = [meta['SecondaryGranule'] for meta in metas] + make_readme(out_product_dir, reference_scenes, secondary_scenes, range_looks, azimuth_looks, water_mask) + unwrapped_phase = f'{product_name}/{product_name}_unw_phase.tif' + make_browse_image(unwrapped_phase, f'{product_name}/{product_name}_unw_phase.png') + if archive: + make_archive(base_name=product_name, format='zip', base_dir=product_name) + + +def merge_tops_bursts(product_directory: Path, filter_strength: float, apply_water_mask: bool): + """Run the full ISCE2 workflow for TOPS burst merging + + Args: + product_directory: The path to the directory containing the UNZIPPED ASF burst product directories + filter_strength: The Goldstein-Werner filter strength + apply_water_mask: Whether or not to apply a water body mask to the coherence file before unwrapping + """ + range_looks, azimuth_looks = get_product_multilook(product_directory) + prepare_products(product_directory) + run_isce2_workflow(range_looks, azimuth_looks, filter_strength=filter_strength, apply_water_mask=apply_water_mask) + package_output(product_directory, f'{range_looks}x{azimuth_looks}', filter_strength, water_mask=apply_water_mask) + + +def main(): + """HyP3 entrypoint for the TOPS burst merging workflow""" + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('directory', type=str, help='Directory where your unzipped burst InSAR products are located') + parser.add_argument( + '--filter-strength', type=float, default=0.6, help='Goldstein-Werner filter strength (between 0 and 1)' + ) + parser.add_argument( + '--apply-water-mask', + type=string_is_true, + default=False, + help='Apply a water body mask to wrapped and unwrapped phase GeoTIFFs (after unwrapping)', + ) + args = parser.parse_args() + product_directory = Path(args.directory) + merge_tops_bursts(product_directory, args.filter_strength, args.apply_water_mask) + + +if __name__ == '__main__': + main() diff --git a/src/hyp3_isce2/metadata/templates/insar_burst/readme.md.txt.j2 b/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_base.md.txt.j2 similarity index 71% rename from src/hyp3_isce2/metadata/templates/insar_burst/readme.md.txt.j2 rename to src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_base.md.txt.j2 index 536c4a1d..a498c2f3 100644 --- a/src/hyp3_isce2/metadata/templates/insar_burst/readme.md.txt.j2 +++ b/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_base.md.txt.j2 @@ -1,58 +1,5 @@ -ASF Sentinel-1 Burst InSAR Data Package (ISCE2) -=============================================== +{% block header %}{% endblock %} -This folder contains burst-based SAR Interferometry (InSAR) products and their associated files. The source data for -these products are Sentinel-1 bursts, extracted from Single Look Complex (SLC) products processed by ESA, -and they were processed using InSAR Scientific Computing Environment version 2 (ISCE2) software. - -Refer to -https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath -for more information on Sentinel-1 bursts. - -This data was processed by ASF DAAC HyP3 {{ processing_date.year }} using the {{ plugin_name }} plugin version -{{ plugin_version }} running {{ processor_name }} release {{ processor_version }}. -Files are projected to {{ projection }}, and the pixel spacing is {{ pixel_spacing|int }} m. - -The source bursts for this InSAR product are: - - Reference: {{ reference_burst_name }} - - Secondary: {{ secondary_burst_name }} - -Processing Date/Time: {{ processing_date.isoformat(timespec='seconds') }} - -The directory name for this product is: {{ product_name }} - -The output directory uses the following naming convention: - -S1_bbbbbb_xxn_yyyymmdd_yyyymmdd_pp_INTzz_cccc - -bbbbbb: Relative burst ID values assigned by ESA. Each value identifies a consistent burst footprint; relative burst - ID values differ from one subswath to the next. - -xxn: Three character combination that represents the acquisition mode and subswath number. The first two characters - xx represent the acquisition mode (IW or EW) and the last character identifies the subswath number - (1-3 for IW, 1-5 for EW). IW mode indicates Interferometric Wideswath, which acquires a 250 km swath composed - of three subswaths. EW mode indicates Extra-Wide Swath, which acquires a 400 km swath composed of 5 subswaths. - -yyyymmdd: Date of acquisition of the reference and secondary images, respectively. - -pp: Two character combination that represents the mode of radar orientation (polarization) for both signal - transmission and reception. The first position represents the transmit orientation mode and the second - position represents the receive orientation mode. - - HH: Horizontal Transmit - Horizontal Receive - HV: Horizontal Transmit - Vertical Receive - VH: Vertical Transmit - Horizontal Receive - VV: Vertical Transmit - Vertical Receive - -INT: The product type (always INT for InSAR). - -zz: The pixel spacing of the output image. - -cccc: 4-character unique product identifier. - -Files contained in the product directory are named using the directory name followed by a tag indicating the file type. - ----------------- ### Pixel Spacing ### When ordering Sentinel-1 Burst InSAR On Demand products, users can choose the number of looks to use in processing, @@ -144,6 +91,12 @@ The files generated in this process include: 8. Water Mask (GeoTIFF) 9. README.md.txt (Text File) +There are also four non-geocoded GeoTIFFs that remain in their native range-doppler coordinates. These four images compose +the image data needed to merge burst InSAR products together. These images include a range-doppler version of the wrapped +interferogram, a two-band range-doppler look vector image in the native ISCE2 format, and latitude/longitude images that +provide the information necessary to map range-doppler images into the geocoded domain. These images files are not +included in merged burst InSAR products. + *See below for detailed descriptions of each of the product files.* ---------------- @@ -277,102 +230,15 @@ package. Users can choose to apply the water mask before phase unwrapping. This mitigates potential errors in the unwrapping process caused by invalid coherence over water bodies. The water mask will also return nodata values in areas of water -in the output products. This product {{ "has" if apply_water_mask else "has not" }} +in the output products. This product {{ "has" if apply_water_mask else "has not" }} had the water mask applied. -The water mask is generated using data from OpenStreetMaps and/or ESA WorldCover depeding on location. Areas within +The water mask is generated using data from OpenStreetMaps and/or ESA WorldCover depending on location. Areas within Canada, Alaska, and Russia are primarily covered by ESA WorldCover data, while the rest of the world is covered by OpenStreetMaps data. Water masks were previously generated from the GSHHG dataset. ************* -# Burst InSAR Processing # - -The basic steps in Sentinel-1 Burst InSAR processing are as follows: - -*Pre-Processing* -1. Download the Sentinel SLC burst data -2. Calculate the intersection area -3. Download the dem covering the intersection area -4. Download the aux calibration file -5. Download the orbit parameter files -6. Set the Configure file for the InSAR process - -*InSAR Processing* -7. Run topsApp steps 'start' and 'preprocess' -8. Switch the reference and secondary bursts as necessary -9. Run topsApp step 'computeBaselines' to 'filter' -10. Optionally apply the water mask to the wrapped image. -11. Run topsApp steps 'unwrap' and 'unwrap2stage' -12. Run step 'geocode' - -*Post-Processing* -13. translate output files to hyp3-gamma format -14. write the README text file -15. write the metadata txt file -16. zip the output files -17. upload the files to s3 - ----------------- -The detailed process, including the calls to ISCE2 software, is as follows: - -The prepare-processing and InSAR processing are combined in the insar_tps_burst function. - -## Pre-processing ## - - download_bursts: Download the bursts according to reference and secondary parameters and re-create a SAFE-style directory - - get_region_of_interest: Calculate the intersection area - - download_dem_for_isce2 : Download the DEM file - - download_aux_cal: Download: the aux cal file - - downloadSentinelOrbitFile: Download the orbit file - - topsapp.TopsappBurstConfig: Set the configuration file - -## InSAR processing ## -The ISCE2 InSAR processing this product uses includes the following ISCE2 topsApp steps: -- startup -- preprocess -- computeBaselines -- verifyDEM -- topo -- subsetoverlaps -- coarseoffsets -- coarseresamp -- overlapifg -- prepesd -- esd -- rangecoreg -- fineoffsets -- fineresamp -- ion -- burstifg -- mergebursts -- filter -- unwrap -- unwrap2stage -- geocode - -These steps are run using these calls within hyp3-isce2: -- topsapp.run_topsapp_burst(start='startup', end='preprocess', config_xml=config_path): Extract the orbits, - IPF (Instrument Processing Facility) version, burst data, and antenna pattern if it is necessary -- topsapp.swap_burst_vrts(): Switch the reference and secondary bursts to use the burst data download from ASF -- topsapp.run_topsapp_burst(start='computeBaselines', end='unwrap2stage', config_xml=config_path): - Run the remaining processing steps including: - - Calculate the perpendicular and parallel baselines - - Verify the DEM file to make sure it covers the bursts - - Map DEM into the radar coordinates of the reference image. This generates the longitude, - latitude, height and LOS angles on a pixel by pixel grid for each burst. - - Estimate the azimuth offsets between the input SLCs (The Enhanced Spectral Diversity (ESD) method is NOT used) - - Estimate the range offsets between the input SLCS - - Coregister the secondary SLC by applying the estimated range and azimuth offsets - - Produce the wrapped phase interferogram - - If --apply-water-mask is selected, create and apply the water mask to the wrapped interferogram. - - Unwrap the wrapped phase interferogram using SNAPHU to produce the wrapped phase interferogram -- topsapp.run_topsapp_burst(start='geocode', end='geocode', config_xml=config_path): Geocode the output products - -## Post-Processing ## - - translate_outputs: Convert the outputs of hyp3-isce2 to hyp3-gamma formatted geotiff files - - make_readme_file: Produce the readme.md.txt file in the product - - make_parameter_file: Produce metadata text file in the product - - make_archive: Zip the output files - - upload_file_to_s3: Upload the files in the product directory to the s3 +{% block burst_insar_processing %}{% endblock %} ************* # The Sentinel-1 mission # diff --git a/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_merge_readme.md.txt.j2 b/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_merge_readme.md.txt.j2 new file mode 100644 index 00000000..9cbc475b --- /dev/null +++ b/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_merge_readme.md.txt.j2 @@ -0,0 +1,125 @@ +{% extends "insar_burst/insar_burst_base.md.txt.j2" %} + +Note to reader: This readme file includes text blocks to extend the insar burst base file. +Only text included in blocks called in the base file will be included in the output readme. +We have simple placeholders for readability in this file to indicate where the base file will have its own sections. + +{% block header %} +ASF Sentinel-1 Burst InSAR Data Package (ISCE2) +=============================================== + +This folder contains merged burst-based SAR Interferometry (InSAR) products and their associated files. The source data for +these products are Sentinel-1 bursts, extracted from Single Look Complex (SLC) products processed by ESA, +and they were processed using InSAR Scientific Computing Environment version 2 (ISCE2) software. + +Refer to +https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath +for more information on Sentinel-1 bursts. + +This data was processed by ASF DAAC HyP3 {{ processing_date.year }} using the {{ plugin_name }} plugin version +{{ plugin_version }} running {{ processor_name }} release {{ processor_version }}. +Files are projected to {{ projection }}, and the pixel spacing is {{ pixel_spacing|int }} m. + +The source bursts for this InSAR product are: + - Reference: {{ reference_burst_name }} + - Secondary: {{ secondary_burst_name }} + +Processing Date/Time: {{ processing_date.isoformat(timespec='seconds') }} + +The directory name for this product is: {{ product_name }} + +The output directory uses the following naming convention: + +S1_rrr__yyyymmdd_yyyymmdd_pp_INTzz_cccc + +rrr: Relative orbit ID values assigned by ESA. Merged burst InSAR products can contain many relative burst IDs, so the + relate orbit ID is used in lieu of relative burst IDs for these products + +yyyymmdd: Date of acquisition of the reference and secondary images, respectively. + +pp: Two character combination that represents the mode of radar orientation (polarization) for both signal + transmission and reception. The first position represents the transmit orientation mode and the second + position represents the receive orientation mode. + + HH: Horizontal Transmit - Horizontal Receive + HV: Horizontal Transmit - Vertical Receive + VH: Vertical Transmit - Horizontal Receive + VV: Vertical Transmit - Vertical Receive + +INT: The product type (always INT for InSAR). + +zz: The pixel spacing of the output image. + +cccc: 4-character unique product identifier. + +Files contained in the product directory are named using the directory name followed by a tag indicating the file type. +{% endblock %} +---------------- +(This is where the base file has the Pixel Spacing section) + +---------------- +(This is where the base file has the Using This Data section) + +*************** +(This is where the base file has parts 1-8 of the Product Contents) + +************* +{% block burst_insar_processing %} +# Burst InSAR Processing # + +The basic steps in Sentinel-1 Burst InSAR processing are as follows: + +*Pre-Processing* +1. Check that the input burst InSAR products are capable of being merged +2. Recreate the ISCE2 post-interferogram generation directory structure +3. Reformat range-doppler burst InSAR product datasets to an ISCE2-compatible format +4. Create ISCE2 Sentinel-1 objects with the correct burst/multilook information + +*InSAR Processing* +5. Run topsApp step 'mergebursts' +6. Optionally apply the water mask to the wrapped image. +7. Run topsApp steps 'unwrap' and 'unwrap2stage' +8. Run step 'geocode' + +*Post-Processing* +9. translate output files to hyp3 format +10. write the README text file +11. write the metadata txt file + +---------------- +The detailed process, including the calls to ISCE2 software, is as follows: + +The prepare-processing and InSAR processing are combined in the insar_tps_burst function. + +## Pre-processing ## + - merge_tops_bursts.check_burst_group_validity:Check that the input burst InSAR products are capable of being merged + - merge_tops_bursts.download_metadata_xmls: Download metadata files for a set of burst InSAR products + - merge_tops_bursts.create_burst_cropped_s1_obj: Create ISCE2 `Sentinel1` objects for the swaths/bursts present + - merge_tops_bursts.spoof_isce2_setup: Recreate an ISCE2 setup post-`burstifg` + - merge_tops_bursts.download_dem_for_multiple_bursts: Download DEM for merge run + +## InSAR processing ## +The ISCE2 InSAR processing this product uses includes the following ISCE2 topsApp steps: +- mergebursts +- filter +- unwrap +- unwrap2stage +- geocode + +These steps are run using these calls within hyp3-isce2: +- merge_tops_bursts.merge_bursts: Merge the wrapped burst interferograms +- merge_tops_bursts.goldstein_werner_filter: Apply the Goldstien-Werner Phase Filter +- merge_tops_bursts.mask_coherence: Optionally mask data before unwrapping +- merge_tops_bursts.snaphu_unwrap: Unwrap the merged interferogram using SNAPHU +- merge_tops_bursts.geocode_products: Geocode the output products + +## Post-Processing ## + - merge_tops_bursts.make_parameter_file: Produce metadata text file in the product + - merge_tops_bursts.translate_outputs: Convert the outputs of hyp3-isce2 to hyp3-gamma formatted geotiff files + - merge_tops_bursts.make_browse_image: Create a browse image for the dataset + - merge_tops_bursts.make_readme: Produce the readme.md.txt file in the product + {% endblock %} + + ----------- + (This is where the base file has a S1 Mission section) + (This is where the base file has a footer) diff --git a/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_readme.md.txt.j2 b/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_readme.md.txt.j2 new file mode 100644 index 00000000..8b963f27 --- /dev/null +++ b/src/hyp3_isce2/metadata/templates/insar_burst/insar_burst_readme.md.txt.j2 @@ -0,0 +1,165 @@ +{% extends "insar_burst/insar_burst_base.md.txt.j2" %} + +Note to reader: This readme file includes text blocks to extend the insar burst base file. +Only text included in blocks called in the base file will be included in the output readme. +We have simple placeholders for readability in this file to indicate where the base file will have its own sections. + +{% block header %} +ASF Sentinel-1 Burst InSAR Data Package (ISCE2) +=============================================== + +This folder contains burst-based SAR Interferometry (InSAR) products and their associated files. The source data for +these products are Sentinel-1 bursts, extracted from Single Look Complex (SLC) products processed by ESA, +and they were processed using InSAR Scientific Computing Environment version 2 (ISCE2) software. + +Refer to +https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath +for more information on Sentinel-1 bursts. + +This data was processed by ASF DAAC HyP3 {{ processing_date.year }} using the {{ plugin_name }} plugin version +{{ plugin_version }} running {{ processor_name }} release {{ processor_version }}. +Files are projected to {{ projection }}, and the pixel spacing is {{ pixel_spacing|int }} m. + +The source bursts for this InSAR product are: + - Reference: {{ reference_burst_name }} + - Secondary: {{ secondary_burst_name }} + +Processing Date/Time: {{ processing_date.isoformat(timespec='seconds') }} + +The directory name for this product is: {{ product_name }} + +The output directory uses the following naming convention: + +S1_bbbbbb_xxn_yyyymmdd_yyyymmdd_pp_INTzz_cccc + +bbbbbb: Relative burst ID values assigned by ESA. Each value identifies a consistent burst footprint; relative burst + ID values differ from one subswath to the next. + +xxn: Three character combination that represents the acquisition mode and subswath number. The first two characters + xx represent the acquisition mode (IW or EW) and the last character identifies the subswath number + (1-3 for IW, 1-5 for EW). IW mode indicates Interferometric Wideswath, which acquires a 250 km swath composed + of three subswaths. EW mode indicates Extra-Wide Swath, which acquires a 400 km swath composed of 5 subswaths. + +yyyymmdd: Date of acquisition of the reference and secondary images, respectively. + +pp: Two character combination that represents the mode of radar orientation (polarization) for both signal + transmission and reception. The first position represents the transmit orientation mode and the second + position represents the receive orientation mode. + + HH: Horizontal Transmit - Horizontal Receive + HV: Horizontal Transmit - Vertical Receive + VH: Vertical Transmit - Horizontal Receive + VV: Vertical Transmit - Vertical Receive + +INT: The product type (always INT for InSAR). + +zz: The pixel spacing of the output image. + +cccc: 4-character unique product identifier. + +Files contained in the product directory are named using the directory name followed by a tag indicating the file type. +{% endblock %} +---------------- +(This is where the base file has the Pixel Spacing section) + +---------------- +(This is where the base file has the Using This Data section) + +*************** +(This is where the base file has parts 1-8 of the Product Contents) + +************* +{% block burst_insar_processing %} +# Burst InSAR Processing # + +The basic steps in Sentinel-1 Burst InSAR processing are as follows: + +*Pre-Processing* +1. Download the Sentinel SLC burst data +2. Calculate the intersection area +3. Download the dem covering the intersection area +4. Download the aux calibration file +5. Download the orbit parameter files +6. Set the Configure file for the InSAR process + +*InSAR Processing* +7. Run topsApp steps 'start' and 'preprocess' +8. Switch the reference and secondary bursts as necessary +9. Run topsApp step 'computeBaselines' to 'filter' +10. Optionally apply the water mask to the wrapped image. +11. Run topsApp steps 'unwrap' and 'unwrap2stage' +12. Run step 'geocode' + +*Post-Processing* +13. translate output files to hyp3-gamma format +14. write the README text file +15. write the metadata txt file +16. zip the output files +17. upload the files to s3 + +---------------- +The detailed process, including the calls to ISCE2 software, is as follows: + +The prepare-processing and InSAR processing are combined in the insar_tps_burst function. + +## Pre-processing ## + - download_bursts: Download the bursts according to reference and secondary parameters and re-create a SAFE-style directory + - get_region_of_interest: Calculate the intersection area + - download_dem_for_isce2 : Download the DEM file + - download_aux_cal: Download: the aux cal file + - downloadSentinelOrbitFile: Download the orbit file + - topsapp.TopsappBurstConfig: Set the configuration file + +## InSAR processing ## +The ISCE2 InSAR processing this product uses includes the following ISCE2 topsApp steps: +- startup +- preprocess +- computeBaselines +- verifyDEM +- topo +- subsetoverlaps +- coarseoffsets +- coarseresamp +- overlapifg +- prepesd +- esd +- rangecoreg +- fineoffsets +- fineresamp +- ion +- burstifg +- mergebursts +- filter +- unwrap +- unwrap2stage +- geocode + +These steps are run using these calls within hyp3-isce2: +- topsapp.run_topsapp_burst(start='startup', end='preprocess', config_xml=config_path): Extract the orbits, + IPF (Instrument Processing Facility) version, burst data, and antenna pattern if it is necessary +- topsapp.swap_burst_vrts(): Switch the reference and secondary bursts to use the burst data download from ASF +- topsapp.run_topsapp_burst(start='computeBaselines', end='unwrap2stage', config_xml=config_path): + Run the remaining processing steps including: + - Calculate the perpendicular and parallel baselines + - Verify the DEM file to make sure it covers the bursts + - Map DEM into the radar coordinates of the reference image. This generates the longitude, + latitude, height and LOS angles on a pixel by pixel grid for each burst. + - Estimate the azimuth offsets between the input SLCs (The Enhanced Spectral Diversity (ESD) method is NOT used) + - Estimate the range offsets between the input SLCS + - Coregister the secondary SLC by applying the estimated range and azimuth offsets + - Produce the wrapped phase interferogram + - If --apply-water-mask is selected, create and apply the water mask to the wrapped interferogram. + - Unwrap the wrapped phase interferogram using SNAPHU to produce the wrapped phase interferogram +- topsapp.run_topsapp_burst(start='geocode', end='geocode', config_xml=config_path): Geocode the output products + +## Post-Processing ## + - translate_outputs: Convert the outputs of hyp3-isce2 to hyp3-gamma formatted geotiff files + - make_readme_file: Produce the readme.md.txt file in the product + - make_parameter_file: Produce metadata text file in the product + - make_archive: Zip the output files + - upload_file_to_s3: Upload the files in the product directory to the s3 bucket + {% endblock %} + + ----------- + (This is where the base file has a S1 Mission section) + (This is where the base file has a footer) \ No newline at end of file diff --git a/src/hyp3_isce2/utils.py b/src/hyp3_isce2/utils.py index 97d1a624..dbbd9899 100644 --- a/src/hyp3_isce2/utils.py +++ b/src/hyp3_isce2/utils.py @@ -2,15 +2,19 @@ import os import shutil import subprocess +from dataclasses import dataclass +from datetime import datetime from pathlib import Path from platform import system -from typing import Tuple +from typing import Optional, Tuple import isceobj import numpy as np from isceobj.Util.ImageUtil.ImageLib import loadImage +from iscesys.Component.ProductManager import ProductManager from osgeo import gdal + gdal.UseExceptions() ESA_HOST = 'dataspace.copernicus.eu' @@ -39,13 +43,103 @@ def __exit__(self, exc_type, exc_val, exc_tb): gdal.SetConfigOption(key, value) +@dataclass +class ParameterFile: + reference_granule: str + secondary_granule: str + reference_orbit_direction: str + reference_orbit_number: str + secondary_orbit_direction: str + secondary_orbit_number: str + baseline: float + utc_time: float + heading: float + spacecraft_height: float + earth_radius_at_nadir: float + slant_range_near: float + slant_range_center: float + slant_range_far: float + range_looks: int + azimuth_looks: int + insar_phase_filter: bool + phase_filter_parameter: float + range_bandpass_filter: bool + azimuth_bandpass_filter: bool + dem_source: str + dem_resolution: int + unwrapping_type: str + speckle_filter: bool + water_mask: bool + radar_n_lines: Optional[int] = None + radar_n_samples: Optional[int] = None + radar_first_valid_line: Optional[int] = None + radar_n_valid_lines: Optional[int] = None + radar_first_valid_sample: Optional[int] = None + radar_n_valid_samples: Optional[int] = None + multilook_azimuth_time_interval: Optional[float] = None + multilook_range_pixel_size: Optional[float] = None + radar_sensing_stop: Optional[datetime] = None + + def __str__(self): + output_strings = [ + f'Reference Granule: {self.reference_granule}\n', + f'Secondary Granule: {self.secondary_granule}\n', + f'Reference Pass Direction: {self.reference_orbit_direction}\n', + f'Reference Orbit Number: {self.reference_orbit_number}\n', + f'Secondary Pass Direction: {self.secondary_orbit_direction}\n', + f'Secondary Orbit Number: {self.secondary_orbit_number}\n', + f'Baseline: {self.baseline}\n', + f'UTC time: {self.utc_time}\n', + f'Heading: {self.heading}\n', + f'Spacecraft height: {self.spacecraft_height}\n', + f'Earth radius at nadir: {self.earth_radius_at_nadir}\n', + f'Slant range near: {self.slant_range_near}\n', + f'Slant range center: {self.slant_range_center}\n', + f'Slant range far: {self.slant_range_far}\n', + f'Range looks: {self.range_looks}\n', + f'Azimuth looks: {self.azimuth_looks}\n', + f'INSAR phase filter: {"yes" if self.insar_phase_filter else "no"}\n', + f'Phase filter parameter: {self.phase_filter_parameter}\n', + f'Range bandpass filter: {"yes" if self.range_bandpass_filter else "no"}\n', + f'Azimuth bandpass filter: {"yes" if self.azimuth_bandpass_filter else "no"}\n', + f'DEM source: {self.dem_source}\n', + f'DEM resolution (m): {self.dem_resolution}\n', + f'Unwrapping type: {self.unwrapping_type}\n', + f'Speckle filter: {"yes" if self.speckle_filter else "no"}\n', + f'Water mask: {"yes" if self.water_mask else "no"}\n', + ] + + # TODO could use a more robust way to check if radar data is present + if self.radar_n_lines: + radar_data = [ + f'Radar n lines: {self.radar_n_lines}\n', + f'Radar n samples: {self.radar_n_samples}\n', + f'Radar first valid line: {self.radar_first_valid_line}\n', + f'Radar n valid lines: {self.radar_n_valid_lines}\n', + f'Radar first valid sample: {self.radar_first_valid_sample}\n', + f'Radar n valid samples: {self.radar_n_valid_samples}\n', + f'Multilook azimuth time interval: {self.multilook_azimuth_time_interval}\n', + f'Multilook range pixel size: {self.multilook_range_pixel_size}\n', + f'Radar sensing stop: {datetime.strftime(self.radar_sensing_stop, "%Y-%m-%dT%H:%M:%S.%f")}\n', + ] + output_strings += radar_data + + return ''.join(output_strings) + + def __repr__(self): + return self.__str__() + + def write(self, out_path: Path): + out_path.write_text(self.__str__()) + + def get_esa_credentials() -> Tuple[str, str]: netrc_name = '_netrc' if system().lower() == 'windows' else '.netrc' netrc_file = Path.home() / netrc_name - if "ESA_USERNAME" in os.environ and "ESA_PASSWORD" in os.environ: - username = os.environ["ESA_USERNAME"] - password = os.environ["ESA_PASSWORD"] + if 'ESA_USERNAME' in os.environ and 'ESA_PASSWORD' in os.environ: + username = os.environ['ESA_USERNAME'] + password = os.environ['ESA_PASSWORD'] return username, password if netrc_file.exists(): @@ -56,8 +150,8 @@ def get_esa_credentials() -> Tuple[str, str]: return username, password raise ValueError( - "Please provide Copernicus Data Space Ecosystem (CDSE) credentials via the " - "ESA_USERNAME and ESA_PASSWORD environment variables, or your netrc file." + 'Please provide Copernicus Data Space Ecosystem (CDSE) credentials via the ' + 'ESA_USERNAME and ESA_PASSWORD environment variables, or your netrc file.' ) @@ -101,14 +195,15 @@ def extent_from_geotransform(geotransform: tuple, x_size: int, y_size: int) -> t def make_browse_image(input_tif: str, output_png: str) -> None: with GDALConfigManager(GDAL_PAM_ENABLED='NO'): stats = gdal.Info(input_tif, format='json', stats=True)['stac']['raster:bands'][0]['stats'] - gdal.Translate(destName=output_png, - srcDS=input_tif, - format='png', - outputType=gdal.GDT_Byte, - width=2048, - strict=True, - scaleParams=[[stats['minimum'], stats['maximum']]], - ) + gdal.Translate( + destName=output_png, + srcDS=input_tif, + format='png', + outputType=gdal.GDT_Byte, + width=2048, + strict=True, + scaleParams=[[stats['minimum'], stats['maximum']]], + ) def oldest_granule_first(g1, g2): @@ -118,7 +213,7 @@ def oldest_granule_first(g1, g2): def load_isce2_image(in_path) -> tuple[isceobj.Image, np.ndarray]: - """ Read an ISCE2 image file and return the image object and array. + """Read an ISCE2 image file and return the image object and array. Args: in_path: The path to the image to resample (not the xml). @@ -129,28 +224,48 @@ def load_isce2_image(in_path) -> tuple[isceobj.Image, np.ndarray]: """ image_obj, _, _ = loadImage(in_path) array = np.fromfile(in_path, image_obj.toNumpyDataType()) + array = np.reshape(array, (-1, image_obj.width)) + if image_obj.bands > 1: + if image_obj.imageType == 'bil': + shape = (image_obj.bands, image_obj.length, image_obj.width) + new_array = np.zeros(shape, dtype=image_obj.toNumpyDataType()) + for i in range(image_obj.bands): + new_array[i, :, :] = array[i:: image_obj.bands] + array = new_array.copy() + else: + raise NotImplementedError('Non-BIL reading is not implemented') return image_obj, array -def write_isce2_image(output_path, array=None, width=None, mode='read', data_type='FLOAT') -> None: - """ Write an ISCE2 image file. +def write_isce2_image(output_path: str, array: np.ndarray) -> None: + """Write a numpy array as an ISCE2 image file. Args: output_path: The path to the output image file. array: The array to write to the file. - width: The width of the image. - mode: The mode to open the image in. - data_type: The data type of the image. """ - if array is not None: - array.tofile(output_path) - width = array.shape[1] - elif width is None: - raise ValueError('Either a width or an input array must be provided') - - out_obj = isceobj.createImage() - out_obj.initImage(output_path, mode, width, data_type) - out_obj.renderHdr() + data_type_dic = {'float32': 'FLOAT', 'float64': 'DOUBLE', 'int32': 'INT', 'complex64': 'CFLOAT', 'int8': 'BYTE'} + + data_type = data_type_dic[str(array.dtype)] + + if array.ndim == 1: + bands = 1 + length = 1 + width = array.shape[0] + elif array.ndim == 2: + bands = 1 + length, width = array.shape + elif array.ndim == 3: + bands, length, width = array.shape + else: + raise NotImplementedError('array with dimension larger than 3 is not implemented') + + image_obj = isceobj.createImage() + image_obj.initImage(output_path, 'write', width, data_type, bands) + image_obj.setLength(length) + image_obj.setImageType('bil') + image_obj.createImage() + write_isce2_image_from_obj(image_obj, array) def get_geotransform_from_dataset(dataset: isceobj.Image) -> tuple: @@ -171,12 +286,7 @@ def get_geotransform_from_dataset(dataset: isceobj.Image) -> tuple: def resample_to_radar( - mask: np.ndarray, - lat: np.ndarray, - lon: np.ndarray, - geotransform: tuple, - data_type: type, - outshape: tuple[int, int] + mask: np.ndarray, lat: np.ndarray, lon: np.ndarray, geotransform: tuple, data_type: type, outshape: tuple[int, int] ) -> np.ndarray: """Resample a geographic image to radar coordinates using a nearest neighbor method. The latin and lonin images are used to map from geographic to radar coordinates. @@ -217,15 +327,16 @@ def resample_to_radar_io(image_to_resample: str, latin: str, lonin: str, output: _, lon = load_isce2_image(lonin) mask = np.reshape(mask, [maskim.coord2.coordSize, maskim.coord1.coordSize]) geotransform = get_geotransform_from_dataset(maskim) - cropped = resample_to_radar(mask=mask, - lat=lat, - lon=lon, - geotransform=geotransform, - data_type=maskim.toNumpyDataType(), - outshape=(latim.coord2.coordSize, latim.coord1.coordSize) - ) + cropped = resample_to_radar( + mask=mask, + lat=lat, + lon=lon, + geotransform=geotransform, + data_type=maskim.toNumpyDataType(), + outshape=(latim.coord2.coordSize, latim.coord1.coordSize), + ) - write_isce2_image(output, array=cropped, data_type=maskim.dataType) + write_isce2_image(output, array=cropped) def isce2_copy(in_path: str, out_path: str): @@ -251,5 +362,113 @@ def image_math(image_a_path: str, image_b_path: str, out_path: str, expression: out_path: The path to the output image. expression: The expression to pass to imageMath.py. """ - cmd = ['imageMath.py', '-e', expression, f'--a={image_a_path}', f'--b={image_b_path}', '-o', out_path] + cmd = ['imageMath.py', f'--a={image_a_path}', f'--b={image_b_path}', '-o', f'{out_path}', '--eval', expression] subprocess.run(cmd, check=True) + + +def load_product(xmlname: str): + """Load an ISCE2 product from an xml file + + Args: + xmlname: The path to the xml file + + Returns: + The ISCE2 product + """ + pm = ProductManager() + pm.configure() + obj = pm.loadProduct(xmlname) + return obj + + +def write_isce2_image_from_obj(image_obj, array): + """Write an ISCE2 image file. + + Args: + image_obj: ISCE2 image object + array: The array to write to the file. + """ + image_obj.renderHdr() + + if image_obj.bands > 1: + if image_obj.imageType == 'bil': + shape = (image_obj.length * image_obj.bands, image_obj.width) + new_array = np.zeros(shape, dtype=image_obj.toNumpyDataType()) + for i in range(image_obj.bands): + new_array[i:: image_obj.bands] = array[i, :, :] + array = new_array.copy() + else: + raise NotImplementedError('Non-BIL writing is not implemented') + + array.tofile(image_obj.filename) + + +def create_image( + out_path: str, + width: Optional[int] = None, + access_mode: str = 'read', + image_subtype: str = 'default', + action: str = 'create', +) -> isceobj.Image: + """Create an ISCE2 image object from a set of parameters + + Args: + out_path: The path to the output image + width: The width of the image + access_mode: The access mode of the image (read or write) + image_subtype: The type of image to create + action: What action to take: + 'create': create a new image object, but don't write metadata files, access_mode='write' + 'finalize': create a new image object based on existed binary file, and write metadata files, + access_mode='read' + 'load': create an image object by loading an existing metadata file, access_mode='read' + + Returns: + The ISCE2 image object + """ + opts = { + 'ifg': (isceobj.createIntImage, 1, 'CFLOAT', 'cpx'), + 'cor': (isceobj.createImage, 1, 'FLOAT', 'cor'), + 'unw': (isceobj.Image.createUnwImage, 2, 'FLOAT', 'unw'), + 'conncomp': (isceobj.createImage, 1, 'BYTE', ''), + 'default': (isceobj.createImage, 1, 'FLOAT', ''), + } + + create_func, bands, dtype, image_type = opts[image_subtype] + image = create_func() + if action == 'load': + image.load(out_path + '.xml') + image.setAccessMode('read') + image.createImage() + return image + + if width is None: + raise ValueError('Width must be specified when the action is create or finalize') + + image.initImage(out_path, access_mode, width, dtype, bands) + image.setImageType(image_type) + if action == 'create': + image.createImage() + elif action == 'finalize': + image.renderVRT() + image.createImage() + image.finalizeImage() + image.renderHdr() + return image + + +def read_product_metadata(meta_file_path: str) -> dict: + """Read the HyP3-generated metadata file for a HyP3 product + + Args: + meta_file_path: The path to the metadata file + Returns: + A dictionary of metadata values + """ + hyp3_meta = {} + with open(meta_file_path) as f: + for line in f: + key, *values = line.strip().replace(' ', '').split(':') + value = ':'.join(values) + hyp3_meta[key] = value + return hyp3_meta diff --git a/src/hyp3_isce2/water_mask.py b/src/hyp3_isce2/water_mask.py index 4b94b14f..43092a38 100644 --- a/src/hyp3_isce2/water_mask.py +++ b/src/hyp3_isce2/water_mask.py @@ -1,5 +1,5 @@ """Create and apply a water body mask""" -from os import system +import subprocess from pathlib import Path from typing import Optional @@ -97,13 +97,13 @@ def create_water_mask(input_image: str, output_image: str, gdal_format='ISCE', t # This is WAY faster than using gdal_merge, because of course it is. if len(tiles) > 1: - build_vrt_command = ' '.join(['gdalbuildvrt', merged_vrt_path] + tiles) - system(build_vrt_command) - translate_command = ' '.join(['gdal_translate', merged_vrt_path, merged_tif_path]) - system(translate_command) + build_vrt_command = ['gdalbuildvrt', merged_vrt_path] + tiles + subprocess.run(build_vrt_command, check=True) + translate_command = ['gdal_translate', merged_vrt_path, merged_tif_path] + subprocess.run(translate_command, check=True) - shapefile_command = ' '.join(['gdaltindex', shape_path, input_image]) - system(shapefile_command) + shapefile_command = ['gdaltindex', shape_path, input_image] + subprocess.run(shapefile_command, check=True) warp_filename = merged_tif_path if len(tiles) > 1 else tiles[0] @@ -119,12 +119,12 @@ def create_water_mask(input_image: str, output_image: str, gdal_format='ISCE', t format='GTiff' ) - flip_values_command = ' '.join([ + flip_values_command = [ 'gdal_calc.py', '-A', merged_warped_path, f'--outfile={output_image}', '--calc="numpy.abs((A.astype(numpy.int16) + 1) - 2)"', # Change 1's to 0's and 0's to 1's. f'--format={gdal_format}' - ]) - system(flip_values_command) + ] + subprocess.run(flip_values_command, check=True) diff --git a/tests/conftest.py b/tests/conftest.py index 3918913a..02080df5 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,10 +1,136 @@ import os +import shutil +from datetime import datetime from pathlib import Path +import lxml.etree as ET import pytest +import hyp3_isce2.burst as burst_utils +import hyp3_isce2.merge_tops_bursts as merge + @pytest.fixture() def test_data_dir(): here = Path(os.path.dirname(__file__)) return here / 'data' + + +@pytest.fixture() +def test_merge_dir(test_data_dir): + merge_dir = test_data_dir / 'merge' + + if not merge_dir.exists(): + print('Unzipping merge test data...') + merge_zip = merge_dir.with_suffix('.zip') + + if not merge_zip.exists(): + raise ValueError('merge data not present, run data/create_merge_test_data.py') + + shutil.unpack_archive(merge_zip, test_data_dir) + + return merge_dir + + +@pytest.fixture() +def annotation_manifest_dirs(tmp_path, test_data_dir): + annotation_dir, manifest_dir = merge.prep_metadata_dirs(tmp_path) + sample_xml = ET.parse(test_data_dir / 'reference_descending.xml').getroot() + + et_args = {'encoding': 'UTF-8', 'xml_declaration': True} + param = burst_utils.BurstParams( + 'S1A_IW_SLC__1SDV_20200604T022251_20200604T022318_032861_03CE65_7C85', 'IW2', 'VV', 1 + ) + burst_metadata = burst_utils.BurstMetadata(sample_xml, param) + ET.ElementTree(burst_metadata.annotation).write(annotation_dir / burst_metadata.annotation_name, **et_args) + ET.ElementTree(burst_metadata.manifest).write(manifest_dir / f'{burst_metadata.safe_name}.xml', **et_args) + return annotation_dir, manifest_dir + + +@pytest.fixture +def burst_product1(test_merge_dir): + product_path1 = list(test_merge_dir.glob('S1_136231*'))[0] + product1 = merge.BurstProduct( + granule='S1A_IW_SLC__1SDV_20200604T022251_20200604T022318_032861_03CE65_7C85', + reference_date=datetime(2020, 6, 4, 2, 23, 12), + secondary_date=datetime(2020, 6, 16, 2, 23, 13), + burst_id=136231, + swath='IW2', + polarization='VV', + burst_number=7, + product_path=product_path1, + n_lines=377, + n_samples=1272, + range_looks=20, + azimuth_looks=4, + first_valid_line=8, + n_valid_lines=363, + first_valid_sample=9, + n_valid_samples=1220, + az_time_interval=0.008222225199999992, + rg_pixel_size=46.59124229430646, + start_utc=datetime(2020, 6, 4, 2, 23, 13, 963847), + stop_utc=datetime(2020, 6, 4, 2, 23, 16, 30988), + relative_orbit=64, + ) + return product1 + + +@pytest.fixture +def burst_product2(test_merge_dir): + product_path2 = list(test_merge_dir.glob('*'))[0] + + product2 = merge.BurstProduct( + granule='S1A_IW_SLC__1SDV_20200604T022251_20200604T022318_032861_03CE65_7C85', + reference_date=datetime(2020, 6, 4, 2, 23, 15), + secondary_date=datetime(2020, 6, 16, 2, 23, 16), + burst_id=136232, + swath='IW2', + polarization='VV', + burst_number=8, + product_path=product_path2, + n_lines=377, + n_samples=1272, + range_looks=20, + azimuth_looks=4, + first_valid_line=8, + n_valid_lines=363, + first_valid_sample=9, + n_valid_samples=1220, + az_time_interval=0.008222225199999992, + rg_pixel_size=46.59124229430646, + start_utc=datetime(2020, 6, 4, 2, 23, 16, 722124), + stop_utc=datetime(2020, 6, 4, 2, 23, 18, 795712), + relative_orbit=64, + ) + return product2 + + +@pytest.fixture +def burst_products(burst_product1, burst_product2): + return [burst_product1, burst_product2] + + +@pytest.fixture +def test_s1_obj(annotation_manifest_dirs, burst_products): + annotation_dir, manifest_dir = annotation_manifest_dirs + s1_obj = merge.load_isce_s1_obj(2, 'VV', annotation_dir.parent) + s1_obj.output = str(annotation_dir.parent / burst_products[0].swath) + s1_obj.select_bursts([x.start_utc for x in burst_products]) + s1_obj.update_burst_properties(burst_products) + return s1_obj + + +@pytest.fixture +def isce2_merge_setup(annotation_manifest_dirs, burst_products): + base_dir = annotation_manifest_dirs[0].parent + s1_obj = merge.create_burst_cropped_s1_obj(2, burst_products, 'VV', base_dir=base_dir) + for product, burst_obj in zip(burst_products, s1_obj.product.bursts): + product.isce2_burst_number = burst_obj.burstNumber + + save_dir = str(base_dir / 'fine_interferogram') + multilooked_swath_obj = merge.modify_for_multilook(burst_products, s1_obj, save_dir) + multilooked_swath_obj.write_xml() + + merge.spoof_isce2_setup(burst_products, s1_obj, base_dir=base_dir) + return base_dir diff --git a/tests/data/create_merge_test_data.py b/tests/data/create_merge_test_data.py new file mode 100644 index 00000000..ef66342a --- /dev/null +++ b/tests/data/create_merge_test_data.py @@ -0,0 +1,62 @@ +import os +import shutil +import subprocess +from pathlib import Path + +import numpy as np +from osgeo import gdal + + +def create_product(tmp_path: Path, out_path: Path, granule1: str, granule2: str): + cwd = Path.cwd() + + tmp_path.mkdir(exist_ok=True) + os.chdir(tmp_path) + + cmd = f'insar_tops_burst {granule1} {granule2}' + subprocess.run(cmd.split(' '), check=True) + result_directory = [x for x in tmp_path.glob('S1_*') if x.is_dir()][0].absolute() + + not_radar = [x for x in result_directory.glob('*') if 'rdr' not in x.name] + also_not_metadata = [x for x in not_radar if not x.name == f'{x.parent.name}.txt'] + for output_file in also_not_metadata: + output_file.unlink() + + os.chdir(cwd) + shutil.copytree(result_directory, out_path / result_directory.name) + shutil.rmtree(result_directory.parent) + + +def replace_geotiff_data(geotiff_path: str, stock_value: float) -> None: + """Load a geotiff file, swap the non-zero data for all ones and resave.""" + ds = gdal.Open(geotiff_path, gdal.GA_Update) + band = ds.GetRasterBand(1) + data = band.ReadAsArray() + data[data != 0 + 0j] = stock_value + band.WriteArray(data) + ds.FlushCache() + ds = None + + +def create_test_products(): + out_dir = Path.cwd() / 'merge' + + pairs = [ + ['S1_136231_IW2_20200604T022312_VV_7C85-BURST', 'S1_136231_IW2_20200616T022313_VV_5D11-BURST'], + ['S1_136232_IW2_20200604T022315_VV_7C85-BURST', 'S1_136232_IW2_20200616T022316_VV_5D11-BURST'], + ] + for granule1, granule2 in pairs: + create_product(Path.cwd() / 'tmp', out_dir, granule1, granule2) + + ifgs = out_dir.glob('./*/*_wrapped_phase_rdr.tif') + value = (1 + 0j) + for ifg in ifgs: + value += (0 + 1j) * np.pi / 2 + replace_geotiff_data(str(ifg), value) + + shutil.make_archive(base_name=out_dir.name, format='zip', base_dir=out_dir.name) + shutil.rmtree(out_dir) + + +if __name__ == '__main__': + create_test_products() diff --git a/tests/data/merge.zip b/tests/data/merge.zip new file mode 100644 index 00000000..cf01af7f Binary files /dev/null and b/tests/data/merge.zip differ diff --git a/tests/test_burst.py b/tests/test_burst.py index 7c789e1b..645a6759 100644 --- a/tests/test_burst.py +++ b/tests/test_burst.py @@ -1,13 +1,17 @@ +"""Tests for the single-burst specific functionality found in burst.py""" +from collections import namedtuple +from datetime import datetime from pathlib import Path from re import match from unittest.mock import patch import asf_search +import numpy as np import pytest from lxml import etree from shapely import geometry -from hyp3_isce2 import burst +from hyp3_isce2 import burst, utils URL_BASE = 'https://datapool.asf.alaska.edu/SLC' @@ -23,6 +27,20 @@ def load_metadata(metadata): return xml +def make_test_image(output_path, array=None): + parent = Path(output_path).parent + if not parent.exists(): + parent.mkdir(parents=True, exist_ok=True) + + if array is None: + array = np.zeros((100, 100), dtype=np.float32) + array[10:90, 10:90] = 1 + + array = array.astype(np.float32) + img_obj = utils.create_image(output_path, array.shape[1], access_mode='write', action='create') + utils.write_isce2_image_from_obj(img_obj, array) + + @pytest.mark.parametrize( 'pattern', ( @@ -90,31 +108,30 @@ def test_get_region_of_interest(tmp_path, orbit): def test_get_product_name(): - - reference_name = "S1_136231_IW2_20200604T022312_VV_7C85-BURST" - secondary_name = "S1_136231_IW2_20200616T022313_VV_5D11-BURST" + reference_name = 'S1_136231_IW2_20200604T022312_VV_7C85-BURST' + secondary_name = 'S1_136231_IW2_20200616T022313_VV_5D11-BURST' name_20m = burst.get_product_name(reference_name, secondary_name, pixel_spacing=20.0) name_80m = burst.get_product_name(reference_name, secondary_name, pixel_spacing=80) - assert match("[A-F0-9]{4}", name_20m[-4:]) is not None - assert match("[A-F0-9]{4}", name_80m[-4:]) is not None + assert match('[A-F0-9]{4}', name_20m[-4:]) is not None + assert match('[A-F0-9]{4}', name_80m[-4:]) is not None assert name_20m.startswith('S1_136231_IW2_20200604_20200616_VV_INT20') assert name_80m.startswith('S1_136231_IW2_20200604_20200616_VV_INT80') def mock_asf_search_results( - slc_name: str, - subswath: str, - polarization: str, - burst_index: int) -> asf_search.ASFSearchResults: + slc_name: str, subswath: str, polarization: str, burst_index: int +) -> asf_search.ASFSearchResults: product = asf_search.ASFProduct() product.umm = {'InputGranules': [slc_name]} - product.properties.update({ - 'burst': {'subswath': subswath, 'burstIndex': burst_index}, - 'polarization': polarization, - }) + product.properties.update( + { + 'burst': {'subswath': subswath, 'burstIndex': burst_index}, + 'polarization': polarization, + } + ) results = asf_search.ASFSearchResults([product]) results.searchComplete = True return results @@ -171,22 +188,131 @@ def test_get_burst_params_multiple_results(): def test_validate_bursts(): - burst.validate_bursts( - 'S1_030349_IW1_20230808T171601_VV_4A37-BURST', - 'S1_030349_IW1_20230820T171602_VV_5AC3-BURST' - ) + burst.validate_bursts('S1_030349_IW1_20230808T171601_VV_4A37-BURST', 'S1_030349_IW1_20230820T171602_VV_5AC3-BURST') with pytest.raises(ValueError, match=r'.*polarizations are not the same.*'): burst.validate_bursts( - 'S1_215032_IW2_20230802T144608_VV_7EE2-BURST', - 'S1_215032_IW2_20230721T144607_HH_B3FA-BURST' + 'S1_215032_IW2_20230802T144608_VV_7EE2-BURST', 'S1_215032_IW2_20230721T144607_HH_B3FA-BURST' ) with pytest.raises(ValueError, match=r'.*burst IDs are not the same.*'): burst.validate_bursts( - 'S1_030349_IW1_20230808T171601_VV_4A37-BURST', - 'S1_030348_IW1_20230820T171602_VV_5AC3-BURST' + 'S1_030349_IW1_20230808T171601_VV_4A37-BURST', 'S1_030348_IW1_20230820T171602_VV_5AC3-BURST' ) with pytest.raises(ValueError, match=r'.*only VV and HH.*'): burst.validate_bursts( - 'S1_030349_IW1_20230808T171601_VH_4A37-BURST', - 'S1_030349_IW1_20230820T171602_VH_5AC3-BURST' + 'S1_030349_IW1_20230808T171601_VH_4A37-BURST', 'S1_030349_IW1_20230820T171602_VH_5AC3-BURST' ) + + +def test_load_burst_position(tmpdir): + product = namedtuple('product', ['bursts']) + bursts = namedtuple( + 'bursts', + [ + 'numberOfLines', + 'numberOfSamples', + 'firstValidLine', + 'numValidLines', + 'firstValidSample', + 'numValidSamples', + 'azimuthTimeInterval', + 'rangePixelSize', + 'sensingStop', + ], + ) + + mock_product = product([bursts(100, 200, 10, 50, 20, 60, 0.1, 0.2, datetime(2020, 1, 1))]) + with patch('hyp3_isce2.burst.load_product') as mock_load_product: + mock_load_product.return_value = mock_product + position = burst.load_burst_position('', 0) + + assert position.n_lines == 100 + assert position.n_samples == 200 + assert position.first_valid_line == 10 + assert position.n_valid_lines == 50 + assert position.first_valid_sample == 20 + assert position.n_valid_samples == 60 + assert position.azimuth_time_interval == 0.1 + assert position.range_pixel_size == 0.2 + assert position.sensing_stop == datetime(2020, 1, 1) + + +def test_evenize(): + length, valid_start, valid_length = burst.evenize(101, 7, 57, 5) + assert length == 100 + assert valid_start == 10 + assert valid_length == 50 + + length, valid_start, valid_length = burst.evenize(100, 5, 55, 5) + assert length == 100 + assert valid_start == 5 + assert valid_length == 55 + + with pytest.raises(ValueError, match=r'.*valid data region.*'): + burst.evenize(20, 6, 20, 5) + + +def test_evenly_subset_position(): + input_pos = burst.BurstPosition(101, 101, 11, 20, 11, 20, 1, 1, datetime(2021, 1, 1, 0, 0, 0)) + ml_params = burst.evenly_subset_position(input_pos, 2, 10) + + assert ml_params.n_lines == 100 + assert ml_params.n_samples == 100 + assert ml_params.n_valid_lines == 10 + assert ml_params.n_valid_samples == 18 + assert ml_params.first_valid_line == 20 + assert ml_params.first_valid_sample == 12 + assert ml_params.azimuth_time_interval == input_pos.azimuth_time_interval + assert ml_params.range_pixel_size == input_pos.range_pixel_size + assert ml_params.sensing_stop == datetime(2020, 12, 31, 23, 59, 59) + + +def test_multilook_position(): + input_pos = burst.BurstPosition(100, 100, 20, 60, 20, 30, 1, 1, datetime(2021, 1, 1, 0, 0, 0)) + output_pos = burst.multilook_position(input_pos, 10, 2) + + assert output_pos.n_lines == 50 + assert output_pos.n_samples == 10 + assert output_pos.first_valid_line == 10 + assert output_pos.n_valid_lines == 30 + assert output_pos.first_valid_sample == 2 + assert output_pos.n_valid_samples == 3 + assert output_pos.azimuth_time_interval == input_pos.azimuth_time_interval * 2 + assert output_pos.range_pixel_size == input_pos.range_pixel_size * 10 + assert output_pos.sensing_stop == input_pos.sensing_stop + + +def test_safely_multilook(tmpdir): + image_path = str(tmpdir / 'image') + make_test_image(image_path) + pos = burst.BurstPosition(100, 100, 20, 60, 20, 60, 0.1, 0.1, datetime(2021, 1, 1, 0, 0, 0)) + burst.safely_multilook(image_path, pos, 5, 5) + _, multilooked_array = utils.load_isce2_image(f'{image_path}.multilooked') + assert multilooked_array.shape == (20, 20) + + golden_array = np.zeros(multilooked_array.shape, dtype=np.float32) + golden_array[4:16, 4:16] = 1 + assert np.all(multilooked_array == golden_array) + + +def test_multilook_radar_merge_inputs(tmpdir): + paths = [ + 'fine_interferogram/IW1/burst_01.int', + 'geom_reference/IW1/lat_01.rdr', + 'geom_reference/IW1/lon_01.rdr', + 'geom_reference/IW1/los_01.rdr', + ] + paths = [Path(tmpdir) / x for x in paths] + for path in paths: + make_test_image(str(path)) + + mock_position = burst.BurstPosition(100, 100, 20, 60, 20, 60, 0.1, 0.1, datetime(2021, 1, 1, 0, 0, 0)) + with patch('hyp3_isce2.burst.load_burst_position') as mock_load_burst_position: + mock_load_burst_position.return_value = mock_position + output = burst.multilook_radar_merge_inputs(1, 5, 2, base_dir=tmpdir) + + assert output.n_lines == 50 + assert output.n_samples == 20 + + multilooked = [x.parent / f'{x.stem}.multilooked{x.suffix}' for x in paths] + for file in multilooked: + assert file.exists() diff --git a/tests/test_merge_tops_bursts.py b/tests/test_merge_tops_bursts.py new file mode 100644 index 00000000..54ff6bb0 --- /dev/null +++ b/tests/test_merge_tops_bursts.py @@ -0,0 +1,442 @@ +"""Tests for hyp3_isce2.merge_tops_bursts module, use single quotes""" +from copy import deepcopy +from dataclasses import dataclass +from datetime import datetime +from pathlib import Path +from unittest.mock import patch + +import asf_search +import isceobj # noqa: I100 +import lxml.etree as ET +import numpy as np +import pytest +from osgeo import gdal, osr +from requests import Session + +import hyp3_isce2.burst as burst_utils +import hyp3_isce2.merge_tops_bursts as merge +from hyp3_isce2 import utils + + +def mock_asf_search_results( + slc_name: str, + subswath: str, + polarization: str, + burst_index: int, + burst_id: int, + path_number: int, +) -> asf_search.ASFSearchResults: + product = asf_search.ASFProduct() + product.umm = {'InputGranules': [slc_name]} + product.properties.update( + { + 'burst': {'subswath': subswath, 'burstIndex': burst_index, 'relativeBurstID': burst_id}, + 'polarization': polarization, + 'url': f'https://foo.com/{slc_name}/baz.zip', + 'startTime': '2020-06-04T02:23:13.963847Z', + 'pathNumber': path_number, + } + ) + results = asf_search.ASFSearchResults([product]) + results.searchComplete = True + return results + + +def create_test_geotiff(output_file, dtype='float', n_bands=1): + """Create a test geotiff for testing""" + opts = {'float': (np.float64, gdal.GDT_Float64), 'cfloat': (np.complex64, gdal.GDT_CFloat32)} + np_dtype, gdal_dtype = opts[dtype] + data = np.ones((10, 10), dtype=np_dtype) + geotransform = [0.0, 1.0, 0.0, 0.0, 0.0, 1.0] + driver = gdal.GetDriverByName('GTiff') + dataset = driver.Create(output_file, 10, 10, n_bands, gdal_dtype) + dataset.SetGeoTransform(geotransform) + srs = osr.SpatialReference() + srs.ImportFromEPSG(4326) + dataset.SetProjection(srs.ExportToWkt()) + for i in range(n_bands): + band = dataset.GetRasterBand(i + 1) + band.WriteArray(data) + dataset = None + + +def test_to_burst_params(burst_product1): + assert burst_product1.to_burst_params() == burst_utils.BurstParams( + 'S1A_IW_SLC__1SDV_20200604T022251_20200604T022318_032861_03CE65_7C85', 'IW2', 'VV', 7 + ) + + +def test_get_burst_metadata(test_merge_dir, burst_product1): + product_path = list(test_merge_dir.glob('S1_136231*'))[0] + + with patch('hyp3_isce2.merge_tops_bursts.asf_search.granule_search') as mock_search: + mock_search.return_value = mock_asf_search_results( + 'S1A_IW_SLC__1SDV_20200604T022251_20200604T022318_032861_03CE65_7C85', 'IW2', 'VV', 7, 136231, 64 + ) + product = merge.get_burst_metadata([product_path])[0] + + assert product == burst_product1 + + +def test_prep_metadata_dirs(tmp_path): + annotation_dir, manifest_dir = merge.prep_metadata_dirs(tmp_path) + assert annotation_dir.is_dir() + assert manifest_dir.is_dir() + + +def test_download_metadata_xmls(monkeypatch, tmp_path, test_data_dir): + params = [burst_utils.BurstParams('foo', 'IW1', 'VV', 1), burst_utils.BurstParams('foo', 'IW2', 'VV', 0)] + sample_xml = ET.parse(test_data_dir / 'reference_descending.xml').getroot() + + with patch('hyp3_isce2.merge_tops_bursts.burst_utils.get_asf_session') as mock_session: + with patch('hyp3_isce2.merge_tops_bursts.burst_utils.download_metadata') as mock_download: + mock_session.return_value = Session() + mock_download.return_value = sample_xml + + merge.download_metadata_xmls(params, tmp_path) + + assert mock_download.call_count == 1 + assert len(list((tmp_path / 'annotation').glob('*.xml'))) == 2 + assert (tmp_path / 'manifest' / 'foo.xml').exists() + + +def test_get_scene_roi(test_s1_obj): + bursts = test_s1_obj.product.bursts + roi = merge.get_scene_roi(bursts) + golden_roi = (53.045079513806, 27.325111859227817, 54.15684468161031, 27.847161580403135) + assert np.all(np.isclose(roi, golden_roi)) + + +def test_load_isce_s1_obj(annotation_manifest_dirs): + annotation_dir, manifest_dir = annotation_manifest_dirs + s1_obj = merge.load_isce_s1_obj(2, 'VV', annotation_dir.parent) + + assert isinstance(s1_obj, merge.Sentinel1BurstSelect) + assert s1_obj.swath == 2 + assert s1_obj.polarization == 'vv' + assert len(s1_obj.tiff) == 1 + assert s1_obj.tiff[0] == '' + + +def test_Sentinel1BurstSelect(annotation_manifest_dirs, tmp_path, burst_product1): + annotation_dir, manifest_dir = annotation_manifest_dirs + s1_obj = merge.load_isce_s1_obj(2, 'VV', annotation_dir.parent) + + # Test select_bursts + test1_obj = deepcopy(s1_obj) + test1_utc = [burst_product1.start_utc] + test1_obj.select_bursts(test1_utc) + assert len(test1_obj.product.bursts) == 1 + assert test1_obj.product.numberOfBursts == 1 + assert test1_obj.product.bursts[0].burstStartUTC == test1_utc[0] + + test2_obj = deepcopy(s1_obj) + test2_utc = [datetime(2020, 6, 4, 2, 22, 57, 414185), datetime(2020, 6, 4, 2, 22, 54, 655908)] + test2_obj.select_bursts(test2_utc) + assert len(test2_obj.product.bursts) == 2 + assert test2_obj.product.bursts[0].burstStartUTC < test2_obj.product.bursts[1].burstStartUTC + + # Test update_burst_properties + test3_obj = deepcopy(test1_obj) + outpath = tmp_path / 'IW2' + test3_obj.output = str(outpath) + test3_obj.update_burst_properties([burst_product1]) + assert test3_obj.product.bursts[0].burstNumber == 1 + assert test3_obj.product.bursts[0].firstValidLine == 8 + assert test3_obj.product.bursts[0].numValidLines == 363 + assert test3_obj.product.bursts[0].firstValidSample == 9 + assert test3_obj.product.bursts[0].numValidSamples == 1220 + assert Path(test3_obj.product.bursts[0].image.filename).name == 'burst_01.slc' + + test4_obj = deepcopy(test1_obj) + bad_product = deepcopy(burst_product1) + bad_product.start_utc = datetime(1999, 1, 1, 1, 0, 0, 0) + with pytest.raises(ValueError, match='.*do not match.*'): + test4_obj.update_burst_properties([bad_product]) + + # Test write_xml + test5_obj = deepcopy(test3_obj) + test5_obj.write_xml() + assert outpath.with_suffix('.xml').exists() + + +def test_create_burst_cropped_s1_obj(annotation_manifest_dirs, burst_product1): + s1_obj = merge.create_burst_cropped_s1_obj(2, [burst_product1], 'VV', base_dir=annotation_manifest_dirs[0].parent) + assert isinstance(s1_obj, merge.Sentinel1BurstSelect) + assert Path(s1_obj.output).with_suffix('.xml').exists() + + +def test_modify_for_multilook(annotation_manifest_dirs, burst_product1): + burst_product = burst_product1 + s1_obj = merge.create_burst_cropped_s1_obj(2, [burst_product], 'VV', base_dir=annotation_manifest_dirs[0].parent) + + pre_burst = s1_obj.product.bursts[0] + assert not pre_burst.numberOfSamples == burst_product.n_samples + assert not pre_burst.numberOfLines == burst_product.n_lines + + burst_product.isce2_burst_number = s1_obj.product.bursts[0].burstNumber + looked_obj = merge.modify_for_multilook([burst_product], s1_obj) + burst = looked_obj.product.bursts[0] + assert burst.numberOfSamples == burst_product.n_samples + assert burst.numberOfLines == burst_product.n_lines + assert burst.firstValidSample == burst_product.first_valid_sample + assert burst.numValidSamples == burst_product.n_valid_samples + assert burst.firstValidLine == burst_product.first_valid_line + assert burst.numValidLines == burst_product.n_valid_lines + assert burst.sensingStop == burst_product.stop_utc + assert burst.azimuthTimeInterval == burst_product.az_time_interval + assert burst.rangePixelSize == burst_product.rg_pixel_size + assert looked_obj.output == 'fine_interferogram/IW2_multilooked' + + bad_product = deepcopy(burst_product) + bad_product.start_utc = datetime(1999, 1, 1, 1, 0, 0, 0) + with pytest.raises(ValueError, match='.*do not match.*'): + looked_obj = merge.modify_for_multilook([bad_product], s1_obj) + + +def test_download_dem_for_multiple_bursts(annotation_manifest_dirs, burst_product1): + base_dir = annotation_manifest_dirs[0].parent + s1_obj = merge.create_burst_cropped_s1_obj(2, [burst_product1], 'VV', base_dir=base_dir) + with patch('hyp3_isce2.merge_tops_bursts.download_dem_for_isce2') as mock_download: + mock_download.return_value = None + merge.download_dem_for_multiple_bursts([s1_obj]) + assert mock_download.call_count == 1 + assert isinstance(mock_download.call_args[0][0], tuple) + assert len(mock_download.call_args[0][0]) == 4 + assert mock_download.call_args[1]['dem_name'] == 'glo_30' + + +@pytest.mark.parametrize('isce_type,dtype,n_bands', [['ifg', 'cfloat', 1], ['lat', 'float', 1], ['los', 'float', 2]]) +def test_translate_image(isce_type, dtype, n_bands, tmp_path): + test_tiff = tmp_path / 'test.tif' + create_test_geotiff(str(test_tiff), dtype, n_bands) + out_path = tmp_path / 'test.bin' + merge.translate_image(str(test_tiff), str(out_path), isce_type) + for ext in ['.xml', '.vrt', '']: + assert (out_path.parent / (out_path.name + ext)).exists() + + opts = {'float': np.float32, 'cfloat': np.complex64} + image, array = utils.load_isce2_image(str(out_path)) + assert np.all(array == np.ones((10, 10), dtype=opts[dtype])) + + +def test_spoof_isce2_setup(annotation_manifest_dirs, burst_product1): + tmp_product = deepcopy(burst_product1) + tmp_product.isce2_burst_number = 1 + base_dir = annotation_manifest_dirs[0].parent + s1_obj = merge.create_burst_cropped_s1_obj(2, [tmp_product], 'VV', base_dir=base_dir) + merge.spoof_isce2_setup([tmp_product], s1_obj, base_dir=base_dir) + + fine_ifg_dir = base_dir / 'fine_interferogram' / 'IW2' + assert fine_ifg_dir.is_dir() + assert len(list(fine_ifg_dir.glob('*'))) == 3 + + geom_ref_dir = base_dir / 'geom_reference' / 'IW2' + assert geom_ref_dir.is_dir() + assert len(list(geom_ref_dir.glob('*'))) == 9 + + +def test_get_swath_list(tmp_path): + test_dir = tmp_path / 'test' + test_dir.mkdir() + assert merge.get_swath_list(str(test_dir)) == [] + + for x in [1, 2, 3]: + (test_dir / f'IW{x}').mkdir() + assert merge.get_swath_list(str(test_dir)) == [1, 2, 3] + + +def test_get_merged_orbit(test_s1_obj): + merged_orbit = merge.get_merged_orbit([test_s1_obj.product]) + assert isinstance(merged_orbit, isceobj.Orbit.Orbit.Orbit) + assert len(merged_orbit.stateVectors) == 17 # This number will change if the test data changes + + +def test_get_frames_and_indexes(isce2_merge_setup): + frames, burst_index = merge.get_frames_and_indexes(isce2_merge_setup / 'fine_interferogram') + assert len(frames) == 1 + assert isinstance(frames[0], isceobj.Sensor.TOPS.TOPSSwathSLCProduct.TOPSSwathSLCProduct) + assert burst_index[0] == [2, 0, 2] + + +# FIX: test_merge_bursts doesn't work due to pathing issue. +# def test_merge_bursts(isce2_merge_setup): +# import os +# os.chdir(isce2_merge_setup) +# merge.merge_bursts(20, 4) + + +def test_goldstein_werner_filter(tmp_path): + in_path = tmp_path / 'test.bin' + coh_path = tmp_path / 'coh.bin' + out_path = tmp_path / 'filtered.bin' + array = np.ones((10, 10), dtype=np.complex64) + utils.write_isce2_image(str(in_path), array) + utils.write_isce2_image(str(coh_path), array.astype(np.float32)) + merge.goldstein_werner_filter(str(in_path), str(out_path), str(coh_path)) + assert out_path.exists() + assert (out_path.parent / f'{out_path.name}.xml').exists() + assert (out_path.parent / f'{out_path.name}.vrt').exists() + + +def test_get_product_name(burst_product1): + product_name = merge.get_product_name(burst_product1, 80) + assert len(product_name) == 39 + assert product_name[:-4] == 'S1_064__20200604_20200616_VV_INT80_' + + +def test_make_parameter_file(test_data_dir, test_merge_dir, test_s1_obj, tmp_path): + ifg_dir = tmp_path / 'fine_interferogram' / 'IW2' + ifg_dir.mkdir(parents=True) + test_s1_obj.output = str(ifg_dir.parent / 'IW2') + test_s1_obj.write_xml() + + metas = merge.get_product_metadata_info(test_merge_dir) + + out_file = tmp_path / 'test.txt' + merge.make_parameter_file(out_file, metas, 20, 4, 0.6, True, base_dir=tmp_path) + assert out_file.exists() + + meta = utils.read_product_metadata(out_file) + assert len(meta['ReferenceGranule'].split(',')) == 2 + assert len(meta['SecondaryGranule'].split(',')) == 2 + assert meta['Rangelooks'] == '20' + assert meta['Azimuthlooks'] == '4' + with pytest.raises(KeyError): + assert meta['Radarnlines'] + + +def test_snaphu_unwrap(test_s1_obj, test_data_dir, tmp_path): + merge_dir = tmp_path / 'merged' + merge_dir.mkdir() + ifg_dir = tmp_path / 'fine_interferogram' / 'IW2' + ifg_dir.mkdir(parents=True, exist_ok=True) + test_s1_obj.output = str(ifg_dir.parent / 'IW2_multilooked') + test_s1_obj.write_xml() + + filt_path = merge_dir / 'filt_topophase.flat' + coh_path = merge_dir / 'coh.bin' + array = np.ones((100, 100), dtype=np.complex64) + utils.write_isce2_image(str(filt_path), array) + utils.write_isce2_image(str(coh_path), array.astype(np.float32)) + merge.snaphu_unwrap(2, 2, str(coh_path), base_dir=merge_dir) + + assert (merge_dir / 'filt_topophase.unw').exists() + assert (merge_dir / 'filt_topophase.unw.xml').exists() + assert (merge_dir / 'filt_topophase.unw.vrt').exists() + + +def test_geocode_products(test_data_dir, tmp_path, test_s1_obj): + merge_dir = tmp_path / 'merged' + merge_dir.mkdir() + ifg_dir = tmp_path / 'fine_interferogram' / 'IW2' + ifg_dir.mkdir(parents=True, exist_ok=True) + test_s1_obj.output = str(ifg_dir.parent / 'IW2') + test_s1_obj.write_xml() + + unw_path = merge_dir / 'filt_topophase.unw' + dem_path = merge_dir / 'dem.bin' + array = np.ones((377, 1272), dtype=np.float32) + utils.write_isce2_image(str(unw_path), array) + utils.write_isce2_image(str(dem_path), array) + merge.geocode_products(1, 1, str(dem_path), base_dir=merge_dir, to_be_geocoded=[str(unw_path)]) + + assert (merge_dir / 'filt_topophase.unw.geo').exists() + assert (merge_dir / 'filt_topophase.unw.geo.xml').exists() + assert (merge_dir / 'filt_topophase.unw.geo.vrt').exists() + + +def test_check_burst_group_validity(): + @dataclass + class Product: + reference_date: datetime + secondary_date: datetime + polarization: str + relative_orbit: int + swath: int + burst_id: int + range_looks: int + azimuth_looks: int + + # Test valid products + good_products = [ + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 1, 111116, 20, 4), + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 1, 111117, 20, 4), + ] + merge.check_burst_group_validity(good_products) + + # Bad polarization + bad_pol = deepcopy(good_products) + bad_pol[0].polarization = 'HH' + with pytest.raises(ValueError, match='All products.*polarization.*'): + merge.check_burst_group_validity(bad_pol) + + # Bad reference date + bad_ref_date = deepcopy(good_products) + bad_ref_date[1].reference_date = datetime(2020, 2, 2) + with pytest.raises(ValueError, match='All products.*reference date.*'): + merge.check_burst_group_validity(bad_ref_date) + + # Bad reference date + bad_sec_date = deepcopy(good_products) + bad_sec_date[1].secondary_date = datetime(2020, 2, 2) + with pytest.raises(ValueError, match='All products.*secondary date.*'): + merge.check_burst_group_validity(bad_sec_date) + + # Bad relative orbit + bad_rel_orbit = deepcopy(good_products) + bad_rel_orbit[0].relative_orbit = 2 + with pytest.raises(ValueError, match='All products.*relative orbit.*'): + merge.check_burst_group_validity(bad_rel_orbit) + + # Bad looks + bad_range_looks = deepcopy(good_products) + bad_range_looks[1].range_looks = 10 + with pytest.raises(ValueError, match='All products.*looks.*'): + merge.check_burst_group_validity(bad_range_looks) + + # Non-contiguous bursts + non_contiguous_swath_products = [ + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 1, 111115, 20, 4), + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 1, 111117, 20, 4), + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 2, 111115, 20, 4), + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 2, 111116, 20, 4), + ] + with pytest.raises(ValueError, match='Products.*swath 1.*contiguous'): + merge.check_burst_group_validity(non_contiguous_swath_products) + + # Non-contiguous swath bursts + non_contiguous_swath_products = [ + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 1, 111116, 20, 4), + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 1, 111117, 20, 4), + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 2, 111114, 20, 4), + Product(datetime(2020, 2, 1), datetime(2020, 1, 1), 'VV', 1, 2, 111113, 20, 4), + ] + with pytest.raises(ValueError, match='Products.*swaths 1 and 2.*overlap'): + merge.check_burst_group_validity(non_contiguous_swath_products) + + +def test_get_product_multilook(tmp_path): + product_dir = tmp_path / 'test' + product_dir.mkdir() + product1 = product_dir / 'S1_111111_IW1_VV_01' + product1.mkdir() + metadata1 = product1 / 'S1_111111_IW1_VV_01.txt' + + metadata1.write_text('Rangelooks: 20\nAzimuthlooks: 4\n') + range_looks, azimuth_looks = merge.get_product_multilook(product_dir) + assert range_looks == 20 + assert azimuth_looks == 4 + + +def test_make_readme(tmp_path): + prod_name = 'foo' + tmp_prod_dir = tmp_path / prod_name + tmp_prod_dir.mkdir(exist_ok=True) + create_test_geotiff(str(tmp_prod_dir / f'{prod_name}_wrapped_phase.tif')) + reference_scenes = ['a_a_a_20200101T000000_a', 'b_b_b_20200101T000000_b'] + secondary_scenes = ['c_c_c_20210101T000000_c', 'd_d_d_20210101T000000_d'] + + merge.make_readme(tmp_prod_dir, reference_scenes, secondary_scenes, 2, 10, True) + out_path = tmp_prod_dir / f'{prod_name}_README.md.txt' + assert out_path.exists() diff --git a/tests/test_utils.py b/tests/test_utils.py index 236735ec..fd4af822 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,37 +1,35 @@ +import filecmp import os +import shutil +from pathlib import Path +from unittest.mock import patch import numpy as np import pytest from osgeo import gdal -from hyp3_isce2.utils import ( - ESA_HOST, - GDALConfigManager, - extent_from_geotransform, - get_esa_credentials, - make_browse_image, - oldest_granule_first, - resample_to_radar, - utm_from_lon_lat, -) +import hyp3_isce2.utils as utils + +import isceobj # noqa + gdal.UseExceptions() def test_utm_from_lon_lat(): - assert utm_from_lon_lat(0, 0) == 32631 - assert utm_from_lon_lat(-179, -1) == 32701 - assert utm_from_lon_lat(179, 1) == 32660 - assert utm_from_lon_lat(27, 89) == 32635 - assert utm_from_lon_lat(182, 1) == 32601 - assert utm_from_lon_lat(-182, 1) == 32660 - assert utm_from_lon_lat(-360, -1) == 32731 + assert utils.utm_from_lon_lat(0, 0) == 32631 + assert utils.utm_from_lon_lat(-179, -1) == 32701 + assert utils.utm_from_lon_lat(179, 1) == 32660 + assert utils.utm_from_lon_lat(27, 89) == 32635 + assert utils.utm_from_lon_lat(182, 1) == 32601 + assert utils.utm_from_lon_lat(-182, 1) == 32660 + assert utils.utm_from_lon_lat(-360, -1) == 32731 def test_extent_from_geotransform(): - assert extent_from_geotransform((0, 1, 0, 0, 0, -1), 1, 1) == (0, 0, 1, -1) - assert extent_from_geotransform((0, 1, 0, 0, 0, -1), 2, 2) == (0, 0, 2, -2) - assert extent_from_geotransform((0, 1, 0, 0, 0, -1), 1, 3) == (0, 0, 1, -3) + assert utils.extent_from_geotransform((0, 1, 0, 0, 0, -1), 1, 1) == (0, 0, 1, -1) + assert utils.extent_from_geotransform((0, 1, 0, 0, 0, -1), 2, 2) == (0, 0, 2, -2) + assert utils.extent_from_geotransform((0, 1, 0, 0, 0, -1), 1, 3) == (0, 0, 1, -3) def test_gdal_config_manager(): @@ -43,7 +41,7 @@ def test_gdal_config_manager(): assert gdal.GetConfigOption('OPTION3') is None assert gdal.GetConfigOption('OPTION4') is None - with GDALConfigManager(OPTION2='CHANGED', OPTION3='VALUE3'): + with utils.GDALConfigManager(OPTION2='CHANGED', OPTION3='VALUE3'): assert gdal.GetConfigOption('OPTION1') == 'VALUE1' assert gdal.GetConfigOption('OPTION2') == 'CHANGED' assert gdal.GetConfigOption('OPTION3') == 'VALUE3' @@ -58,17 +56,17 @@ def test_gdal_config_manager(): def test_oldest_granule_first(): - oldest = "S1_249434_IW1_20230511T170732_VV_07DE-BURST" - latest = "S1_249434_IW1_20230523T170733_VV_8850-BURST" - assert oldest_granule_first(oldest, latest) == (oldest, latest) - assert oldest_granule_first(latest, oldest) == (oldest, latest) + oldest = 'S1_249434_IW1_20230511T170732_VV_07DE-BURST' + latest = 'S1_249434_IW1_20230523T170733_VV_8850-BURST' + assert utils.oldest_granule_first(oldest, latest) == (oldest, latest) + assert utils.oldest_granule_first(latest, oldest) == (oldest, latest) def test_make_browse_image(): - input_tif = "tests/data/test_geotiff.tif" - output_png = "tests/data/test_browse_image2.png" - make_browse_image(input_tif, output_png) - assert open(output_png, "rb").read() == open("tests/data/test_browse_image.png", "rb").read() + input_tif = 'tests/data/test_geotiff.tif' + output_png = 'tests/data/test_browse_image2.png' + utils.make_browse_image(input_tif, output_png) + assert open(output_png, 'rb').read() == open('tests/data/test_browse_image.png', 'rb').read() os.remove(output_png) @@ -89,7 +87,7 @@ def check_correctness_of_resample(mask, lat, lon, geotransform, data_type, outsh lat[row, col] = y + row * mask_y_res lon[row, col] = x + col * mask_x_res - resampled_image = resample_to_radar(mask, lat, lon, geotransform, data_type, outshape) + resampled_image = utils.resample_to_radar(mask, lat, lon, geotransform, data_type, outshape) lon_lat_complex = lon + 1j * lat @@ -129,14 +127,39 @@ def test_resample_to_radar(): resample_with_different_case(30, 10, 10, 10, geotransform) +def test_resample_to_radar_io(tmp_path, test_merge_dir): + out_paths = [] + array = np.ones((10, 10), dtype=np.float32) + for image_name in ['input', 'lat', 'lon']: + out_path = str(tmp_path / image_name) + image = utils.create_image(str(out_path), 10, access_mode='write') + image.coord1.coordSize = 10 + image.coord2.coordSize = 10 + utils.write_isce2_image_from_obj(image, array) + out_paths.append(out_path) + + image_to_resample, latin, lonin = out_paths + output = str(tmp_path / 'output') + + with patch('hyp3_isce2.utils.resample_to_radar') as mock_resample: + mock_resample.return_value = array + utils.resample_to_radar_io(image_to_resample, latin, lonin, output) + + assert Path(output).is_file() + + latim, lat = utils.load_isce2_image(latin) + outputim, outputarray = utils.load_isce2_image(output) + assert np.all(outputarray == array) + + def test_get_esa_credentials_env(tmp_path, monkeypatch): with monkeypatch.context() as m: m.setenv('ESA_USERNAME', 'foo') m.setenv('ESA_PASSWORD', 'bar') m.setenv('HOME', str(tmp_path)) - (tmp_path / '.netrc').write_text(f'machine {ESA_HOST} login netrc_username password netrc_password') + (tmp_path / '.netrc').write_text(f'machine {utils.ESA_HOST} login netrc_username password netrc_password') - username, password = get_esa_credentials() + username, password = utils.get_esa_credentials() assert username == 'foo' assert password == 'bar' @@ -146,9 +169,9 @@ def test_get_esa_credentials_netrc(tmp_path, monkeypatch): m.delenv('ESA_USERNAME', raising=False) m.delenv('ESA_PASSWORD', raising=False) m.setenv('HOME', str(tmp_path)) - (tmp_path / '.netrc').write_text(f'machine {ESA_HOST} login foo password bar') + (tmp_path / '.netrc').write_text(f'machine {utils.ESA_HOST} login foo password bar') - username, password = get_esa_credentials() + username, password = utils.get_esa_credentials() assert username == 'foo' assert password == 'bar' @@ -161,7 +184,7 @@ def test_get_esa_credentials_missing(tmp_path, monkeypatch): (tmp_path / '.netrc').write_text('') msg = 'Please provide.*' with pytest.raises(ValueError, match=msg): - get_esa_credentials() + utils.get_esa_credentials() with monkeypatch.context() as m: m.setenv('ESA_USERNAME', 'env_username') @@ -170,4 +193,165 @@ def test_get_esa_credentials_missing(tmp_path, monkeypatch): (tmp_path / '.netrc').write_text('') msg = 'Please provide.*' with pytest.raises(ValueError, match=msg): - get_esa_credentials() + utils.get_esa_credentials() + + +def test_create_image(tmp_path): + def _check_create_image(path: str, image_subtype: str = 'default'): + # create an isce image (includes binary, .vrt, and .xml files) + array = np.array(range(150), dtype=np.float32) + array = array.reshape(15, 10) + length, width = array.shape + out_path = str(tmp_path / 'isce_image_2d') + utils.write_isce2_image(out_path, array=array) + + # test ifg in create, finalize, and load modes + path_c = path + '/img_via_create' + img_c = utils.create_image( + path_c, width=width, access_mode='write', image_subtype=image_subtype, action='create' + ) + assert Path(img_c.getFilename()).is_file() + + path_f = path + '/img_via_finalize' + shutil.copy(out_path, path_f) + img_f = utils.create_image( + path_f, width=width, access_mode='read', image_subtype=image_subtype, action='finalize' + ) + assert Path(img_f.getFilename()).is_file() + assert Path(img_f.getFilename() + '.vrt').is_file() + assert Path(img_f.getFilename() + '.xml').is_file() + + path_l = path + '/img_via_load' + shutil.copy(out_path, path_l) + shutil.copy(f'{out_path}.vrt', f'{path_l}.vrt') + shutil.copy(f'{out_path}.xml', f'{path_l}.xml') + + img_l = utils.create_image(path_l, access_mode='load', image_subtype=image_subtype, action='load') + assert Path(img_l.getFilename()).is_file() + assert Path(img_f.getFilename() + '.vrt').is_file() + assert Path(img_f.getFilename() + '.xml').is_file() + + _check_create_image(str(tmp_path), image_subtype='ifg') + _check_create_image(str(tmp_path), image_subtype='cor') + _check_create_image(str(tmp_path), image_subtype='unw') + _check_create_image(str(tmp_path), image_subtype='conncomp') + _check_create_image(str(tmp_path), image_subtype='default') + + +def test_write_isce2_image(tmp_path): + array = np.array(range(150), dtype=np.float32) + array = array.reshape(15, 10) + out_path = str(tmp_path / 'isce_image_2d') + utils.write_isce2_image(out_path, array=array) + assert Path(out_path).is_file() + + +def test_write_isce2_image_from_obj(tmp_path): + def _check_write_isce2_image_from_obj(out_path, bands, length, width): + image = isceobj.createImage() + image.initImage(out_path, 'write', width, 'FLOAT', bands) + image.setLength(length) + image.setImageType('bil') + image.createImage() + utils.write_isce2_image_from_obj(image, array=array) + assert Path(image.filename).is_file() + assert Path(image.filename + '.vrt').is_file() + assert Path(image.filename + '.xml').is_file() + + # 1D array, it shape(width), band=1, length=1 + array = np.array(range(150), dtype=np.float32) + bands = 1 + length = 1 + width = array.shape[0] + out_path = str(tmp_path / 'isce_image_1d') + _check_write_isce2_image_from_obj(out_path, bands, length, width) + + # 2D array, is shape(length, width), band=1 + array = np.array(range(150), dtype=np.float32) + array = array.reshape(15, 10) + bands = 1 + length, width = array.shape + out_path = str(tmp_path / 'isce_image_2d') + _check_write_isce2_image_from_obj(out_path, bands, length, width) + + # multi-D array, its shape(band,length, width) + array = np.array(range(150), dtype=np.float32) + array = array.reshape(3, 5, 10) + bands, length, width = array.shape + out_path = str(tmp_path / 'isce_image_md') + _check_write_isce2_image_from_obj(out_path, bands, length, width) + + +def test_load_isce2_image(tmp_path): + in_path = str(tmp_path / 'isce_image_md') + arrayin = np.array(range(150), dtype=np.float32) + arrayin = arrayin.reshape(3, 5, 10) + utils.write_isce2_image(in_path, array=arrayin) + image_obj, arrayout = utils.load_isce2_image(in_path) + assert isinstance(image_obj, isceobj.Image.Image.Image) + assert np.array_equal(arrayin, arrayout) + + in_path = str(tmp_path / 'isce_image_2d') + + arrayin = np.array(range(150), dtype=np.float32) + arrayin = arrayin.reshape(15, 10) + utils.write_isce2_image(in_path, array=arrayin) + image_obj, arrayout = utils.load_isce2_image(in_path) + assert isinstance(image_obj, isceobj.Image.Image.Image) + assert np.array_equal(arrayin, arrayout) + + in_path = str(tmp_path / 'isce_image_1d') + arrayin = np.array(range(150), dtype=np.float32) + arrayin = arrayin.reshape(1, 150) + utils.write_isce2_image(in_path, array=arrayin) + image_obj, arrayout = utils.load_isce2_image(in_path) + assert isinstance(image_obj, isceobj.Image.Image.Image) + assert np.array_equal(arrayin, arrayout) + + +def test_get_geotransform_from_dataset(tmp_path): + image_path = tmp_path / 'test' + image = utils.create_image(str(image_path), 100, access_mode='write') + + image.coord1.coordStart = 1 + image.coord1.coordDelta = 10 + image.coord2.coordStart = 11 + image.coord2.coordDelta = 100 + assert utils.get_geotransform_from_dataset(image) == (1, 10, 0, 11, 0, 100) + + +def test_isce2_copy(tmp_path): + in_path = str(tmp_path / 'isce_image_2d') + arrayin = np.array(range(150), dtype=np.float32) + arrayin = arrayin.reshape(15, 10) + utils.write_isce2_image(in_path, array=arrayin) + out_path = str(tmp_path / 'isce_image_2d_copy') + utils.isce2_copy(in_path, out_path) + assert filecmp.cmp(in_path, out_path) + + +def test_image_math(tmp_path): + in_path1 = str(tmp_path / 'isce_image_2d_1') + array1 = np.array(range(150), dtype=np.float32) + array1 = array1.reshape(15, 10) + utils.write_isce2_image(in_path1, array=array1) + + in_path2 = str(tmp_path / 'isce_image_2d_2') + array2 = np.array(range(3, 153), dtype=np.float32) + array2 = array2.reshape(15, 10) + utils.write_isce2_image(in_path2, array=array2) + + out_path = str(tmp_path / 'isce_image_2d_out') + utils.image_math(in_path1, in_path2, out_path, 'a + b') + image_obj_out, arrayout = utils.load_isce2_image(out_path) + assert np.array_equal(array1 + array2, arrayout) + + +def test_read_product_metadata(test_merge_dir): + metafile = list(test_merge_dir.glob('S1_136232*/*.txt'))[0] + metas = utils.read_product_metadata(metafile) + assert metas['ReferenceGranule'] == 'S1_136232_IW2_20200604T022315_VV_7C85-BURST' + assert metas['SecondaryGranule'] == 'S1_136232_IW2_20200616T022316_VV_5D11-BURST' + assert np.isclose(float(metas['Baseline']), -66.10716474087386) + assert int(metas['ReferenceOrbitNumber']) == 32861 + assert int(metas['SecondaryOrbitNumber']) == 33036