diff --git a/CHANGELOG.md b/CHANGELOG.md index 281f3d5d..9fbf5616 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ 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.7.0] + +## Added +* `opera_rgb.py`, a workflow for creating OPERA RTC RGB decompositions. ## [0.6.0] diff --git a/environment.yml b/environment.yml index 7646b3e7..76188652 100644 --- a/environment.yml +++ b/environment.yml @@ -20,6 +20,7 @@ dependencies: - pytest-cov # For running - astropy + - asf_search - boto3 - fiona - gdal>=3.7 diff --git a/pyproject.toml b/pyproject.toml index 4a5f0959..73ebc51b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ classifiers=[ ] dependencies = [ "astropy", + "asf_search", "fiona", "gdal>=3.3", "numpy", @@ -41,6 +42,7 @@ make_composite = "asf_tools.composite:main" water_map = "asf_tools.hydrosar.water_map:main" calculate_hand = "asf_tools.hydrosar.hand.calculate:main" flood_map = "asf_tools.hydrosar.flood_map:main" +opera_rgb = "asf_tools.opera_rgb:main" [project.entry-points.hyp3] water_map = "asf_tools.hydrosar.water_map:hyp3" diff --git a/src/asf_tools/opera_rgb.py b/src/asf_tools/opera_rgb.py new file mode 100644 index 00000000..aa3c1be2 --- /dev/null +++ b/src/asf_tools/opera_rgb.py @@ -0,0 +1,158 @@ +"""Create a georefernced RGB decomposition RTC image from two input OPERA RTCs""" +import argparse +from pathlib import Path +from typing import Optional + +import asf_search +import numpy as np +from osgeo import gdal, gdalconst + + +gdal.UseExceptions() + + +def prep_data(granule: str): + """Download and prepare the data needed to create an RGB decomposition image. + + Args: + granule: OPERA granule name. + + Returns: + Tuple of co-pol, cross-pol, and mask filenames, None for each if not available. + """ + result = asf_search.granule_search([granule])[0] + urls = [result.properties['url']] + others = [x for x in result.properties['additionalUrls'] if 'tif' in x] + urls += others + + names = [Path(x).name for x in urls] + copol, crosspol, mask = [None, None, None] + + for name in names: + image_type = name.split('_')[-1].split('.')[0] + if image_type in ['VV', 'HH']: + copol = name + + if image_type in ['VH', 'HV']: + crosspol = name + + if image_type == 'mask': + mask = name + + if copol is None or crosspol is None: + raise ValueError('Both co-pol AND cross-pol data must be available to create an RGB decomposition') + + asf_search.download_urls(urls, path='.') + return copol, crosspol, mask + + +def normalize_browse_image_band( + band_image: np.ndarray, min_percentile: int = 3, max_percentile: int = 97 +) -> np.ndarray: + """Normalize a single band of a browse image to remove outliers. + + Args: + band_image: A single band numpy array to normalize. + min_percentile: Lower percentile to clip outliers to. + max_percentile: Upper percentile to clip outliers to. + + Returns: + A normalized numpy array. + """ + vmin = np.nanpercentile(band_image, min_percentile) + vmax = np.nanpercentile(band_image, max_percentile) + + # gamma correction: 0.5 + is_not_negative = band_image - vmin >= 0 + is_negative = band_image - vmin < 0 + band_image[is_not_negative] = np.sqrt((band_image[is_not_negative] - vmin) / (vmax - vmin)) + band_image[is_negative] = 0 + band_image = np.clip(band_image, 0, 1) + return band_image + + +def create_decomposition_rgb( + copol_path: str, crosspol_path: str, output_path: str, alpha_path: Optional[str] = None +) -> None: + """Create a georeferenced RGB decomposition RTC image from co-pol and cross-pol RTCs. + + Args: + copol_path: Path to co-pol RTC. + crosspol_path: Path to cross-pol RTC. + output_path: Path to save resulting RGB image to. + alpha_path: Path to alpha band image. If not provided, the alpha band will be the valid data mask. + """ + band_list = [None, None, None] + + for filename, is_copol in ((copol_path, True), (crosspol_path, False)): + gdal_ds = gdal.Open(filename, gdal.GA_ReadOnly) + + gdal_band = gdal_ds.GetRasterBand(1) + band_image = np.asarray(gdal_band.ReadAsArray(), dtype=np.float32) + + if is_copol: + # Using copol as template for output + is_valid = np.isfinite(band_image) + if alpha_path is None: + alpha_band = np.asarray(is_valid, dtype=np.float32) + else: + alpha_ds = gdal.Open(alpha_path, gdal.GA_ReadOnly) + alpha_band = np.asarray(alpha_ds.GetRasterBand(1).ReadAsArray(), dtype=np.float32) + + geotransform = gdal_ds.GetGeoTransform() + projection = gdal_ds.GetProjection() + shape = band_image.shape + band_list_index = [0, 2] + else: + band_list_index = [1] + + normalized_image = normalize_browse_image_band(band_image) + for index in band_list_index: + band_list[index] = normalized_image + + image = np.dstack((band_list[0], band_list[1], band_list[2], alpha_band)) + + driver = gdal.GetDriverByName('GTiff') + output_ds = driver.Create(output_path, shape[1], shape[0], 4, gdal.GDT_Float32) + output_ds.SetGeoTransform(geotransform) + output_ds.SetProjection(projection) + for i in range(4): + output_ds.GetRasterBand(i + 1).WriteArray(image[:, :, i]) + + +def create_gibs_formatted_image(input_file_name: str, output_file_name: str) -> str: + translated_file_name = 'translated.tif' + + warp_config = gdal.WarpOptions( + dstSRS='EPSG:4326', xRes=0.000274658203125, yRes=0.000274658203125, srcNodata=0, dstAlpha=True + ) + + translate_config = gdal.TranslateOptions( + resampleAlg='average', + format='COG', + creationOptions=['OVERVIEWS=NONE'], + outputType=gdalconst.GDT_Byte, + noData=0, + ) + + gdal.Translate(translated_file_name, input_file_name, options=translate_config) + gdal.Warp(output_file_name, translated_file_name, options=warp_config) + + return output_file_name + + +def create_gibs_rgb(granule, outpath): + copol, crosspol, mask = prep_data(granule) + unprojected_name = 'rgb_unprojected.tif' + create_decomposition_rgb(copol, crosspol, unprojected_name, mask) + create_gibs_formatted_image(unprojected_name, outpath) + + +def main(): + """Entrypoint for OPERA RTC decomposition image creation.""" + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) + parser.add_argument('granule', nargs=1, default=None, help='OPERA granule name') + parser.add_argument('--outpath', default='rgb.tif', help='Path to save resulting RGB image to') + args = parser.parse_args() + + create_gibs_rgb(args.granule[0], args.outpath)