From 3e0f9a5a2a5597a3c5a84f48abe53981eebb61fd Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Thu, 9 May 2024 22:35:57 -0700 Subject: [PATCH] Fix multilooking bug for ancilliary layers --- tools/ARIAtools/extractProduct.py | 38 +++++++++++++------------------ tools/ARIAtools/util/dem.py | 17 +++++--------- tools/ARIAtools/util/mask.py | 12 ++++------ tools/ARIAtools/util/vrt.py | 13 +---------- 4 files changed, 27 insertions(+), 53 deletions(-) diff --git a/tools/ARIAtools/extractProduct.py b/tools/ARIAtools/extractProduct.py index 8ffe96d4..a9e9671b 100644 --- a/tools/ARIAtools/extractProduct.py +++ b/tools/ARIAtools/extractProduct.py @@ -721,6 +721,10 @@ def handle_epoch_layers( sec_workdir = copy.deepcopy(workdir) ref_workdir = copy.deepcopy(workdir) + # Set output res + if multilooking is not None: + arrres = [arrres[0] * multilooking, arrres[1] * multilooking] + # If specified workdir doesn't exist, create it all_workdirs = [workdir, sec_workdir, ref_workdir] all_workdirs = list(set(all_workdirs)) @@ -873,7 +877,7 @@ def handle_epoch_layers( finalize_metadata( j[1][:-4], bounds, arrres, dem_bounds, prods_TOTbbox, dem, lat, lon, hgt_field, prod_ver_list, is_nisar_file, - mask, outputFormat, verbose=verbose) + outputFormat, verbose=verbose) # Track consistency of dimensions if j[0] == 0: @@ -954,8 +958,8 @@ def export_product_worker( # Interpolate/intersect with DEM before cropping finalize_metadata( outname, bounds, arrres, dem_bounds, prods_TOTbbox, dem_expanded, - lat, lon, hgt_field, product, is_nisar_file, mask, - outputFormatPhys, verbose=verbose) + lat, lon, hgt_field, product, is_nisar_file, outputFormatPhys, + verbose=verbose) # Extract/crop full res layers, except for "unw" and "conn_comp" # which requires advanced stitching @@ -1021,11 +1025,7 @@ def export_product_worker( osgeo.gdal.Open(j + '.vrt').ReadAsArray() update_file.GetRasterBand(1).WriteArray(mask_arr) - if layer != 'unwrappedPhase' and layer != 'connectedComponents' and \ - not any(":/science/grids/imagingGeometry" - in s for s in product) and \ - not any(":/science/grids/corrections" - in s for s in product): + if layer != 'unwrappedPhase' and layer != 'connectedComponents': # If necessary, resample raster if multilooking is not None: @@ -1262,8 +1262,14 @@ def export_products( prog_bar = ARIAtools.util.misc.ProgressBar( maxValue=len(product_dict[0]), prefix='Generating: ' + key + ' - ') + # Set output res + if multilooking is not None: + iono_arrres = [arrres[0] * multilooking, arrres[1] * multilooking] + else: + iono_arrres = arrres + lyr_input_dict = dict(input_iono_files=None, - arrres=arrres, + arrres=iono_arrres, epsg=epsg_code, output_iono=None, output_format=outputFormat, @@ -1390,8 +1396,7 @@ def export_products( def finalize_metadata(outname, bbox_bounds, arrres, dem_bounds, prods_TOTbbox, dem, lat, lon, hgt_field, prod_list, is_nisar_file=False, - mask=None, outputFormat='ENVI', verbose=None, - num_threads='2'): + outputFormat='ENVI', verbose=None, num_threads='2'): """Interpolate and extract 2D metadata layer. 2D metadata layer is derived by interpolating and then intersecting 3D layers with a DEM. @@ -1520,17 +1525,6 @@ def finalize_metadata(outname, bbox_bounds, arrres, dem_bounds, prods_TOTbbox, osgeo.gdal.Translate( outname + '.vrt', outname, options=translate_options) - # Apply mask (if specified) - if mask is not None: - out_interpolated = ( - mask.ReadAsArray() * osgeo.gdal.Open(outname).ReadAsArray()) - - # Update VRT with new raster - update_file = osgeo.gdal.Open(outname, osgeo.gdal.GA_Update) - update_file.GetRasterBand(1).WriteArray(out_interpolated) - update_file = None - out_interpolated = None - data_array = None diff --git a/tools/ARIAtools/util/dem.py b/tools/ARIAtools/util/dem.py index 64668bb7..bd8199b3 100644 --- a/tools/ARIAtools/util/dem.py +++ b/tools/ARIAtools/util/dem.py @@ -24,7 +24,7 @@ def prep_dem(demfilename, bbox_file, prods_TOTbbox, prods_TOTbbox_metadatalyr, proj, arrres=None, workdir='./', outputFormat='ENVI', num_threads='2', dem_name: str = 'glo_90', - multilooking=None, rankedResampling=False): + multilooking=1, rankedResampling=False): """ Function to load and export DEM, lat, lon arrays. If "Download" flag is specified, DEM will be downloaded on the fly. @@ -45,6 +45,9 @@ def prep_dem(demfilename, bbox_file, prods_TOTbbox, prods_TOTbbox_metadatalyr, "Cannot proceed with VRT format, using ENVI format instead") outputFormat = 'ENVI' + # Set output res + arrres = [arrres[0] * multilooking, arrres[1] * multilooking] + if demfilename.lower() == 'download': if dem_name not in dem_stitcher.datasets.DATASETS: raise ValueError( @@ -90,22 +93,14 @@ def prep_dem(demfilename, bbox_file, prods_TOTbbox, prods_TOTbbox_metadatalyr, 'Applied cutline to produce 3 arc-sec SRTM DEM: %s', aria_dem) - # Apply multilooking, if specified - if multilooking is not None: - ARIAtools.util.vrt.resampleRaster( - aria_dem, multilooking, bounds, prods_TOTbbox, rankedResampling, - outputFormat=outputFormat, num_threads=num_threads) - ds_aria = osgeo.gdal.Open(aria_dem) - # Load DEM and setup lat and lon arrays # pass expanded DEM for metadata field interpolation bounds = list( ARIAtools.util.shp.open_shp(prods_TOTbbox_metadatalyr).bounds) - gt = ds_aria.GetGeoTransform() gdal_warp_kwargs = { - 'format': outputFormat, 'outputBounds': bounds, 'xRes': abs(gt[1]), - 'yRes': abs(gt[-1]), 'targetAlignedPixels': True, 'multithread': True, + 'format': outputFormat, 'outputBounds': bounds, 'xRes': arrres[0], + 'yRes': arrres[1], 'targetAlignedPixels': True, 'multithread': True, 'options': ['NUM_THREADS=%s' % (num_threads) + ' -overwrite']} demfile_expanded = aria_dem.replace('.dem', '_expanded.dem') ds_aria_expanded = osgeo.gdal.Warp( diff --git a/tools/ARIAtools/util/mask.py b/tools/ARIAtools/util/mask.py index 9d7454c7..e4fa4707 100644 --- a/tools/ARIAtools/util/mask.py +++ b/tools/ARIAtools/util/mask.py @@ -27,7 +27,7 @@ def prep_mask( product_dict, maskfilename, bbox_file, prods_TOTbbox, proj, amp_thresh=None, arrres=None, workdir='./', outputFormat='ENVI', - num_threads='2', multilooking=None, rankedResampling=False): + num_threads='2', multilooking=1, rankedResampling=False): """ Function to load and export mask file with tile_mate """ @@ -45,6 +45,9 @@ def prep_mask( if outputFormat == 'VRT': outputFormat = 'ENVI' + # Set output res + arrres = [arrres[0] * multilooking, arrres[1] * multilooking] + # Download mask if maskfilename.lower() == 'download' or \ maskfilename.lower() in tile_mate.stitcher.DATASET_SHORTNAMES: @@ -181,13 +184,6 @@ def prep_mask( osgeo.gdal.Translate( maskfilename + '.vrt', maskfilename, options=translate_options) - # Apply multilooking, if specified - if multilooking is not None: - ARIAtools.util.vrt.resampleRaster( - maskfilename, multilooking, bounds, prods_TOTbbox, - rankedResampling, outputFormat=outputFormat, - num_threads=num_threads) - # remove temporary file for j in glob.glob(ref_file + '*'): if os.path.isfile(j): diff --git a/tools/ARIAtools/util/vrt.py b/tools/ARIAtools/util/vrt.py index 3e1154fe..a70f42e6 100644 --- a/tools/ARIAtools/util/vrt.py +++ b/tools/ARIAtools/util/vrt.py @@ -90,20 +90,9 @@ def resampleRaster( geotrans = ds.GetGeoTransform() proj = ds.GetProjection() - # Must resample mask with nearest-neighbor - if fname.split('/')[-2] == 'mask': - - # Resample raster - warp_options = osgeo.gdal.WarpOptions( - format=outputFormat, cutlineDSName=prods_TOTbbox, - outputBounds=bounds, xRes=arrres[0], yRes=arrres[1], - targetAlignedPixels=True, resampleAlg='near', multithread=True, - options=['NUM_THREADS=%s' % (num_threads) + ' -overwrite']) - osgeo.gdal.Warp(fname, inputname, options=warp_options) - # Use pixel function to downsample connected components/unw files # based off of frequency of connected components in each window - elif fname.split('/')[-2] == 'connectedComponents' \ + if fname.split('/')[-2] == 'connectedComponents' \ or fname.split('/')[-2] == 'unwrappedPhase': # Resample unw phase based off of mode of connected components