From 12d2a8e3f6677cd0ccc5fe320b0acb74d8aaa3a2 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Thu, 13 Oct 2022 19:33:20 -0700 Subject: [PATCH 01/21] update Dockerfile to Python 3.9 --- docker/Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docker/Dockerfile b/docker/Dockerfile index 2fe6ebd..9736255 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -1,5 +1,5 @@ FROM ubuntu:20.04 -FROM python:3.8 +FROM python:3.9 # Set default UID and GID ENV USER_ID 1000 From 13e5b9daa34a94b19c0489b3479111e81702b86e Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Wed, 19 Oct 2022 16:54:59 -0700 Subject: [PATCH 02/21] Fix DIAG no data value (from 34464 to 65535) --- src/proteus/dswx_hls.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 7991a03..0a4d44a 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -76,7 +76,7 @@ 'qa': 'Fmask'} DIAGNOSTIC_LAYER_NO_DATA_DECIMAL = 0b100000 -DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR = 100000 +DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR = 65535 interpreted_dswx_band_dict = { 0b00000 : 0, # (Not Water) @@ -3003,7 +3003,12 @@ def _get_binary_representation(diagnostic_layer_decimal, nbits=6): for i in range(nbits): diagnostic_layer_decimal, bit_array = \ np.divmod(diagnostic_layer_decimal, 2) - diagnostic_layer_binary += bit_array * (10 ** i) + if i < 5: + diagnostic_layer_binary += bit_array * (10 ** i) + else: + # UInt16 max value is 65535 + diagnostic_layer_binary[bit_array] = \ + DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR return diagnostic_layer_binary From bc2a83c9ef9e4509393c1f218611d64db2cd1f2a Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Sat, 19 Nov 2022 10:27:07 -0800 Subject: [PATCH 03/21] add cast shadow masking --- bin/dswx_hls.py | 1 + src/proteus/defaults/dswx_hls.yaml | 3 + src/proteus/dswx_hls.py | 104 +++++++++++++++++++++++++---- src/proteus/schemas/dswx_hls.yaml | 3 + tests/test_dswx_hls_workflow.py | 1 + 5 files changed, 99 insertions(+), 13 deletions(-) diff --git a/bin/dswx_hls.py b/bin/dswx_hls.py index 7383165..146cae0 100755 --- a/bin/dswx_hls.py +++ b/bin/dswx_hls.py @@ -78,6 +78,7 @@ def main(): flag_use_otsu_terrain_masking=args.flag_use_otsu_terrain_masking, min_slope_angle = args.min_slope_angle, max_sun_local_inc_angle=args.max_sun_local_inc_angle, + apply_cast_shadow_masking=args.apply_cast_shadow_masking, mask_adjacent_to_cloud_mode=args.mask_adjacent_to_cloud_mode, flag_debug=args.flag_debug) diff --git a/src/proteus/defaults/dswx_hls.yaml b/src/proteus/defaults/dswx_hls.yaml index 08ab603..ddc20da 100644 --- a/src/proteus/defaults/dswx_hls.yaml +++ b/src/proteus/defaults/dswx_hls.yaml @@ -61,6 +61,9 @@ runconfig: # Maximum sun local-incidence angle in degrees for terrain masking max_sun_local_inc_angle: 40 + # Apply cast shadow masking + apply_cast_shadow_masking: False + # Define how areas adjacent to cloud/cloud-shadow should be handled mask_adjacent_to_cloud_mode: 'mask' diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 0a4d44a..ced4525 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -12,6 +12,7 @@ from osgeo.gdalconst import GDT_Float32 from osgeo import gdal, osr from scipy.ndimage import binary_dilation +import scipy from proteus.core import save_as_cog @@ -319,6 +320,8 @@ class RunConfigConstants: Minimum slope angle max_sun_local_inc_angle: float Maximum local-incidence angle + apply_cast_shadow_masking: bool + Apply cast shadow masking mask_adjacent_to_cloud_mode: str Define how areas adjacent to cloud/cloud-shadow should be handled. Options: "mask", "ignore", and "cover" @@ -332,6 +335,7 @@ def __init__(self): self.flag_use_otsu_terrain_masking = None self.min_slope_angle = None self.max_sun_local_inc_angle = None + self.apply_cast_shadow_masking = None self.mask_adjacent_to_cloud_mode = None self.browse_image_height = None self.browse_image_width = None @@ -531,6 +535,12 @@ def get_dswx_hls_cli_parser(): type=float, help='Maximum local-incidence angle') + parser.add_argument('--apply-cast-shadow-masking', + dest='apply_cast_shadow_masking', + action='store_true', + default=None, + help='Apply cast shadow masking') + parser.add_argument('--mask-adjacent-to-cloud-mode', dest='mask_adjacent_to_cloud_mode', type=str, @@ -2916,6 +2926,7 @@ def _compute_hillshade(dem_file, scratch_dir, sun_azimuth_angle, def _compute_opera_shadow_layer(dem, sun_azimuth_angle, sun_elevation_angle, min_slope_angle, max_sun_local_inc_angle, + apply_cast_shadow_masking, pixel_spacing_x = 30, pixel_spacing_y = 30): """Compute hillshade using new OPERA shadow masking @@ -2933,6 +2944,8 @@ def _compute_opera_shadow_layer(dem, sun_azimuth_angle, sun_elevation_angle, Minimum slope angle max_sun_local_inc_angle: float Maximum local-incidence angle + apply_cast_shadow_masking: bool + Apply cast shadow masking pixel_spacing_x: float (optional) Pixel spacing in the X direction pixel_spacing_y: float (optional) @@ -2947,11 +2960,17 @@ def _compute_opera_shadow_layer(dem, sun_azimuth_angle, sun_elevation_angle, sun_zenith_degrees = 90 - sun_elevation_angle sun_zenith = np.radians(sun_zenith_degrees) + # vector coordinates: (x, y, z) target_to_sun_unit_vector = [np.sin(sun_azimuth) * np.sin(sun_zenith), np.cos(sun_azimuth) * np.sin(sun_zenith), np.cos(sun_zenith)] + # numpy computes the gradient as (y, x) gradient_h = np.gradient(dem) + + # The terrain normal vector is calculated as: + # N = [-dh/dx, -dh/dy, 1] where dh/dx is the west-east slope + # and dh/dy is the south-north slope wrt to the DEM grid terrain_normal_vector = [-gradient_h[1] / pixel_spacing_x, -gradient_h[0] / - abs(pixel_spacing_y), 1] @@ -2975,6 +2994,52 @@ def _compute_opera_shadow_layer(dem, sun_azimuth_angle, sun_elevation_angle, low_sun_inc_angle_mask = sun_inc_angle_degrees <= max_sun_local_inc_angle shadow_mask = (low_sun_inc_angle_mask | (~ backslope_mask)) + + if apply_cast_shadow_masking: + + # approx_sun_distance_to_earth_km = 1.50e8 + # 1.5e8 meters rather than 1.5e8 km + fake_sun_distance_to_earth = 1.50e6 + target_to_sun_unit_vector = np.asarray(target_to_sun_unit_vector, dtype=np.double) + sun_position = fake_sun_distance_to_earth * target_to_sun_unit_vector + length = dem.shape[0] + width = dem.shape[1] + y_dist_vect = np.linspace(- length // 2, length // 2 - 1, length) * -abs(pixel_spacing_y) + x_dist_vect = np.linspace(- width // 2, width // 2 - 1, width) * pixel_spacing_x + x_dist, y_dist = np.meshgrid(x_dist_vect, y_dist_vect) + distance_from_sun = np.sqrt((sun_position[0] - x_dist) ** 2 + + (sun_position[1] - y_dist) ** 2 + + (sun_position[2] - dem) ** 2) + + sun_inc_angle = np.degrees(np.arccos(np.absolute(sun_position[2] - dem) / distance_from_sun)) + + n_shadow_pixels = None + cast_shadow = None + shift_factor = 0 + + # max shift_factor is 500m + CAST_SHADOW_LONGEST_DIST = 500 + + while (n_shadow_pixels is None or (n_shadow_pixels > 0 and shift_factor < CAST_SHADOW_LONGEST_DIST)): + shift_factor += 5 + + positive_shift = [shift_factor * np.cos(sun_azimuth), + - shift_factor * np.sin(sun_azimuth)] + + sun_inc_angle_positive_shift = scipy.ndimage.shift(sun_inc_angle, positive_shift) + + # cast_shadow values: 1: shadow 0: non-zero + new_cast_shadow = sun_inc_angle < sun_inc_angle_positive_shift + + if cast_shadow is None: + cast_shadow = new_cast_shadow + else: + cast_shadow |= new_cast_shadow + + n_shadow_pixels = np.count_nonzero(new_cast_shadow) + + shadow_mask = (shadow_mask & (~cast_shadow)) + return shadow_mask @@ -3064,6 +3129,7 @@ def generate_dswx_layers(input_list, flag_use_otsu_terrain_masking=None, min_slope_angle=None, max_sun_local_inc_angle=None, + apply_cast_shadow_masking=None, mask_adjacent_to_cloud_mode=None, flag_debug=False): """Compute the DSWx-HLS product @@ -3135,6 +3201,8 @@ def generate_dswx_layers(input_list, Minimum slope angle max_sun_local_inc_angle: float (optional) Maximum local-incidence angle + apply_cast_shadow_masking: bool + Apply cast shadow masking flag_debug: bool (optional) Flag to indicate if execution is for debug purposes. If so, only a subset of the image will be loaded into memory @@ -3148,13 +3216,15 @@ def generate_dswx_layers(input_list, Flag success indicating if execution was successful """ - flag_read_runconfig_constants = (hls_thresholds is None or - flag_use_otsu_terrain_masking is None or - min_slope_angle is None or - max_sun_local_inc_angle is None or - mask_adjacent_to_cloud_mode is None or - browse_image_height is None or - browse_image_width is None) + flag_read_runconfig_constants = \ + any([p is None for p in [hls_thresholds, + flag_use_otsu_terrain_masking, + min_slope_angle, + max_sun_local_inc_angle, + apply_cast_shadow_masking, + mask_adjacent_to_cloud_mode, + browse_image_height, + browse_image_width]]) if flag_read_runconfig_constants: runconfig_constants = parse_runconfig_file() @@ -3166,6 +3236,8 @@ def generate_dswx_layers(input_list, min_slope_angle = runconfig_constants.min_slope_angle if max_sun_local_inc_angle is None: max_sun_local_inc_angle = runconfig_constants.max_sun_local_inc_angle + if apply_cast_shadow_masking is None: + apply_cast_shadow_masking = runconfig_constants.apply_cast_shadow_masking if mask_adjacent_to_cloud_mode is None: mask_adjacent_to_cloud_mode = runconfig_constants.mask_adjacent_to_cloud_mode if browse_image_height is None: @@ -3194,10 +3266,16 @@ def generate_dswx_layers(input_list, logger.info(f' product version: {product_version}') logger.info(f' software version: {SOFTWARE_VERSION}') logger.info(f'processing parameters:') - logger.info(f' flag_use_otsu_terrain_masking: {flag_use_otsu_terrain_masking}') - logger.info(f' min_slope_angle: {min_slope_angle}') - logger.info(f' max_sun_local_inc_angle: {max_sun_local_inc_angle}') - logger.info(f' mask_adjacent_to_cloud_mode: {mask_adjacent_to_cloud_mode}') + if flag_use_otsu_terrain_masking: + terrain_masking_algorighm = 'Otsu' + else: + terrain_masking_algorighm = 'sun local incidence angle' + if not flag_use_otsu_terrain_masking: + logger.info(f' terrain masking algorithm: {terrain_masking_algorighm}') + logger.info(f' min. slope angle: {min_slope_angle}') + logger.info(f' max. sun local inc. angle: {max_sun_local_inc_angle}') + logger.info(f' apply cast shadow masking: {apply_cast_shadow_masking}') + logger.info(f' mask adjacent cloud/cloud-shadow mode: {mask_adjacent_to_cloud_mode}') if output_browse_image: logger.info(f'browse image:') logger.info(f' browse_image_height: {browse_image_height}') @@ -3308,8 +3386,8 @@ def generate_dswx_layers(input_list, # new OPERA shadow masking shadow_layer_with_margin = _compute_opera_shadow_layer( dem_with_margin, sun_azimuth_angle, sun_elevation_angle, - min_slope_angle = min_slope_angle, - max_sun_local_inc_angle = max_sun_local_inc_angle) + min_slope_angle, max_sun_local_inc_angle, + apply_cast_shadow_masking) # remove extra margin from shadow_layer shadow_layer = _crop_2d_array_all_sides(shadow_layer_with_margin, diff --git a/src/proteus/schemas/dswx_hls.yaml b/src/proteus/schemas/dswx_hls.yaml index 094eba6..f8015e8 100644 --- a/src/proteus/schemas/dswx_hls.yaml +++ b/src/proteus/schemas/dswx_hls.yaml @@ -62,6 +62,9 @@ runconfig: # Maximum sun local-incidence angle in degrees for terrain masking max_sun_local_inc_angle: num(min=-180, max=180, required=False) + # Apply cast shadow masking + apply_cast_shadow_masking: bool(required=False) + # Define how areas adjacent to cloud/cloud-shadow should be handled mask_adjacent_to_cloud_mode: enum('mask', 'ignore', 'cover', required=False) diff --git a/tests/test_dswx_hls_workflow.py b/tests/test_dswx_hls_workflow.py index a455534..dcacbe8 100644 --- a/tests/test_dswx_hls_workflow.py +++ b/tests/test_dswx_hls_workflow.py @@ -90,6 +90,7 @@ def test_workflow(): flag_use_otsu_terrain_masking=args.flag_use_otsu_terrain_masking, min_slope_angle = args.min_slope_angle, max_sun_local_inc_angle=args.max_sun_local_inc_angle, + apply_cast_shadow_masking=args.apply_cast_shadow_masking, mask_adjacent_to_cloud_mode=args.mask_adjacent_to_cloud_mode, flag_debug=args.flag_debug) From 4a2e9e7d242343e9925c303926f056d77e51b98a Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Mon, 5 Dec 2022 08:30:47 -0800 Subject: [PATCH 04/21] use tempfile.NamedTemporaryFile() for temporary variables --- src/proteus/dswx_hls.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index ced4525..908a32d 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -814,9 +814,10 @@ def create_landcover_mask(copernicus_landcover_file, logger.info(f'copernicus landcover 100 m file: {copernicus_landcover_file}') logger.info(f'worldcover 10 m file: {worldcover_file}') - # Reproject Copernicus land cover - copernicus_landcover_reprojected_file = os.path.join( - scratch_dir, 'copernicus_reprojected.tif') + # Reproject Copernicus land cover + copernicus_landcover_reprojected_file = tempfile.NamedTemporaryFile( + suffix='copernicus_reprojected', dir=scratch_dir, suffix='.tif').name + copernicus_landcover_array = _warp(copernicus_landcover_file, geotransform, projection, length, width, scratch_dir, resample_algorithm='nearest', @@ -830,8 +831,9 @@ def create_landcover_mask(copernicus_landcover_file, geotransform_up_3[5] = geotransform[5] / 3 # dy / 3 length_up_3 = 3 * length width_up_3 = 3 * width - worldcover_reprojected_up_3_file = os.path.join( - scratch_dir, 'worldcover_reprojected_up_3.tif') + worldcover_reprojected_up_3_file = tempfile.NamedTemporaryFile( + suffix='worldcover_reprojected_up_3', + dir=scratch_dir, suffix='.tif').name worldcover_array_up_3 = _warp(worldcover_file, geotransform_up_3, projection, length_up_3, width_up_3, scratch_dir, resample_algorithm='nearest', From 1405f9309b7b52727828154185c9f7f4212687d9 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Wed, 7 Dec 2022 18:42:37 -0800 Subject: [PATCH 05/21] update the naming of uncollapsed water classes 3 and 4 --- src/proteus/dswx_hls.py | 42 +++++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 18 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 908a32d..12ad13c 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -124,8 +124,8 @@ WTR_COLLAPSED_PARTIAL_SURFACE_WATER = 2 WTR_UNCOLLAPSED_HIGH_CONF_WATER = 1 WTR_UNCOLLAPSED_MODERATE_CONF_WATER = 2 -WTR_UNCOLLAPSED_POTENTIAL_WETLAND = 3 -WTR_UNCOLLAPSED_LOW_CONF_WATER = 4 +WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_CONSERVATIVE = 3 +WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_AGGRESSIVE = 4 FIRST_UNCOLLAPSED_WATER_CLASS = 1 LAST_UNCOLLAPSED_WATER_CLASS = 4 @@ -150,8 +150,8 @@ 1. High-confidence water; 2. Moderate-confidence water; -3. Potential wetland; -4. Low-confidence water or wetland. +3. Partial surface water conservative (previous potential wetland); +4. Partial surface water aggressive (previous low-confidence water or wetland). These classes are collapsed into 2 classes when DSWx-HLS WTR layers are saved: @@ -162,8 +162,10 @@ WTR_NOT_WATER: WTR_NOT_WATER, WTR_UNCOLLAPSED_HIGH_CONF_WATER: WTR_COLLAPSED_OPEN_WATER, WTR_UNCOLLAPSED_MODERATE_CONF_WATER: WTR_COLLAPSED_OPEN_WATER, - WTR_UNCOLLAPSED_POTENTIAL_WETLAND: WTR_COLLAPSED_PARTIAL_SURFACE_WATER, - WTR_UNCOLLAPSED_LOW_CONF_WATER: WTR_COLLAPSED_PARTIAL_SURFACE_WATER, + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_CONSERVATIVE: \ + WTR_COLLAPSED_PARTIAL_SURFACE_WATER, + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_AGGRESSIVE: \ + WTR_COLLAPSED_PARTIAL_SURFACE_WATER, WTR_CLOUD_MASKED_SNOW: WTR_CLOUD_MASKED_SNOW, WTR_CLOUD_MASKED: WTR_CLOUD_MASKED, UINT8_FILL_VALUE: UINT8_FILL_VALUE @@ -190,8 +192,8 @@ WTR_NOT_WATER: CONF_NOT_WATER, WTR_UNCOLLAPSED_HIGH_CONF_WATER: 95, WTR_UNCOLLAPSED_MODERATE_CONF_WATER: 80, - WTR_UNCOLLAPSED_POTENTIAL_WETLAND: 70, - WTR_UNCOLLAPSED_LOW_CONF_WATER: 60, + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_CONSERVATIVE: 70, + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_AGGRESSIVE: 60, WTR_CLOUD_MASKED_SNOW: CONF_CLOUD_MASKED_SNOW, WTR_CLOUD_MASKED: CONF_CLOUD_MASKED, UINT8_FILL_VALUE: UINT8_FILL_VALUE @@ -1056,16 +1058,20 @@ def _apply_landcover_and_shadow_masks(interpreted_layer, nir, to_mask_ind = np.where( _is_landcover_class_evergreen(landcover_mask) & (nir > hls_thresholds.lcmask_nir) & - ((interpreted_layer == WTR_UNCOLLAPSED_POTENTIAL_WETLAND) | - (interpreted_layer == WTR_UNCOLLAPSED_LOW_CONF_WATER))) + ((interpreted_layer == + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_CONSERVATIVE) | + (interpreted_layer == + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_AGGRESSIVE))) landcover_shadow_masked_dswx[to_mask_ind] = WTR_NOT_WATER # Check landcover (low intensity developed) to_mask_ind = np.where( _is_landcover_class_low_intensity_developed(landcover_mask) & (nir > hls_thresholds.lcmask_nir) & - ((interpreted_layer == WTR_UNCOLLAPSED_POTENTIAL_WETLAND) | - (interpreted_layer == WTR_UNCOLLAPSED_LOW_CONF_WATER))) + ((interpreted_layer == + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_CONSERVATIVE) | + (interpreted_layer == + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_AGGRESSIVE))) landcover_shadow_masked_dswx[to_mask_ind] = WTR_NOT_WATER # Check landcover (high intensity developed) @@ -1116,12 +1122,12 @@ def _get_interpreted_dswx_ctable( # Light blue - Water (moderate conf.) dswx_ctable.SetColorEntry(WTR_UNCOLLAPSED_MODERATE_CONF_WATER, (0, 127, 255)) - # Dark green - Potential wetland - dswx_ctable.SetColorEntry(WTR_UNCOLLAPSED_POTENTIAL_WETLAND, - (0, 127, 0)) - # Green - Low confidence water or wetland - dswx_ctable.SetColorEntry(WTR_UNCOLLAPSED_LOW_CONF_WATER, - (0, 255, 0)) + # Dark green - Partial surface water conservative + dswx_ctable.SetColorEntry( + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_CONSERVATIVE, (0, 127, 0)) + # Green - Partial surface water aggressive + dswx_ctable.SetColorEntry( + WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_AGGRESSIVE, (0, 255, 0)) # Gray - QA masked (Cloud/cloud-shadow) dswx_ctable.SetColorEntry(WTR_CLOUD_MASKED, (127, 127, 127)) From 5e35d6452b9de9298d59ac00e87845537f6688d2 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Wed, 7 Dec 2022 18:43:01 -0800 Subject: [PATCH 06/21] simplify code --- src/proteus/core.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/proteus/core.py b/src/proteus/core.py index fac5556..113718c 100644 --- a/src/proteus/core.py +++ b/src/proteus/core.py @@ -28,7 +28,7 @@ def save_as_cog(filename, scratch_dir = '.', logger = None, logger = logging.getLogger('proteus') logger.info('COG step 1: add overviews') - gdal_ds = gdal.Open(filename, 1) + gdal_ds = gdal.Open(filename, gdal.GA_Update) gdal_dtype = gdal_ds.GetRasterBand(1).DataType dtype_name = gdal.GetDataTypeName(gdal_dtype).lower() @@ -62,7 +62,6 @@ def save_as_cog(filename, scratch_dir = '.', logger = None, if flag_compress: gdal_translate_options += ['COMPRESS=DEFLATE'] - is_integer = 'byte' in dtype_name or 'int' in dtype_name if is_integer: gdal_translate_options += ['PREDICTOR=2'] else: From 47b084eaa5bea4c8b8e6a714c94b9c77a408fe6d Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Wed, 7 Dec 2022 18:55:01 -0800 Subject: [PATCH 07/21] add placeholder for runconfig parameter `check_ancillary_inputs_coverage` --- src/proteus/defaults/dswx_hls.yaml | 4 ++++ src/proteus/dswx_hls.py | 1 - src/proteus/schemas/dswx_hls.yaml | 4 ++++ 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/proteus/defaults/dswx_hls.yaml b/src/proteus/defaults/dswx_hls.yaml index ddc20da..81c58fb 100644 --- a/src/proteus/defaults/dswx_hls.yaml +++ b/src/proteus/defaults/dswx_hls.yaml @@ -52,6 +52,10 @@ runconfig: product_version: processing: + + # Check if ancillary input cover entirely output products + check_ancillary_inputs_coverage: True + # Use terrain-masking based on the Otsu method flag_use_otsu_terrain_masking: False diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 12ad13c..8d4d815 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -119,7 +119,6 @@ WTR_NOT_WATER = 0 # Water classes - WTR_COLLAPSED_OPEN_WATER = 1 WTR_COLLAPSED_PARTIAL_SURFACE_WATER = 2 WTR_UNCOLLAPSED_HIGH_CONF_WATER = 1 diff --git a/src/proteus/schemas/dswx_hls.yaml b/src/proteus/schemas/dswx_hls.yaml index f8015e8..6c7e113 100644 --- a/src/proteus/schemas/dswx_hls.yaml +++ b/src/proteus/schemas/dswx_hls.yaml @@ -53,6 +53,10 @@ runconfig: product_version: num(required=False) processing: + + # Check if ancillary input cover entirely output products + check_ancillary_inputs_coverage: bool(required=False) + # Use terrain-masking based on the Otsu method flag_use_otsu_terrain_masking: bool(required=False) From 948d66ce5fb071e5677ad04267e4d837bf3ac5e4 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Sun, 11 Dec 2022 10:19:25 -0800 Subject: [PATCH 08/21] add 4.9-km HLS margin to MGRS 100 km x 100 km tiles --- src/proteus/core.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/proteus/core.py b/src/proteus/core.py index 113718c..1947ba1 100644 --- a/src/proteus/core.py +++ b/src/proteus/core.py @@ -89,7 +89,7 @@ def save_as_cog(filename, scratch_dir = '.', logger = None, f' optimized GeoTIFF!') -def get_geographic_boundaries_from_mgrs_tile(mgrs_tile_name, verbose=False): +def get_hls_geographic_boundaries_from_mgrs_tile(mgrs_tile_name, verbose=False): import mgrs mgrs_obj = mgrs.MGRS() @@ -121,8 +121,10 @@ def get_geographic_boundaries_from_mgrs_tile(mgrs_tile_name, verbose=False): for offset_x_multiplier in range(2): for offset_y_multiplier in range(2): - x = x_min + offset_x_multiplier * 109.8 * 1000 - y = y_min + offset_y_multiplier * 109.8 * 1000 + # We are using MGRS 100km x 100km tiles + # HLS tiles have 4.9 km of margin => width/length = 109.8 km + x = x_min - 4.9 * 1000 + offset_x_multiplier * 109.8 * 1000 + y = y_min - 4.9 * 1000 + offset_y_multiplier * 109.8 * 1000 lat, lon, z = transformation.TransformPoint(x, y, elevation) if verbose: From 89c36b272e51ece14af9596110ffa3a8fe3ff9c8 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Fri, 16 Dec 2022 09:09:26 -0800 Subject: [PATCH 09/21] fix `check_ancillary_inputs_coverage` comments --- src/proteus/defaults/dswx_hls.yaml | 2 +- src/proteus/dswx_hls.py | 3 ++- src/proteus/schemas/dswx_hls.yaml | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/proteus/defaults/dswx_hls.yaml b/src/proteus/defaults/dswx_hls.yaml index 81c58fb..b9766e1 100644 --- a/src/proteus/defaults/dswx_hls.yaml +++ b/src/proteus/defaults/dswx_hls.yaml @@ -53,7 +53,7 @@ runconfig: processing: - # Check if ancillary input cover entirely output products + # Check if ancillary inputs cover entirely the output product check_ancillary_inputs_coverage: True # Use terrain-masking based on the Otsu method diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 8d4d815..7a63803 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -616,7 +616,8 @@ def compare_dswx_hls_products(file_1, file_2): image_1 = gdal_band_1.ReadAsArray() image_2 = gdal_band_2.ReadAsArray() flag_bands_are_equal = np.allclose( - image_1, image_2, atol=COMPARE_DSWX_HLS_PRODUCTS_ERROR_TOLERANCE) + image_1, image_2, atol=COMPARE_DSWX_HLS_PRODUCTS_ERROR_TOLERANCE, + equal_nan=True) flag_bands_are_equal_str = _get_prefix_str(flag_bands_are_equal, flag_all_ok) print(f'{flag_bands_are_equal_str} Band {b} -' diff --git a/src/proteus/schemas/dswx_hls.yaml b/src/proteus/schemas/dswx_hls.yaml index 6c7e113..8f6c049 100644 --- a/src/proteus/schemas/dswx_hls.yaml +++ b/src/proteus/schemas/dswx_hls.yaml @@ -54,7 +54,7 @@ runconfig: processing: - # Check if ancillary input cover entirely output products + # Check if ancillary inputs cover entirely the output product check_ancillary_inputs_coverage: bool(required=False) # Use terrain-masking based on the Otsu method From 2f379eb95c539d3773b805fb7386b6c08ae46880 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Fri, 16 Dec 2022 09:11:44 -0800 Subject: [PATCH 10/21] use tempfile.NamedTemporaryFile() for temporary variables (2) --- src/proteus/dswx_hls.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 7a63803..31331e3 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -818,7 +818,7 @@ def create_landcover_mask(copernicus_landcover_file, # Reproject Copernicus land cover copernicus_landcover_reprojected_file = tempfile.NamedTemporaryFile( - suffix='copernicus_reprojected', dir=scratch_dir, suffix='.tif').name + dir=scratch_dir, suffix='.tif').name copernicus_landcover_array = _warp(copernicus_landcover_file, geotransform, projection, @@ -834,7 +834,6 @@ def create_landcover_mask(copernicus_landcover_file, length_up_3 = 3 * length width_up_3 = 3 * width worldcover_reprojected_up_3_file = tempfile.NamedTemporaryFile( - suffix='worldcover_reprojected_up_3', dir=scratch_dir, suffix='.tif').name worldcover_array_up_3 = _warp(worldcover_file, geotransform_up_3, projection, length_up_3, width_up_3, scratch_dir, From 15d4c2fd01eb18a8d5b8109cb6a1ac62014d9073 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Fri, 16 Dec 2022 09:23:23 -0800 Subject: [PATCH 11/21] rename HLS QA band to FMask and DSWx-HLS QA to CLOUD mask --- src/proteus/dswx_hls.py | 56 ++++++++++++++++++++--------------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 31331e3..86a7c03 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -50,7 +50,7 @@ 'nir': 'band05', 'swir1': 'band06', 'swir2': 'band07', - 'qa': 'QA'} + 'fmask': 'QA'} s30_v1_band_dict = {'blue': 'band02', 'green': 'band03', @@ -58,7 +58,7 @@ 'nir': 'band8A', 'swir1': 'band11', 'swir2': 'band12', - 'qa': 'QA'} + 'fmask': 'QA'} l30_v2_band_dict = {'blue': 'B02', 'green': 'B03', @@ -66,7 +66,7 @@ 'nir': 'B05', 'swir1': 'B06', 'swir2': 'B07', - 'qa': 'Fmask'} + 'fmask': 'Fmask'} s30_v2_band_dict = {'blue': 'B02', 'green': 'B03', @@ -74,7 +74,7 @@ 'nir': 'B8A', 'swir1': 'B11', 'swir2': 'B12', - 'qa': 'Fmask'} + 'fmask': 'Fmask'} DIAGNOSTIC_LAYER_NO_DATA_DECIMAL = 0b100000 DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR = 65535 @@ -1128,10 +1128,10 @@ def _get_interpreted_dswx_ctable( dswx_ctable.SetColorEntry( WTR_UNCOLLAPSED_PARTIAL_SURFACE_WATER_AGGRESSIVE, (0, 255, 0)) - # Gray - QA masked (Cloud/cloud-shadow) + # Gray - CLOUD masked (Cloud/cloud-shadow) dswx_ctable.SetColorEntry(WTR_CLOUD_MASKED, (127, 127, 127)) - # Cyan - QA masked (Snow) + # Cyan - CLOUD masked (Snow) dswx_ctable.SetColorEntry(WTR_CLOUD_MASKED_SNOW, (0, 255, 255)) # Black - Fill value @@ -1317,11 +1317,11 @@ def _get_binary_water_layer(interpreted_water_layer): for class_value in range(1, 5): binary_water_layer[interpreted_water_layer == class_value] = BWTR_WATER - # Q/A masked (cloud/snow) + # CLOUD masked (cloud/snow) binary_water_layer[interpreted_water_layer == WTR_CLOUD_MASKED_SNOW] = \ WTR_CLOUD_MASKED_SNOW - # Q/A masked (cloud/cloud-shadow) + # CLOUD masked (cloud/cloud-shadow) binary_water_layer[interpreted_water_layer == WTR_CLOUD_MASKED] = \ WTR_CLOUD_MASKED @@ -1437,7 +1437,7 @@ def _compute_diagnostic_tests(blue, green, red, nir, swir1, swir2, def _compute_mask_and_filter_interpreted_layer( - unmasked_interpreted_water_layer, qa_band, + unmasked_interpreted_water_layer, fmask, mask_adjacent_to_cloud_mode): """Compute cloud/cloud-shadow mask and filter interpreted water layer @@ -1445,8 +1445,8 @@ def _compute_mask_and_filter_interpreted_layer( ---------- unmasked_interpreted_water_layer: numpy.ndarray Cloud-unmasked interpreted water layer - qa_band: numpy ndarray - HLS Q/A band + fmask: numpy ndarray + HLS FMask mask_adjacent_to_cloud_mode: str Define how areas adjacent to cloud/cloud-shadow should be handled. Options: "mask", "ignore", and "cover" @@ -1461,7 +1461,7 @@ def _compute_mask_and_filter_interpreted_layer( mask = np.zeros(shape, dtype = np.uint8) ''' - QA band - Landsat 8 + HLS FMask BITS: 0 - Cirrus (reserved but not used) 1 - Cloud (*1) @@ -1489,26 +1489,26 @@ def _compute_mask_and_filter_interpreted_layer( logger.info(f'mask adjacent to cloud/cloud-shadow mode:' f' {mask_adjacent_to_cloud_mode}') - # Check QA cloud shadow bit (3) => bit 0 - mask[np.bitwise_and(qa_band, 2**3) == 2**3] = 1 + # Check FMask cloud shadow bit (3) => bit 0 + mask[np.bitwise_and(fmask, 2**3) == 2**3] = 1 if mask_adjacent_to_cloud_mode == 'mask': - # Check QA adjacent to cloud/shadow bit (2) => bit 0 - mask[np.bitwise_and(qa_band, 2**2) == 2**2] = 1 + # Check FMask adjacent to cloud/shadow bit (2) => bit 0 + mask[np.bitwise_and(fmask, 2**2) == 2**2] = 1 - # Check QA cloud bit (1) => bit 2 - mask[np.bitwise_and(qa_band, 2**1) == 2**1] += 4 + # Check FMask cloud bit (1) => bit 2 + mask[np.bitwise_and(fmask, 2**1) == 2**1] += 4 # If cloud (1) or cloud shadow (3), mark WTR as WTR_CLOUD_MASKED masked_interpreted_water_layer[mask != 0] = WTR_CLOUD_MASKED - # Check QA snow bit (4) => bit 1 - snow_mask = np.bitwise_and(qa_band, 2**4) == 2**4 + # Check FMask snow bit (4) => bit 1 + snow_mask = np.bitwise_and(fmask, 2**4) == 2**4 # Cover areas marked as adjacent to cloud/shadow if mask_adjacent_to_cloud_mode == 'cover': # Dilate snow mask over areas adjacent to cloud/shadow - adjacent_to_cloud_mask = np.bitwise_and(qa_band, 2**2) == 2**2 + adjacent_to_cloud_mask = np.bitwise_and(fmask, 2**2) == 2**2 areas_to_dilate = (adjacent_to_cloud_mask) & (mask == 0) snow_mask = binary_dilation(snow_mask, iterations=10, @@ -1586,7 +1586,7 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, hls_dataset_name = hls_dataset_name.replace(f'.{band_suffix}', '') image_dict['hls_dataset_name'] = hls_dataset_name - if key == 'qa': + if key == 'fmask': if flag_debug: logger.info('reading in debug mode') image_dict[key] = layer_gdal_dataset.ReadAsArray( @@ -1851,9 +1851,9 @@ def _get_binary_water_ctable(): binary_water_ctable.SetColorEntry(WTR_NOT_WATER, (255, 255, 255)) # Blue - Water binary_water_ctable.SetColorEntry(BWTR_WATER, (0, 0, 255)) - # Cyan - QA masked (snow) + # Cyan - CLOUD masked (snow) binary_water_ctable.SetColorEntry(WTR_CLOUD_MASKED_SNOW, (0, 255, 255)) - # Gray - QA masked (cloud/cloud-shadow) + # Gray - CLOUD masked (cloud/cloud-shadow) binary_water_ctable.SetColorEntry(WTR_CLOUD_MASKED, (127, 127, 127)) # Black (transparent) - Fill value binary_water_ctable.SetColorEntry(UINT8_FILL_VALUE, (0, 0, 0, 255)) @@ -1883,11 +1883,11 @@ def _get_confidence_layer_ctable(): # White - Not water confidence_layer_ctable.SetColorEntry(CONF_NOT_WATER, (255, 255, 255)) - # Cyan - QA masked (snow) + # Cyan - CLOUD masked (snow) confidence_layer_ctable.SetColorEntry(CONF_CLOUD_MASKED_SNOW, (0, 255, 255)) - # Gray - QA masked (cloud/cloud-shadow) + # Gray - CLOUD masked (cloud/cloud-shadow) confidence_layer_ctable.SetColorEntry(CONF_CLOUD_MASKED, (127, 127, 127)) # Black - Fill value @@ -3342,7 +3342,7 @@ def generate_dswx_layers(input_list, nir = image_dict['nir'] swir1 = image_dict['swir1'] swir2 = image_dict['swir2'] - qa = image_dict['qa'] + fmask = image_dict['fmask'] geotransform = image_dict['geotransform'] projection = image_dict['projection'] @@ -3510,7 +3510,7 @@ def generate_dswx_layers(input_list, output_files_list=build_vrt_list) cloud, masked_dswx_band = _compute_mask_and_filter_interpreted_layer( - landcover_shadow_masked_dswx, qa, + landcover_shadow_masked_dswx, fmask, mask_adjacent_to_cloud_mode) if invalid_ind is not None: From 6e6e56faa85ee861d3152aaa42c46d3c083f9f7e Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Fri, 16 Dec 2022 09:36:32 -0800 Subject: [PATCH 12/21] rename `fill_data` to `fill_value`; mask out points where Fmask == fill value --- src/proteus/dswx_hls.py | 43 ++++++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 18 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 86a7c03..e83d7a3 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -1446,7 +1446,7 @@ def _compute_mask_and_filter_interpreted_layer( unmasked_interpreted_water_layer: numpy.ndarray Cloud-unmasked interpreted water layer fmask: numpy ndarray - HLS FMask + HLS Fmask mask_adjacent_to_cloud_mode: str Define how areas adjacent to cloud/cloud-shadow should be handled. Options: "mask", "ignore", and "cover" @@ -1461,7 +1461,7 @@ def _compute_mask_and_filter_interpreted_layer( mask = np.zeros(shape, dtype = np.uint8) ''' - HLS FMask + HLS Fmask BITS: 0 - Cirrus (reserved but not used) 1 - Cloud (*1) @@ -1489,20 +1489,20 @@ def _compute_mask_and_filter_interpreted_layer( logger.info(f'mask adjacent to cloud/cloud-shadow mode:' f' {mask_adjacent_to_cloud_mode}') - # Check FMask cloud shadow bit (3) => bit 0 + # Check Fmask cloud shadow bit (3) => bit 0 mask[np.bitwise_and(fmask, 2**3) == 2**3] = 1 if mask_adjacent_to_cloud_mode == 'mask': - # Check FMask adjacent to cloud/shadow bit (2) => bit 0 + # Check Fmask adjacent to cloud/shadow bit (2) => bit 0 mask[np.bitwise_and(fmask, 2**2) == 2**2] = 1 - # Check FMask cloud bit (1) => bit 2 + # Check Fmask cloud bit (1) => bit 2 mask[np.bitwise_and(fmask, 2**1) == 2**1] += 4 # If cloud (1) or cloud shadow (3), mark WTR as WTR_CLOUD_MASKED masked_interpreted_water_layer[mask != 0] = WTR_CLOUD_MASKED - # Check FMask snow bit (4) => bit 1 + # Check Fmask snow bit (4) => bit 1 snow_mask = np.bitwise_and(fmask, 2**4) == 2**4 # Cover areas marked as adjacent to cloud/shadow @@ -1578,7 +1578,7 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, if layer_gdal_dataset is None: return None band = layer_gdal_dataset.GetRasterBand(1) - fill_data = band.GetNoDataValue() + fill_value = band.GetNoDataValue() if 'hls_dataset_name' not in image_dict.keys(): hls_dataset_name = os.path.splitext(os.path.basename(filename))[0] @@ -1586,6 +1586,8 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, hls_dataset_name = hls_dataset_name.replace(f'.{band_suffix}', '') image_dict['hls_dataset_name'] = hls_dataset_name + metadata = layer_gdal_dataset.GetMetadata() + if key == 'fmask': if flag_debug: logger.info('reading in debug mode') @@ -1593,13 +1595,12 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, xoff=0, yoff=0, xsize=1000, ysize=1000) else: image_dict[key] = layer_gdal_dataset.ReadAsArray() + image_dict['fmask_fill_value'] = metadata['_FillValue'] return True offset = 0.0 scale_factor = 1. - metadata = layer_gdal_dataset.GetMetadata() - if 'SPACECRAFT_NAME' not in dswx_metadata_dict.keys(): for k, v in metadata.items(): if k.upper() in METADATA_FIELDS_TO_COPY_FROM_HLS_LIST: @@ -1670,12 +1671,12 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, else: image = layer_gdal_dataset.ReadAsArray() - if fill_data is None and '_FillValue' in metadata.keys(): - fill_data = float(metadata['_FillValue']) - elif fill_data is None: - fill_data = -9999 + if fill_value is None and '_FillValue' in metadata.keys(): + fill_value = float(metadata['_FillValue']) + elif fill_value is None: + fill_value = -9999 - invalid_ind = np.where(image == fill_data) + invalid_ind = np.where(image == fill_value) if FLAG_CLIP_NEGATIVE_REFLECTANCE: image = np.clip(image, 1, None) if flag_offset_and_scale_inputs: @@ -1683,7 +1684,7 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, offset) image[invalid_ind] == np.nan elif FLAG_CLIP_NEGATIVE_REFLECTANCE: - image[invalid_ind] = fill_data + image[invalid_ind] = fill_value image_dict[key] = image @@ -1692,7 +1693,7 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, layer_gdal_dataset.GetGeoTransform() image_dict['projection'] = \ layer_gdal_dataset.GetProjection() - image_dict['fill_data'] = fill_data + image_dict['fill_value'] = fill_value image_dict['length'] = image_dict[key].shape[0] image_dict['width'] = image_dict[key].shape[1] @@ -3343,6 +3344,10 @@ def generate_dswx_layers(input_list, swir1 = image_dict['swir1'] swir2 = image_dict['swir2'] fmask = image_dict['fmask'] + fill_value = image_dict['fill_value'] + fmask_fill_value = image_dict['fmask_fill_value'] + print('*** fill_value:', fill_value) + print('*** fmask_fill_value:', fmask_fill_value) geotransform = image_dict['geotransform'] projection = image_dict['projection'] @@ -3436,9 +3441,11 @@ def generate_dswx_layers(input_list, # Set array of invalid pixels if not flag_offset_and_scale_inputs: - invalid_ind = np.where(blue == image_dict['fill_data']) + invalid_ind = np.where((blue == fill_value) & + (fmask == fmask_fill_value)) else: - invalid_ind = np.where(np.isnan(blue)) + invalid_ind = np.where((np.isnan(blue)) & + (fmask == fmask_fill_value)) if output_rgb_file: _save_output_rgb_file(red, green, blue, output_rgb_file, From 25ff66106e990bb1e818d0574e32c19ba0af974e Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Fri, 16 Dec 2022 09:44:35 -0800 Subject: [PATCH 13/21] remove debug messages --- src/proteus/dswx_hls.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index e83d7a3..4626dec 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -3346,8 +3346,6 @@ def generate_dswx_layers(input_list, fmask = image_dict['fmask'] fill_value = image_dict['fill_value'] fmask_fill_value = image_dict['fmask_fill_value'] - print('*** fill_value:', fill_value) - print('*** fmask_fill_value:', fmask_fill_value) geotransform = image_dict['geotransform'] projection = image_dict['projection'] From 8ea0ccd4f8a2b3b1210e885ab700db51aefa1692 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Fri, 16 Dec 2022 10:13:24 -0800 Subject: [PATCH 14/21] mask out points where Fmask == fill value (2) --- src/proteus/dswx_hls.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 4626dec..8f66478 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -111,7 +111,8 @@ 0b10000 : 4, 0b10001 : 4, 0b10010 : 4, - 0b10100 : 4 + 0b10100 : 4, + DIAGNOSTIC_LAYER_NO_DATA_DECIMAL: UINT8_FILL_VALUE } @@ -1595,7 +1596,7 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, xoff=0, yoff=0, xsize=1000, ysize=1000) else: image_dict[key] = layer_gdal_dataset.ReadAsArray() - image_dict['fmask_fill_value'] = metadata['_FillValue'] + image_dict['fmask_fill_value'] = int(metadata['_FillValue']) return True offset = 0.0 @@ -3439,10 +3440,10 @@ def generate_dswx_layers(input_list, # Set array of invalid pixels if not flag_offset_and_scale_inputs: - invalid_ind = np.where((blue == fill_value) & + invalid_ind = np.where((blue == fill_value) | (fmask == fmask_fill_value)) else: - invalid_ind = np.where((np.isnan(blue)) & + invalid_ind = np.where((np.isnan(blue)) | (fmask == fmask_fill_value)) if output_rgb_file: @@ -3468,12 +3469,10 @@ def generate_dswx_layers(input_list, diagnostic_layer_decimal = _compute_diagnostic_tests( blue, green, red, nir, swir1, swir2, hls_thresholds) - diagnostic_layer_decimal[invalid_ind] = DIAGNOSTIC_LAYER_NO_DATA_DECIMAL interpreted_dswx_band = generate_interpreted_layer( diagnostic_layer_decimal) - diagnostic_layer = _get_binary_representation(diagnostic_layer_decimal) del diagnostic_layer_decimal From 53d723491fd9e94bd241b4eb7af48973ced827c8 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Fri, 16 Dec 2022 18:25:37 -0800 Subject: [PATCH 15/21] add ctable to diagnostic layer --- src/proteus/dswx_hls.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 8f66478..71ba773 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -1896,6 +1896,40 @@ def _get_confidence_layer_ctable(): confidence_layer_ctable.SetColorEntry(UINT8_FILL_VALUE, (0, 0, 0, 255)) return confidence_layer_ctable + +def _get_diagnostic_layer_ctable(diagnostic_layer): + """ + Get diagnostic layer RGB color table + + Returns + ------- + diagnostic_layer_ctable : gdal.ColorTable + Diagnostic layer color table + """ + # create color table + diagnostic_layer_ctable = gdal.ColorTable() + + n_bits = 5 + unique_elements = np.unique(diagnostic_layer) + for element in unique_elements: + n_positive_tests = 0 + element_floor = element + for i in range(n_bits): + element_floor, remainder = np.divmod(element_floor, 10) + n_positive_tests += remainder + positive_tests_255 = int(float(n_positive_tests) * 255 // 5) + + diagnostic_layer_ctable.SetColorEntry( + int(element), (255 - positive_tests_255, + 255 - positive_tests_255, + 255)) + + # Transparent black - Fill value + diagnostic_layer_ctable.SetColorEntry(DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR, + (0, 0, 0, 255)) + return diagnostic_layer_ctable + + def _collapse_wtr_classes(interpreted_layer): """ Collapse interpreted layer classes onto final DSWx-HLS @@ -3477,12 +3511,15 @@ def generate_dswx_layers(input_list, del diagnostic_layer_decimal if output_diagnostic_layer: + diagnostic_layer_ctable = _get_diagnostic_layer_ctable( + diagnostic_layer) _save_array(diagnostic_layer, output_diagnostic_layer, dswx_metadata_dict, geotransform, projection, description=band_description_dict['DIAG'], scratch_dir=scratch_dir, output_files_list=build_vrt_list, output_dtype=gdal.GDT_UInt16, + ctable=diagnostic_layer_ctable, no_data_value=DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR) From 7d277bff61310585bca7f289b7bce18b61f71530 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Sat, 17 Dec 2022 09:25:28 -0800 Subject: [PATCH 16/21] fix assignment of DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR onto diagnostic_layer_binary --- src/proteus/dswx_hls.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 71ba773..9aff702 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -3107,7 +3107,6 @@ def _get_binary_representation(diagnostic_layer_decimal, nbits=6): diagnostic_layer_binary = np.zeros_like(diagnostic_layer_decimal, dtype=np.uint16) - for i in range(nbits): diagnostic_layer_decimal, bit_array = \ np.divmod(diagnostic_layer_decimal, 2) @@ -3115,7 +3114,7 @@ def _get_binary_representation(diagnostic_layer_decimal, nbits=6): diagnostic_layer_binary += bit_array * (10 ** i) else: # UInt16 max value is 65535 - diagnostic_layer_binary[bit_array] = \ + diagnostic_layer_binary[np.where(bit_array)] = \ DIAGNOSTIC_LAYER_NO_DATA_BINARY_REPR return diagnostic_layer_binary From d412e04865857f261c52ea141aa9c4fdfd9d1a06 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Mon, 19 Dec 2022 08:40:45 -0800 Subject: [PATCH 17/21] expose `copernicus_forest_classes` to the runconfig --- bin/dswx_hls.py | 1 + src/proteus/defaults/dswx_hls.yaml | 3 +++ src/proteus/dswx_hls.py | 32 ++++++++++++++++++++++++++---- src/proteus/schemas/dswx_hls.yaml | 3 +++ tests/test_dswx_hls_workflow.py | 1 + 5 files changed, 36 insertions(+), 4 deletions(-) diff --git a/bin/dswx_hls.py b/bin/dswx_hls.py index 146cae0..37db08c 100755 --- a/bin/dswx_hls.py +++ b/bin/dswx_hls.py @@ -80,6 +80,7 @@ def main(): max_sun_local_inc_angle=args.max_sun_local_inc_angle, apply_cast_shadow_masking=args.apply_cast_shadow_masking, mask_adjacent_to_cloud_mode=args.mask_adjacent_to_cloud_mode, + copernicus_forest_classes=args.copernicus_forest_classes, flag_debug=args.flag_debug) diff --git a/src/proteus/defaults/dswx_hls.yaml b/src/proteus/defaults/dswx_hls.yaml index b9766e1..44080a4 100644 --- a/src/proteus/defaults/dswx_hls.yaml +++ b/src/proteus/defaults/dswx_hls.yaml @@ -71,6 +71,9 @@ runconfig: # Define how areas adjacent to cloud/cloud-shadow should be handled mask_adjacent_to_cloud_mode: 'mask' + # Copernicus CGLS Land Cover 100m forest classes + copernicus_forest_classes: [20, 111, 113, 115, 116, 121, 123, 125, 126] + save_wtr: True # Layer 1 - WTR save_bwtr: True # Layer 2 - BWTR save_conf: True # Layer 3 - CONF diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 9aff702..40bf420 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -327,6 +327,8 @@ class RunConfigConstants: mask_adjacent_to_cloud_mode: str Define how areas adjacent to cloud/cloud-shadow should be handled. Options: "mask", "ignore", and "cover" + copernicus_forest_classes: list(int) + Copernicus CGLS Land Cover 100m forest classes browse_image_height: int Height in pixels of the browse image PNG browse_image_width: int @@ -339,9 +341,9 @@ def __init__(self): self.max_sun_local_inc_angle = None self.apply_cast_shadow_masking = None self.mask_adjacent_to_cloud_mode = None + self.copernicus_forest_classes = None self.browse_image_height = None self.browse_image_width = None - def get_dswx_hls_cli_parser(): @@ -551,6 +553,11 @@ def get_dswx_hls_cli_parser(): ' should be handled. Options: "mask", "ignore", and' ' "cover"') + parser.add_argument('--copernicus-forest-classes', + dest='copernicus_forest_classes', + type=list, + help='Copernicus CGLS Land Cover 100m forest classes') + parser.add_argument('--debug', dest='flag_debug', action='store_true', @@ -770,6 +777,7 @@ def _update_landcover_array(conglomerate_array, agg_sum, threshold, def create_landcover_mask(copernicus_landcover_file, worldcover_file, output_file, scratch_dir, mask_type, geotransform, projection, length, width, + copernicus_forest_classes, dswx_metadata_dict = None, output_files_list = None, temp_files_list = None): """ @@ -797,6 +805,8 @@ def create_landcover_mask(copernicus_landcover_file, DSWx-HLS product's length (number of lines) width: int DSWx-HLS product's width (number of columns) + copernicus_forest_classes: list(int) + Copernicus CGLS Land Cover 100m forest classes dswx_metadata_dict: dict (optional) Metadata dictionary that will store band metadata output_files_list: list (optional) @@ -876,8 +886,20 @@ def create_landcover_mask(copernicus_landcover_file, size_y, size_x) del tree_binary_mask - tree_aggregate_sum = np.where(copernicus_landcover_array == 111, - tree_aggregate_sum, 0) + copernicus_forest = np.zeros_like(tree_aggregate_sum, dtype=np.uint8) + if copernicus_forest_classes is None: + copernicus_forest_classes = [20, 111, 113, 115, 116, 121, 123, 125, + 126] + + logger.info(' CGLS Land Cover 100m forest classes:' + f' {copernicus_forest_classes}') + + for copernicus_forest_class in copernicus_forest_classes: + copernicus_forest |= (copernicus_landcover_array == + copernicus_forest_class) + + tree_aggregate_sum = np.where(copernicus_forest, tree_aggregate_sum, 0) + del copernicus_forest logger.info(f'combining masks') # create array filled with 30000 @@ -3173,6 +3195,7 @@ def generate_dswx_layers(input_list, max_sun_local_inc_angle=None, apply_cast_shadow_masking=None, mask_adjacent_to_cloud_mode=None, + copernicus_forest_classes=None, flag_debug=False): """Compute the DSWx-HLS product @@ -3468,7 +3491,8 @@ def generate_dswx_layers(input_list, landcover_mask = create_landcover_mask( landcover_file, worldcover_file, output_landcover, scratch_dir, landcover_mask_type, geotransform, projection, - length, width, dswx_metadata_dict = dswx_metadata_dict, + length, width, copernicus_forest_classes, + dswx_metadata_dict = dswx_metadata_dict, output_files_list=build_vrt_list, temp_files_list=temp_files_list) # Set array of invalid pixels diff --git a/src/proteus/schemas/dswx_hls.yaml b/src/proteus/schemas/dswx_hls.yaml index 8f6c049..64adc70 100644 --- a/src/proteus/schemas/dswx_hls.yaml +++ b/src/proteus/schemas/dswx_hls.yaml @@ -72,6 +72,9 @@ runconfig: # Define how areas adjacent to cloud/cloud-shadow should be handled mask_adjacent_to_cloud_mode: enum('mask', 'ignore', 'cover', required=False) + # Copernicus CGLS Land Cover 100m forest classes + copernicus_forest_classes: list(int(), min=0, required=False) + # Layer 1 - WTR save_wtr: bool(required=False) diff --git a/tests/test_dswx_hls_workflow.py b/tests/test_dswx_hls_workflow.py index dcacbe8..d1b8e50 100644 --- a/tests/test_dswx_hls_workflow.py +++ b/tests/test_dswx_hls_workflow.py @@ -92,6 +92,7 @@ def test_workflow(): max_sun_local_inc_angle=args.max_sun_local_inc_angle, apply_cast_shadow_masking=args.apply_cast_shadow_masking, mask_adjacent_to_cloud_mode=args.mask_adjacent_to_cloud_mode, + copernicus_forest_classes=args.copernicus_forest_classes, flag_debug=args.flag_debug) ref_files = glob.glob(os.path.join(ref_dir, '*')) From 264ed13821fb94aa00d35e1b60dc7f86b893578a Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Mon, 19 Dec 2022 10:57:38 -0800 Subject: [PATCH 18/21] update SOFTWARE_VERSION to v0.5.2 --- src/proteus/dswx_hls.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 40bf420..a80a4e5 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -16,7 +16,7 @@ from proteus.core import save_as_cog -SOFTWARE_VERSION = '0.5' +SOFTWARE_VERSION = '0.5.2' ''' Internally, DSWx-HLS SAS stores 4 water classes. The flag below enables or From 5d0cb1560d0f6e9c61691b7cf28b9fae79c6dcc3 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Mon, 19 Dec 2022 11:51:46 -0800 Subject: [PATCH 19/21] update PROCESSING_DATETIME and SPACECRAFT_NAME to follow DSWx-HLS product specs --- src/proteus/dswx_hls.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index a80a4e5..761fe02 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -1638,8 +1638,9 @@ def _load_hls_from_file(filename, image_dict, offset_dict, scale_dict, # HLS Sentinel metadata contain attribute SPACECRAFT_NAME if 'SPACECRAFT_NAME' in metadata: - spacecraft_name = metadata['SPACECRAFT_NAME'].upper() - if 'SENTINEL' not in spacecraft_name and 'LANDSAT' not in spacecraft_name: + spacecraft_name = metadata['SPACECRAFT_NAME'] + if ('SENTINEL' not in spacecraft_name.upper() and + 'LANDSAT' not in spacecraft_name.upper()): logger.info(f'ERROR the platform "{spacecraft_name}" is not supported') return False @@ -2797,7 +2798,7 @@ def _get_dswx_metadata_dict(product_id, product_version): # save datetime 'YYYY-MM-DD HH:MM:SS' dswx_metadata_dict['PROCESSING_DATETIME'] = \ - datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%S") + datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%SZ") return dswx_metadata_dict From 8a764b7b6e890d4c92b7a530c2b49aef7bb6fbfd Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Mon, 19 Dec 2022 12:36:49 -0800 Subject: [PATCH 20/21] store software version in a different file: src/proteus/version.py --- setup.py | 18 +++++++++++++++++- src/proteus/dswx_hls.py | 4 ++-- src/proteus/version.py | 1 + 3 files changed, 20 insertions(+), 3 deletions(-) create mode 100644 src/proteus/version.py diff --git a/setup.py b/setup.py index b3c86d7..b0c0ed5 100755 --- a/setup.py +++ b/setup.py @@ -2,7 +2,23 @@ from setuptools import setup from setuptools import Command -__version__ = version = VERSION = '0.1' + +def _get_version(): + """Returns the PROTEUS software version + + Returns + ------- + version : str + PROTEUS software version + """ + + version_file = 'src/proteus/version.py' + version = open(version_file).read().split('=')[1] + version = version.replace("'", '').replace('"', '').strip() + + return version + +__version__ = version = VERSION = _get_version() class CleanCommand(Command): diff --git a/src/proteus/dswx_hls.py b/src/proteus/dswx_hls.py index 761fe02..653a8a9 100644 --- a/src/proteus/dswx_hls.py +++ b/src/proteus/dswx_hls.py @@ -15,8 +15,7 @@ import scipy from proteus.core import save_as_cog - -SOFTWARE_VERSION = '0.5.2' +from proteus.version import VERSION as SOFTWARE_VERSION ''' Internally, DSWx-HLS SAS stores 4 water classes. The flag below enables or @@ -3319,6 +3318,7 @@ def generate_dswx_layers(input_list, elif product_id is None: product_id = 'dswx_hls' + logger.info(f'PROTEUS software version: {SOFTWARE_VERSION}') logger.info('input parameters:') logger.info(' file(s):') for input_file in input_list: diff --git a/src/proteus/version.py b/src/proteus/version.py new file mode 100644 index 0000000..c0741e1 --- /dev/null +++ b/src/proteus/version.py @@ -0,0 +1 @@ +VERSION = '0.5.2' \ No newline at end of file From 46f96cbf1b21eca48e941ff8c6cb1b128fe38e56 Mon Sep 17 00:00:00 2001 From: Gustavo Shiroma Date: Mon, 19 Dec 2022 12:40:40 -0800 Subject: [PATCH 21/21] store software version in a different file: src/proteus/version.py (2) --- src/proteus/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/proteus/version.py b/src/proteus/version.py index c0741e1..fe53b9c 100644 --- a/src/proteus/version.py +++ b/src/proteus/version.py @@ -1 +1 @@ -VERSION = '0.5.2' \ No newline at end of file +VERSION = '0.5.2'