From 6351b46fc059ca7a90c8caea80075f251ccab6e6 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Mon, 24 Feb 2025 09:52:35 -0600 Subject: [PATCH 1/5] add script for converting outputs to geotiffs --- src/hyp3_srg/convert_srg_outputs.py | 182 ++++++++++++++++++++++++++++ 1 file changed, 182 insertions(+) create mode 100644 src/hyp3_srg/convert_srg_outputs.py diff --git a/src/hyp3_srg/convert_srg_outputs.py b/src/hyp3_srg/convert_srg_outputs.py new file mode 100644 index 0000000..b48d2f3 --- /dev/null +++ b/src/hyp3_srg/convert_srg_outputs.py @@ -0,0 +1,182 @@ +import argparse +from pathlib import Path + +import numpy as np +from osgeo import gdal + + +gdal.UseExceptions() + + +def parse_rsc(rsc_file: Path) -> dict: + """Parse the metadata from an SRG RSC file + + Args: + rsc_file: Path to the RSC file + + Returns: + Dictionary with the metadata from the RSC file + """ + types_dict = { + 'WIDTH': int, + 'FILE_LENGTH': int, + 'X_FIRST': float, + 'Y_FIRST': float, + 'X_STEP': float, + 'Y_STEP': float, + 'X_UNIT': str, + 'Y_UNIT': str, + 'Z_OFFSET': float, + 'Z_SCALE': float, + 'PROJECTION': str, + 'xstart': float, + 'ystart': float, + 'xend': float, + 'yend': float, + 'xsize': int, + 'ysize': int, + } + with open(rsc_file, 'r') as f: + lines = f.readlines() + rsc_dict = {} + for line in lines: + key, value = line.strip().split() + rsc_dict[key] = types_dict[key](value) + return rsc_dict + + +def parse_sbas_list(sbas_list: Path) -> tuple[list[str], list[list[str]]]: + """Parse the list of SBAS pairs from a SBAS list file + + Args: + sbas_list: Path to the SBAS list file + + Returns: + List with the names of the files in the SBAS list + List with the pairs of files in the SBAS list + """ + with open(sbas_list, 'r') as f: + lines = f.readlines() + pairs = [[Path(x).name for x in line.strip().split()[:2]] for line in lines] + files = sorted(list(set([x for pair in pairs for x in pair]))) + return files, pairs + + +def write_geotiff(data: np.ndarray, info: dict, filename: Path, band_names: list[str] = None) -> None: + """Write a numpy array to a GeoTIFF file + + Args: + data: Numpy array with the data to be written + info: Dictionary with the metadata from the DEM RSC file + filename: Path to the output GeoTIFF file + band_names: List with the names of the bands in the output GeoTIFF + """ + if data.dtype == np.float32: + dtype = gdal.GDT_Float32 + elif data.dtype == np.int16: + dtype = gdal.GDT_Int16 + elif data.dtype == np.int32: + dtype = gdal.GDT_Int32 + else: + raise ValueError('Data type not supported') + + if len(data.shape) == 3: + n_bands, n_rows, n_cols = data.shape + else: + n_bands = 1 + n_rows, n_cols = data.shape + + driver = gdal.GetDriverByName('GTiff') + ds = driver.Create(str(filename), n_cols, n_rows, n_bands, dtype) + ds.SetGeoTransform((info['X_FIRST'], info['X_STEP'], 0, info['Y_FIRST'], 0, info['Y_STEP'])) + ds.SetProjection('EPSG:4326') + if n_bands == 1: + ds.GetRasterBand(1).WriteArray(data) + else: + for i in range(n_bands): + band = ds.GetRasterBand(i + 1) + band.WriteArray(data[i, :, :]) + if band_names is not None: + band.SetDescription(band_names[i]) + band.FlushCache() + del band + ds = None + + +def convert_single_band(file: Path, info: dict, dtype: np.dtype) -> None: + """Convert a single-band SRG file to a GeoTIFF + + Args: + file: Path to the SRG file + info: Dictionary with the metadata from the DEM RSC file + dtype: Numpy data type of the file + """ + data = np.fromfile(file, dtype=dtype).reshape(info['FILE_LENGTH'], -1) + write_geotiff(data, info, file.with_suffix('.tif')) + + +def convert_displacement(displacement_file: Path, info: dict, band_names: list[str]) -> None: + """Convert the velocity file to separate GeoTIFFs for amplitude and displacement + Save each date as a separate band + + Args: + displacement_file: Path to the SRG displacement file + info: Dictionary with the metadata from the DEM RSC file + band_names: List with the dates of the SBAS pairs + """ + data = np.fromfile(displacement_file, dtype=np.float32).reshape(-1, info['WIDTH']) + + amplitude = data[0::2, :] + amplitude = amplitude.reshape(len(band_names), info['FILE_LENGTH'], info['WIDTH']) + write_geotiff(amplitude, info, displacement_file.parent / 'amplitude.tif', band_names=band_names) + + displacement = data[1::2, :] + displacement = displacement.reshape(len(band_names), info['FILE_LENGTH'], info['WIDTH']) + write_geotiff(displacement, info, displacement_file.parent / 'displacement.tif', band_names=band_names) + + +def convert_velocity(velocity_file: Path, info: dict, band_names: list[str]) -> None: + """Convert the velocity file to a GeoTIFF + Save each date as a separate band + + Args: + velocity_file: Path to the SRG velocity file + info: Dictionary with the metadata from the DEM RSC file + band_names: List with the dates of the SBAS pairs + """ + data = np.fromfile(velocity_file, dtype=np.float32).reshape(-1, info['WIDTH']) + velocity = data.reshape(len(band_names), info['FILE_LENGTH'], info['WIDTH']) + write_geotiff(velocity, info, velocity_file.parent / 'velocity.tif', band_names=band_names) + + +def convert_files(sbas_dir: Path) -> None: + """Convert the contents of a SRG SBAS directory to GeoTIFFs + + Args: + sbas_dir: Path to the SBAS directory for a SRG run + """ + info = parse_rsc(sbas_dir / 'dem.rsc') + files, pairs = parse_sbas_list(sbas_dir / 'sbas_list') + band_names = [x[17:25] for x in files][1:] + convert_single_band(sbas_dir / 'dem', info, dtype=np.int16) + convert_single_band(sbas_dir / 'npts', info, dtype=np.int32) + convert_single_band(sbas_dir / 'stacktime', info, dtype=np.int32) + convert_displacement(sbas_dir / 'displacement', info, band_names) + convert_velocity(sbas_dir / 'velocity', info, band_names) + + +def main(): + """CLI for the convert_srg_outputs.py script + + Example: + python convert_srg_outputs.py ./path/to/srg_sbas_directory + """ + parser = argparse.ArgumentParser(description='Convert the contents of a SRG SBAS directory to GeoTIFF') + parser.add_argument('directory', type=str, help='SRG SBAS directory') + args = parser.parse_args() + args.directory = Path(args.directory) + convert_files(args.directory) + + +if __name__ == '__main__': + main() From c9de6bec1858908443804499e9eac07ac012f17c Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Mon, 24 Feb 2025 09:53:26 -0600 Subject: [PATCH 2/5] update changelog --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 372aae0..1290625 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ 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). +## [0.9.4] + +### Added +* Utility that enables converting outputs to geotiffs. + ## [0.9.3] ### Removed From 2fe1f22c588f5c0709693c98bc22118fd790d67f Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Tue, 26 Aug 2025 15:19:29 -0500 Subject: [PATCH 3/5] move geotiff functionality --- {src/hyp3_srg => scripts}/convert_srg_outputs.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {src/hyp3_srg => scripts}/convert_srg_outputs.py (100%) diff --git a/src/hyp3_srg/convert_srg_outputs.py b/scripts/convert_srg_outputs.py similarity index 100% rename from src/hyp3_srg/convert_srg_outputs.py rename to scripts/convert_srg_outputs.py From 2e3bee48a0e8a850a47cd191aafab8466bb1145d Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Tue, 26 Aug 2025 16:40:31 -0500 Subject: [PATCH 4/5] create image viewer --- scripts/view_image.py | 102 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100755 scripts/view_image.py diff --git a/scripts/view_image.py b/scripts/view_image.py new file mode 100755 index 0000000..4e7c8b2 --- /dev/null +++ b/scripts/view_image.py @@ -0,0 +1,102 @@ +import argparse +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +from convert_srg_outputs import parse_rsc + + +def load_image(file_path: Path, info: dict): + suffix = file_path.suffix + if suffix == '.geo': + ftype = 'slc' + data = np.fromfile(file_path, dtype=np.float32).reshape(info['FILE_LENGTH'], -1) + real = data[:, 0::2] + imag = data[:, 1::2] + data = np.abs(real + 1j * imag) + elif suffix == '.amp': + ftype = 'amp' + data = np.fromfile(file_path, dtype=np.float32).reshape(info['FILE_LENGTH'], -1) + elif suffix == '.int': + ftype = 'int' + data = np.fromfile(file_path, dtype=np.float32).reshape(info['FILE_LENGTH'], -1) + real = data[:, 0::2] + imag = data[:, 1::2] + data = np.angle(real + 1j * imag) + elif suffix == '.unw': + ftype = 'unw' + data = np.fromfile(file_path, dtype=np.float32).reshape(info['FILE_LENGTH'], -1) + data = data[:, data.shape[1] // 2 :] + elif suffix == '.cc': + ftype = 'cc' + data = np.fromfile(file_path, dtype=np.float32).reshape(info['FILE_LENGTH'], -1) + data = data[:, data.shape[1] // 2 :] + elif suffix == '': + name = file_path.name + if name == 'dem': + ftype = 'dem' + data = np.fromfile(file_path, dtype=np.int16).reshape(info['FILE_LENGTH'], -1) + elif name == 'npts': + ftype = 'npts' + data = np.fromfile(file_path, dtype=np.int32).reshape(info['FILE_LENGTH'], -1) + elif name == 'stacktime': + ftype = 'stacktime' + data = np.fromfile(file_path, dtype=np.int32).reshape(info['FILE_LENGTH'], -1) + elif name == 'velocity': + raise NotImplementedError('Loading of velocity files not implemented.') + elif name == 'displacement': + raise NotImplementedError('Loading of displacement files not implemented.') + else: + raise ValueError(f'File {file_path} not recognized') + else: + raise ValueError(f'File {file_path} not recognized') + return ftype, data + + +def plot_image(file_path: Path, rsc_path: Path, file_path2: Path | None = None): + info = parse_rsc(rsc_path) + ftype, data = load_image(file_path, info) + title = file_path.name + + if file_path2 is not None: + ftype2, data2 = load_image(file_path2, info) + if ftype != ftype2: + raise ValueError('Files must be of the same type to plot the difference.') + data = data - data2 + title = f'Difference: {file_path.name} - {file_path2.name}' + + cmap = 'gray' + if ftype in ['amp', 'slc']: + data = 10 * np.log10(data) + if ftype == 'int': + cmap = 'hsv' + + f, ax = plt.subplots(1, 1, figsize=(10, 10)) + im = ax.imshow(data, cmap=cmap, vmin=float(np.percentile(data, 1)), vmax=float(np.percentile(data, 99))) + f.colorbar(im, ax=ax, shrink=0.5) + ax.set_title(title) + plt.tight_layout() + plt.show() + + +def main(): + """CLI for the convert_srg_outputs.py script + + Example: + python convert_srg_outputs.py ./path/to/srg_sbas_directory + """ + parser = argparse.ArgumentParser(description='Convert the contents of a SRG SBAS directory to GeoTIFF') + parser.add_argument('file', type=Path, help='Path to file to display.') + parser.add_argument('rsc', type=Path, help='Path to companion RSC file.') + parser.add_argument( + 'file2', + type=Path, + help='Path a second file. Image displayed will be the difference between the two files.', + nargs='?', + ) + args = parser.parse_args() + plot_image(args.file, args.rsc, args.file2) + + +if __name__ == '__main__': + main() From be66768d908ed070c10264097011709809921a83 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 4 Sep 2025 16:02:29 -0500 Subject: [PATCH 5/5] update viewing scripts --- scripts/view_image.py | 67 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 55 insertions(+), 12 deletions(-) diff --git a/scripts/view_image.py b/scripts/view_image.py index 4e7c8b2..4d7561a 100755 --- a/scripts/view_image.py +++ b/scripts/view_image.py @@ -53,25 +53,16 @@ def load_image(file_path: Path, info: dict): return ftype, data -def plot_image(file_path: Path, rsc_path: Path, file_path2: Path | None = None): +def plot_image(file_path: Path, rsc_path: Path): info = parse_rsc(rsc_path) ftype, data = load_image(file_path, info) title = file_path.name - if file_path2 is not None: - ftype2, data2 = load_image(file_path2, info) - if ftype != ftype2: - raise ValueError('Files must be of the same type to plot the difference.') - data = data - data2 - title = f'Difference: {file_path.name} - {file_path2.name}' - cmap = 'gray' - if ftype in ['amp', 'slc']: - data = 10 * np.log10(data) if ftype == 'int': cmap = 'hsv' - f, ax = plt.subplots(1, 1, figsize=(10, 10)) + f, ax = plt.subplots(1, 1, figsize=(15, 10)) im = ax.imshow(data, cmap=cmap, vmin=float(np.percentile(data, 1)), vmax=float(np.percentile(data, 99))) f.colorbar(im, ax=ax, shrink=0.5) ax.set_title(title) @@ -79,6 +70,55 @@ def plot_image(file_path: Path, rsc_path: Path, file_path2: Path | None = None): plt.show() +def plot_hist(ax, array, color, label): + ax.hist( + array.flatten(), + bins=50, + color=color, + label=label, + alpha=0.25, + range=(np.percentile(array, 1), np.percentile(array, 99)), + ) + + +def plot_diff(file_path: Path, file_path2: Path, rsc_path: Path, save: bool = False): + info = parse_rsc(rsc_path) + ftype, data = load_image(file_path, info) + title = file_path.name + + ftype2, data2 = load_image(file_path2, info) + if ftype != ftype2: + raise ValueError('Files must be of the same type to plot the difference.') + diff = np.abs(data - data2) + title = f'{file_path}\n-\n{file_path2}' + + cmap = 'gray' + if ftype == 'int': + cmap = 'hsv' + + f, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 10)) + im = ax1.imshow(diff, cmap=cmap, vmin=np.percentile(diff, 1), vmax=np.percentile(diff, 99)) + plot_hist(ax2, diff, 'black', 'Difference') + f.colorbar(im, ax=ax1, shrink=0.5) + ax1.set_title(title) + plt.tight_layout() + if not save: + plt.show() + else: + out_path = file_path.parent / f'{file_path.stem}_minus_{file_path2.stem}.png' + plt.savefig(out_path, dpi=300) + print(f'Saved figure to {out_path}') + plt.close(f) + + +def plots(): + image1 = Path('stanford/raw/S1A_IW_RAW__0SDV_20241006T161640_20241006T161717_055984_06D8A0_FC33.geo') + image2s = list(Path('hawaii_old').glob('*.geo')) + rsc = Path('stanford/elevation.dem.rsc') + for image2 in image2s: + plot_diff(image1, image2, rsc, save=True) + + def main(): """CLI for the convert_srg_outputs.py script @@ -95,7 +135,10 @@ def main(): nargs='?', ) args = parser.parse_args() - plot_image(args.file, args.rsc, args.file2) + if args.file2 is None: + plot_image(args.file, args.rsc) + else: + plot_diff(args.file, args.file2, args.rsc) if __name__ == '__main__':