From cd756d6371c2b03de912cea6c299ed11f9b654d8 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 17 Jan 2024 10:00:44 -0600 Subject: [PATCH 1/7] create first version of RGB workflow --- environment.yml | 1 + pyproject.toml | 2 + src/asf_tools/opera_rgb.py | 90 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 93 insertions(+) create mode 100644 src/asf_tools/opera_rgb.py 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..7ec910f1 --- /dev/null +++ b/src/asf_tools/opera_rgb.py @@ -0,0 +1,90 @@ +"""Create a georefernced RGB decomposition RTC image from two input OPERA RTCs""" +import argparse +from pathlib import Path + +import asf_search +import numpy as np +from osgeo import gdal + + +gdal.UseExceptions() + +BROWSE_IMAGE_MIN_PERCENTILE = 3 +BROWSE_IMAGE_MAX_PERCENTILE = 97 + + +def normalize_browse_image_band(band_image): + vmin = np.nanpercentile(band_image, BROWSE_IMAGE_MIN_PERCENTILE) + vmax = np.nanpercentile(band_image, BROWSE_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_browse_imagery(copol_path, crosspol_path, output_path, alpha_channel=None): + 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_channel is None: + alpha_channel = np.asarray(is_valid, 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_channel)) + + 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 prep_data(granule): + 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 + asf_search.download_urls(urls, path='.') + names = [Path(x).name for x in urls] + for name in names: + image_type = name.split('_')[-1].split('.')[0] + if image_type in ['VV', 'HH']: + copol = name + elif image_type in ['VH', 'HV']: + crosspol = name + elif image_type == 'mask': + mask = name + return copol, crosspol, mask + + +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() + + copol, crosspol, _ = prep_data(args.granule[0]) + create_browse_imagery(Path(copol), Path(crosspol), Path(args.outpath)) From eb54cce5d7d42f81977f8f97129cd6817e826e0a Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 17 Jan 2024 10:03:17 -0600 Subject: [PATCH 2/7] update changelog --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) 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] From 660cf08fbb50919cfcacf23ba46a1723b127a7a3 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 17 Jan 2024 12:43:05 -0600 Subject: [PATCH 3/7] add ability to grab mask data from OPERA RTC products --- src/asf_tools/opera_rgb.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/src/asf_tools/opera_rgb.py b/src/asf_tools/opera_rgb.py index 7ec910f1..4d390599 100644 --- a/src/asf_tools/opera_rgb.py +++ b/src/asf_tools/opera_rgb.py @@ -26,7 +26,7 @@ def normalize_browse_image_band(band_image): return band_image -def create_browse_imagery(copol_path, crosspol_path, output_path, alpha_channel=None): +def create_browse_imagery(copol_path, crosspol_path, output_path, alpha_path=None): band_list = [None, None, None] for filename, is_copol in ((copol_path, True), (crosspol_path, False)): @@ -38,8 +38,12 @@ def create_browse_imagery(copol_path, crosspol_path, output_path, alpha_channel= if is_copol: # Using copol as template for output is_valid = np.isfinite(band_image) - if alpha_channel is None: - alpha_channel = np.asarray(is_valid, dtype=np.float32) + 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 @@ -51,7 +55,7 @@ def create_browse_imagery(copol_path, crosspol_path, output_path, alpha_channel= for index in band_list_index: band_list[index] = normalized_image - image = np.dstack((band_list[0], band_list[1], band_list[2], alpha_channel)) + 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) @@ -66,16 +70,25 @@ def prep_data(granule): urls = [result.properties['url']] others = [x for x in result.properties['additionalUrls'] if 'tif' in x] urls += others - asf_search.download_urls(urls, path='.') + 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 - elif image_type in ['VH', 'HV']: + + if image_type in ['VH', 'HV']: crosspol = name - elif image_type == 'mask': + + 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 @@ -86,5 +99,5 @@ def main(): parser.add_argument('--outpath', default='rgb.tif', help='Path to save resulting RGB image to') args = parser.parse_args() - copol, crosspol, _ = prep_data(args.granule[0]) - create_browse_imagery(Path(copol), Path(crosspol), Path(args.outpath)) + copol, crosspol, mask = prep_data(args.granule[0]) + create_browse_imagery(copol, crosspol, args.outpath, mask) From 516fe03d201080fefd6c5311cb0af2c0f79e4bf9 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 17 Jan 2024 12:48:26 -0600 Subject: [PATCH 4/7] add docstrings and typing --- src/asf_tools/opera_rgb.py | 35 +++++++++++++++++++++++++++++++---- 1 file changed, 31 insertions(+), 4 deletions(-) diff --git a/src/asf_tools/opera_rgb.py b/src/asf_tools/opera_rgb.py index 4d390599..7e5ed0c3 100644 --- a/src/asf_tools/opera_rgb.py +++ b/src/asf_tools/opera_rgb.py @@ -1,6 +1,7 @@ """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 @@ -13,7 +14,15 @@ BROWSE_IMAGE_MAX_PERCENTILE = 97 -def normalize_browse_image_band(band_image): +def normalize_browse_image_band(band_image: np.ndarray) -> np.ndarray: + """Normalize a single band of a browse image to remove outliers. + + Args: + band_image: A single band numpy array to normalize. + + Returns: + A normalized numpy array. + """ vmin = np.nanpercentile(band_image, BROWSE_IMAGE_MIN_PERCENTILE) vmax = np.nanpercentile(band_image, BROWSE_IMAGE_MAX_PERCENTILE) @@ -26,7 +35,17 @@ def normalize_browse_image_band(band_image): return band_image -def create_browse_imagery(copol_path, crosspol_path, output_path, alpha_path=None): +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)): @@ -65,7 +84,15 @@ def create_browse_imagery(copol_path, crosspol_path, output_path, alpha_path=Non output_ds.GetRasterBand(i + 1).WriteArray(image[:, :, i]) -def prep_data(granule): +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] @@ -100,4 +127,4 @@ def main(): args = parser.parse_args() copol, crosspol, mask = prep_data(args.granule[0]) - create_browse_imagery(copol, crosspol, args.outpath, mask) + create_decomposition_rgb(copol, crosspol, args.outpath, mask) From c11d8db4c0e1c3102b94e7a96051d5bdc19f68f3 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 25 Jan 2024 13:34:31 -0600 Subject: [PATCH 5/7] remove un-needed variables --- src/asf_tools/opera_rgb.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/asf_tools/opera_rgb.py b/src/asf_tools/opera_rgb.py index 7e5ed0c3..88c485d1 100644 --- a/src/asf_tools/opera_rgb.py +++ b/src/asf_tools/opera_rgb.py @@ -10,9 +10,6 @@ gdal.UseExceptions() -BROWSE_IMAGE_MIN_PERCENTILE = 3 -BROWSE_IMAGE_MAX_PERCENTILE = 97 - def normalize_browse_image_band(band_image: np.ndarray) -> np.ndarray: """Normalize a single band of a browse image to remove outliers. @@ -23,8 +20,8 @@ def normalize_browse_image_band(band_image: np.ndarray) -> np.ndarray: Returns: A normalized numpy array. """ - vmin = np.nanpercentile(band_image, BROWSE_IMAGE_MIN_PERCENTILE) - vmax = np.nanpercentile(band_image, BROWSE_IMAGE_MAX_PERCENTILE) + vmin = np.nanpercentile(band_image, 3) + vmax = np.nanpercentile(band_image, 97) # gamma correction: 0.5 is_not_negative = band_image - vmin >= 0 From 11406222496ce12defeb5d8938557cc301d771f3 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 25 Jan 2024 14:01:37 -0600 Subject: [PATCH 6/7] make percentiles function args --- src/asf_tools/opera_rgb.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/asf_tools/opera_rgb.py b/src/asf_tools/opera_rgb.py index 88c485d1..2030104e 100644 --- a/src/asf_tools/opera_rgb.py +++ b/src/asf_tools/opera_rgb.py @@ -11,17 +11,21 @@ gdal.UseExceptions() -def normalize_browse_image_band(band_image: np.ndarray) -> np.ndarray: +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, 3) - vmax = np.nanpercentile(band_image, 97) + 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 From 9d02e90fb54b1a8a68ba4165c25b8cb428f5b810 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 25 Jan 2024 15:01:41 -0600 Subject: [PATCH 7/7] add code to make GIBS-compatible tifs --- src/asf_tools/opera_rgb.py | 87 +++++++++++++++++++++++++------------- 1 file changed, 57 insertions(+), 30 deletions(-) diff --git a/src/asf_tools/opera_rgb.py b/src/asf_tools/opera_rgb.py index 2030104e..aa3c1be2 100644 --- a/src/asf_tools/opera_rgb.py +++ b/src/asf_tools/opera_rgb.py @@ -5,12 +5,47 @@ import asf_search import numpy as np -from osgeo import gdal +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: @@ -85,39 +120,32 @@ def create_decomposition_rgb( output_ds.GetRasterBand(i + 1).WriteArray(image[:, :, i]) -def prep_data(granule: str): - """Download and prepare the data needed to create an RGB decomposition image. +def create_gibs_formatted_image(input_file_name: str, output_file_name: str) -> str: + translated_file_name = 'translated.tif' - Args: - granule: OPERA granule name. + warp_config = gdal.WarpOptions( + dstSRS='EPSG:4326', xRes=0.000274658203125, yRes=0.000274658203125, srcNodata=0, dstAlpha=True + ) - 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 + translate_config = gdal.TranslateOptions( + resampleAlg='average', + format='COG', + creationOptions=['OVERVIEWS=NONE'], + outputType=gdalconst.GDT_Byte, + noData=0, + ) - names = [Path(x).name for x in urls] - copol, crosspol, mask = [None, None, None] + gdal.Translate(translated_file_name, input_file_name, options=translate_config) + gdal.Warp(output_file_name, translated_file_name, options=warp_config) - for name in names: - image_type = name.split('_')[-1].split('.')[0] - if image_type in ['VV', 'HH']: - copol = name + return output_file_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 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(): @@ -127,5 +155,4 @@ def main(): parser.add_argument('--outpath', default='rgb.tif', help='Path to save resulting RGB image to') args = parser.parse_args() - copol, crosspol, mask = prep_data(args.granule[0]) - create_decomposition_rgb(copol, crosspol, args.outpath, mask) + create_gibs_rgb(args.granule[0], args.outpath)