Skip to content

Commit

Permalink
Merge pull request #377 from aria-tools/sss_lksbugfix
Browse files Browse the repository at this point in the history 
Fix multilooking bug
  • Loading branch information
@sssangha
sssangha authored Aug 8, 2023
2 parents 6c5eea5 + ef48728 commit 0de2354
Show file tree
Hide file tree
Showing 4 changed files with 119 additions and 122 deletions.
129 changes: 75 additions & 54 deletions tools/ARIAtools/extractProduct.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
"""
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand All @@ -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)

Expand All @@ -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,
Expand All @@ -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:
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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,
Expand Down
15 changes: 12 additions & 3 deletions tools/ARIAtools/mask_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,16 +18,17 @@
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__)
logger.setLevel(logging.DEBUG)


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
Expand Down Expand Up @@ -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)
Expand All @@ -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)

Expand Down
59 changes: 32 additions & 27 deletions tools/ARIAtools/tsSetup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand Down
38 changes: 0 additions & 38 deletions tools/ARIAtools/vrtmanager.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down

0 comments on commit 0de2354

Please sign in to comment.