From ef48728c4303fec1a7e6e1e355f8c378428d537c Mon Sep 17 00:00:00 2001 From: Simran S Sangha Date: Sat, 5 Aug 2023 19:39:19 -0700 Subject: [PATCH] Fix multilooking bug --- tools/ARIAtools/extractProduct.py | 129 +++++++++++++++++------------- tools/ARIAtools/mask_util.py | 15 +++- tools/ARIAtools/tsSetup.py | 59 +++++++------- tools/ARIAtools/vrtmanager.py | 38 --------- 4 files changed, 119 insertions(+), 122 deletions(-) diff --git a/tools/ARIAtools/extractProduct.py b/tools/ARIAtools/extractProduct.py index 1d226168..2e1d7cfa 100755 --- a/tools/ARIAtools/extractProduct.py +++ b/tools/ARIAtools/extractProduct.py @@ -28,7 +28,7 @@ from ARIAtools.unwrapStitching import product_stitch_overlap, \ product_stitch_2stage from ARIAtools.vrtmanager import renderVRT, resampleRaster, layerCheck, \ - get_basic_attrs, ancillaryLooks, dim_check + get_basic_attrs, dim_check from ARIAtools.sequential_stitching import product_stitch_sequential import pyproj from pyproj import CRS, Transformer @@ -366,7 +366,8 @@ def __run__(self): 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'): + outputFormat='ENVI', num_threads='2', dem_name: str = 'glo_90', + multilooking=None, rankedResampling=False): """Function to load and export DEM, lat, lon arrays. If "Download" flag is specified, DEM will be downloaded on the fly. """ @@ -422,6 +423,14 @@ def prep_dem(demfilename, bbox_file, prods_TOTbbox, prods_TOTbbox_metadatalyr, ds_aria = gdal.Translate(f'{aria_dem}.vrt', aria_dem, format='VRT') log.info('Applied cutline to produce 3 arc-sec SRTM DEM: %s', aria_dem) + # Apply multilooking, if specified + if multilooking is not None: + resampleRaster(aria_dem, multilooking, + bounds, prods_TOTbbox, rankedResampling, + outputFormat=outputFormat, + num_threads=num_threads) + ds_aria = gdal.Open(aria_dem, gdal.GA_ReadOnly) + # Load DEM and setup lat and lon arrays # pass expanded DEM for metadata field interpolation bounds = list(open_shapefile(prods_TOTbbox_metadatalyr, 0, 0).bounds) @@ -936,11 +945,6 @@ def handle_epoch_layers(layers, prods_TOTbbox, dem, lat, lon, hgt_field, prod_ver_list, mask, outputFormat, verbose=verbose) - # If necessary, resample raster - if multilooking is not None: - resampleRaster(j[1][:-4], multilooking, bounds, - prods_TOTbbox, rankedResampling, - outputFormat=outputFormat, num_threads=num_threads) # Track consistency of dimensions if j[0] == 0: @@ -1149,7 +1153,7 @@ def export_products(full_product_dict, bbox_file, prods_TOTbbox, layers, in s for s in i[1]): # make VRT pointing to metadata layers in standard product hgt_field, model_name, outname = prep_metadatalayers(outname, - i[1], dem, key, layers) + i[1], dem, key, layers) # Interpolate/intersect with DEM before cropping finalize_metadata(outname, bounds, dem_bounds, @@ -1179,14 +1183,6 @@ def export_products(full_product_dict, bbox_file, prods_TOTbbox, layers, gdal.Translate(outname+'.vrt', outname, options=gdal.TranslateOptions(format="VRT")) - # Apply mask (if specified). - if mask is not None: - update_file = gdal.Open(outname, gdal.GA_Update) - mask_arr = mask.ReadAsArray() * \ - gdal.Open(outname + '.vrt').ReadAsArray() - update_file.GetRasterBand(1).WriteArray(mask_arr) - del update_file, mask_arr - # Extract/crop phs and conn_comp layers else: # get connected component input files @@ -1209,7 +1205,7 @@ def export_products(full_product_dict, bbox_file, prods_TOTbbox, layers, prods_TOTbbox, outFileUnw=outFilePhs, outFileConnComp=outFileConnComp, - mask=mask, + #mask=mask, outputFormat=outputFormatPhys, verbose=verbose) @@ -1221,7 +1217,7 @@ def export_products(full_product_dict, bbox_file, prods_TOTbbox, layers, prods_TOTbbox, outFileUnw=outFilePhs, outFileConnComp=outFileConnComp, - mask=mask, + #mask=mask, outputFormat=outputFormatPhys, verbose=verbose) @@ -1233,7 +1229,7 @@ def export_products(full_product_dict, bbox_file, prods_TOTbbox, layers, clip_json=prods_TOTbbox, output_unw=outFilePhs, output_conn=outFileConnComp, - mask_file=mask, # str filename + #mask_file=mask, # str filename output_format=outputFormatPhys, range_correction=True, save_fig=False, @@ -1246,15 +1242,35 @@ def export_products(full_product_dict, bbox_file, prods_TOTbbox, layers, prods_TOTbbox, rankedResampling, outputFormat=outputFormatPhys, num_threads=num_threads) - - # If necessary, resample raster - if multilooking is not None and \ - key != 'unwrappedPhase' and \ - key != 'connectedComponents': - resampleRaster(outname, multilooking, bounds, prods_TOTbbox, - rankedResampling, - outputFormat=outputFormatPhys, - num_threads=num_threads) + # Apply mask (if specified) + if mask is not None: + for j in [outFileConnComp, outFilePhs]: + update_file = gdal.Open(j, gdal.GA_Update) + mask_arr = mask.ReadAsArray() * \ + gdal.Open(j + '.vrt').ReadAsArray() + update_file.GetRasterBand(1).WriteArray(mask_arr) + del update_file, mask_arr + + if key != 'unwrappedPhase' and \ + key != 'connectedComponents' and \ + not any(":/science/grids/imagingGeometry" + in s for s in i[1]) and \ + not any(":/science/grids/corrections" + in s for s in i[1]): + # If necessary, resample raster + if multilooking is not None: + resampleRaster(outname, multilooking, bounds, + prods_TOTbbox, + rankedResampling, + outputFormat=outputFormatPhys, + num_threads=num_threads) + # Apply mask (if specified) + if mask is not None: + update_file = gdal.Open(outname, gdal.GA_Update) + mask_arr = mask.ReadAsArray() * \ + gdal.Open(outname + '.vrt').ReadAsArray() + update_file.GetRasterBand(1).WriteArray(mask_arr) + del update_file, mask_arr # Track consistency of dimensions if key_ind == 0: @@ -1880,24 +1896,43 @@ def main(inps=None): if 'amplitude' in d: for item in list(set(d['amplitude'])): amplitude_products.append(item) - inps.mask = prep_mask(amplitude_products, inps.mask, - standardproduct_info.bbox_file, - prods_TOTbbox, proj, amp_thresh=inps.amp_thresh, - arrres=arrres, - workdir=inps.workdir, - outputFormat=inps.outputFormat, - num_threads=inps.num_threads) + # mask parms + mask_dict = { + 'product_dict': amplitude_products, + 'maskfilename': inps.mask, + 'bbox_file': standardproduct_info.bbox_file, + 'prods_TOTbbox': prods_TOTbbox, + 'proj': proj, + 'amp_thresh': inps.amp_thresh, + 'arrres': arrres, + 'workdir': inps.workdir, + 'outputFormat': inps.outputFormat, + 'num_threads': inps.num_threads, + 'multilooking': inps.multilooking, + 'rankedResampling': inps.rankedResampling + } + inps.mask = prep_mask(**mask_dict) # Download/Load DEM & Lat/Lon arrays, providing bbox, # expected DEM shape, and output dir as input. if inps.demfile is not None: print('Download/cropping DEM') + # DEM parms + dem_dict = { + 'demfilename': inps.demfile, + 'bbox_file': standardproduct_info.bbox_file, + 'prods_TOTbbox': prods_TOTbbox, + 'prods_TOTbbox_metadatalyr': prods_TOTbbox_metadatalyr, + 'proj': proj, + 'arrres': arrres, + 'workdir': inps.workdir, + 'outputFormat': inps.outputFormat, + 'num_threads': inps.num_threads, + 'multilooking': inps.multilooking, + 'rankedResampling': inps.rankedResampling + } # Pass DEM-filename, loaded DEM array, and lat/lon arrays - inps.demfile, demfile, Latitude, Longitude = prep_dem( - inps.demfile, standardproduct_info.bbox_file, - prods_TOTbbox, prods_TOTbbox_metadatalyr, proj, - arrres=arrres, workdir=inps.workdir, - outputFormat=inps.outputFormat, num_threads=inps.num_threads) + inps.demfile, demfile, Latitude, Longitude = prep_dem(**dem_dict) else: demfile, Latitude, Longitude = None, None, None @@ -1927,20 +1962,6 @@ def main(inps=None): # Extract user expected layers arrshape = export_products(**export_dict) - # If necessary, resample DEM/mask AFTER they have been used to extract - ancillary_dict = { - 'mask': inps.mask, - 'dem': demfile, - 'arrshape': arrshape, - 'standardproduct_info': standardproduct_info, - 'multilooking': inps.multilooking, - 'prods_TOTbbox': prods_TOTbbox, - 'rankedResampling': inps.rankedResampling, - 'outputFormat': inps.outputFormat, - 'num_threads': inps.num_threads - } - ancillaryLooks(**ancillary_dict) - # Perform GACOS-based tropospheric corrections (if specified). if inps.gacos_products: gacos_correction(standardproduct_info.products, inps.gacos_products, diff --git a/tools/ARIAtools/mask_util.py b/tools/ARIAtools/mask_util.py index dbb42fd8..f78d8a81 100755 --- a/tools/ARIAtools/mask_util.py +++ b/tools/ARIAtools/mask_util.py @@ -18,7 +18,7 @@ from osgeo import gdal, ogr, gdalconst from ARIAtools.shapefile_util import open_shapefile -from ARIAtools.vrtmanager import rasterAverage +from ARIAtools.vrtmanager import rasterAverage, resampleRaster from ARIAtools.logger import logger log = logging.getLogger(__name__) @@ -26,8 +26,9 @@ def prep_mask(product_dict, maskfilename, bbox_file, prods_TOTbbox, proj, - amp_thresh=None, arrres=None, workdir='./', - outputFormat='ENVI', num_threads='2'): + amp_thresh=None, arrres=None, workdir='./', + outputFormat='ENVI', num_threads='2', + multilooking=None, rankedResampling=False): """Function to load and export mask file. If "Download" flag, GSHHS water mask will be donwloaded on the fly. If full resolution NLCD landcover data is given (NLCD...img) it gets cropped @@ -168,6 +169,7 @@ def prep_mask(product_dict, maskfilename, bbox_file, prods_TOTbbox, proj, xRes=arrres[0], yRes=arrres[1], targetAlignedPixels=True, multithread=True, options=[f'NUM_THREADS={num_threads} -overwrite']) + mask = gdal.Open(maskfilename, gdal.GA_Update) mask.SetProjection(proj) mask.SetDescription(maskfilename) @@ -179,6 +181,13 @@ def prep_mask(product_dict, maskfilename, bbox_file, prods_TOTbbox, proj, gdal.Translate(maskfilename+'.vrt', maskfilename, options=gdal.TranslateOptions(format="VRT")) + # Apply multilooking, if specified + if multilooking is not None: + resampleRaster(maskfilename, multilooking, + bounds, prods_TOTbbox, rankedResampling, + outputFormat=outputFormat, + num_threads=num_threads) + # pass maskfile object mask = gdal.Open(maskfilename) diff --git a/tools/ARIAtools/tsSetup.py b/tools/ARIAtools/tsSetup.py index eb9a63ee..abbeb21e 100644 --- a/tools/ARIAtools/tsSetup.py +++ b/tools/ARIAtools/tsSetup.py @@ -29,7 +29,7 @@ from ARIAtools.mask_util import prep_mask from ARIAtools.shapefile_util import open_shapefile from ARIAtools.vrtmanager import resampleRaster, layerCheck, \ - get_basic_attrs, ancillaryLooks, dim_check + get_basic_attrs, dim_check from ARIAtools.extractProduct import merged_productbbox, prep_dem, \ export_products, gacos_correction @@ -427,12 +427,22 @@ def main(inps=None): # expected DEM shape, and output dir as input. if inps.demfile is not None: print('Download/cropping DEM') + # DEM parms + dem_dict = { + 'demfilename': inps.demfile, + 'bbox_file': standardproduct_info.bbox_file, + 'prods_TOTbbox': prods_TOTbbox, + 'prods_TOTbbox_metadatalyr': prods_TOTbbox_metadatalyr, + 'proj': proj, + 'arrres': arrres, + 'workdir': inps.workdir, + 'outputFormat': inps.outputFormat, + 'num_threads': inps.num_threads, + 'multilooking': inps.multilooking, + 'rankedResampling': inps.rankedResampling + } # Pass DEM-filename, loaded DEM array, and lat/lon arrays - inps.demfile, demfile, Latitude, Longitude = prep_dem( - inps.demfile, standardproduct_info.bbox_file, - prods_TOTbbox, prods_TOTbbox_metadatalyr, proj, - arrres=arrres, workdir=inps.workdir, - outputFormat=inps.outputFormat, num_threads=inps.num_threads) + inps.demfile, demfile, Latitude, Longitude = prep_dem(**dem_dict) # Load or download mask (if specified). if inps.mask is not None: @@ -442,13 +452,22 @@ def main(inps=None): if 'amplitude' in d: for item in list(set(d['amplitude'])): amplitude_products.append(item) - inps.mask = prep_mask(amplitude_products, inps.mask, - standardproduct_info.bbox_file, - prods_TOTbbox, proj, amp_thresh=inps.amp_thresh, - arrres=arrres, - workdir=inps.workdir, - outputFormat=inps.outputFormat, - num_threads=inps.num_threads) + # mask parms + mask_dict = { + 'product_dict': amplitude_products, + 'maskfilename': inps.mask, + 'bbox_file': standardproduct_info.bbox_file, + 'prods_TOTbbox': prods_TOTbbox, + 'proj': proj, + 'amp_thresh': inps.amp_thresh, + 'arrres': arrres, + 'workdir': inps.workdir, + 'outputFormat': inps.outputFormat, + 'num_threads': inps.num_threads, + 'multilooking': inps.multilooking, + 'rankedResampling': inps.rankedResampling + } + inps.mask = prep_mask(**mask_dict) # Extract # aria_extract default parms @@ -529,20 +548,6 @@ def main(inps=None): # Track consistency of dimensions dim_check(ref_arr_record, prod_arr_record) - # If necessary, resample DEM/mask AFTER they have been used to extract - ancillary_dict = { - 'mask': inps.mask, - 'dem': demfile, - 'arrshape': ref_arr_record, - 'standardproduct_info': standardproduct_info, - 'multilooking': inps.multilooking, - 'prods_TOTbbox': prods_TOTbbox, - 'rankedResampling': inps.rankedResampling, - 'outputFormat': inps.outputFormat, - 'num_threads': inps.num_threads - } - ancillaryLooks(**ancillary_dict) - # Perform GACOS-based tropospheric corrections (if specified). if inps.gacos_products: gacos_correction(standardproduct_info.products, inps.gacos_products, diff --git a/tools/ARIAtools/vrtmanager.py b/tools/ARIAtools/vrtmanager.py index 2a9763ab..18325d7a 100755 --- a/tools/ARIAtools/vrtmanager.py +++ b/tools/ARIAtools/vrtmanager.py @@ -263,44 +263,6 @@ def resampleRaster(fname, multilooking, bounds, prods_TOTbbox, return -###Check ancillary layers' -def ancillaryLooks(mask, dem, arrshape, standardproduct_info, multilooking, - prods_TOTbbox, rankedResampling, outputFormat='ENVI', - num_threads='2'): - """Check and apply looks on ancillary layers, if necessary""" - # Cannot proceed with VRT format. Defaulting to ENVI format. - if outputFormat=='VRT': - outputFormat='ENVI' - - # If necessary, resample DEM/mask AFTER they have been used to extract - # metadata layers and mask output layers, respectively - if multilooking is not None: - ref_height = arrshape[0] - ref_wid = arrshape[1] - bounds = open_shapefile(standardproduct_info.bbox_file, 0, 0).bounds - # Resample mask - if mask is not None: - prod_wid, prod_height, ref_geotrans, \ - _, _ = get_basic_attrs(mask.GetDescription()) - if (ref_wid != prod_wid) or (ref_height != prod_height): - resampleRaster(mask.GetDescription(), multilooking, - bounds, prods_TOTbbox, rankedResampling, - outputFormat=outputFormat, - num_threads=num_threads) - # Resample DEM - if dem is not None: - prod_wid, prod_height, prod_geotrans, \ - _, _ = get_basic_attrs(dem.GetDescription()) - if (ref_wid != prod_wid) or (ref_height != prod_height) or \ - (ref_geotrans != prod_geotrans): - resampleRaster(dem.GetDescription(), multilooking, - bounds, prods_TOTbbox, rankedResampling, - outputFormat=outputFormat, - num_threads=num_threads) - - return - - ###Average rasters def rasterAverage(outname, product_dict, bounds, prods_TOTbbox, arrres, outputFormat='ENVI', thresh=None):