Fix multilooking bug for ancilliary layers

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
 
 

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(

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):

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