From 566c3072c5cbbb5abf20969c098353ed0b05e523 Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 11 Apr 2022 17:09:32 -0700 Subject: [PATCH 1/7] Unwrapping script to use isce3 modules for Phass and Snaphu v2.0.4 This PR allows using isce3 modules for unwrapping with Phass and Snaphu algorithms. Script unwrap_isce3.py contains functions: 1. unwrap_phass - prepares and runs unwrapping with Phass, as Phass takes wrapped interferogram input in radians, function vrt_cpx2rad creates vrt file that extracts phase values on-the-fly using pixel function: phase 2. unwrap_snaphu - prepares and runs unwrapping with Snaphu version 2.0.4, using tiling option useful for unwrapping high-resolution interferograms. Added options related to tiling, such as number of tiles per row and column, number of processes, minimum region and etc.... and deformation_threshold option which I find useful when snaphu ends up in the infite loop due to a cost overflow (an unexpected increase in the total solution cost), examples Ridgecrest earthquake coseismic pairs. snaphu manpage suggests increasing the value of DEFOTHRESHFACTOR or decreasing the value of DEFOCONST. Snaphu standalone version gives this warning but it appers Python bindings for snaphu, does not, they just end up in very long with segmentation fault at the end 3. unw_vrt_xml - function that creates .xml and .vrt files for output in the similar manner as isceobj.createImage() but without dependence on the isce2 --- python/unwrap_isce3.py | 437 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 437 insertions(+) create mode 100644 python/unwrap_isce3.py diff --git a/python/unwrap_isce3.py b/python/unwrap_isce3.py new file mode 100644 index 0000000..ac60778 --- /dev/null +++ b/python/unwrap_isce3.py @@ -0,0 +1,437 @@ +#! /usr/bin/env python3 + +# Author: Marin Govorcin +# Caltech JPL, Dec 05, 2021 + +import os +import argparse +import isce3 +import numpy as np +from isce3.io.gdal import Raster +import xml.etree.ElementTree as ET +from xml.dom import minidom + +class ParseKwargs(argparse.Action): + def __call__(self, parser, namespace, values, option_string=None): + setattr(namespace, self.dest, dict()) + for value in values: + key, value = value.split('=') + getattr(namespace, self.dest)[key] = value + +def cmdLineParser(): + ''' + Command Line Parser + Script to use Phass/Snaphu in the isce3 environment + + USAGE: + Phass + unwrap_isce3.py -i 20190610_20190716.int -c tcorr_ds_ps.bin -k corr_thresh=0.3, good_corr=0.7 + unwrap_isce3.py -i 20190610_20190716.int -c tcorr_ds_ps.bin -k corr_thresh=0.3, good_corr=0.7 min_region=200 + + Snaphu: + unwrap_isce3.py -i 20190610_20190716.int -c tcorr_ds_ps.bin -u snaphu -k nlooks=27 tiles=[4,4,500,500] nproc=8 min_region=200 + ''' + unw_options = ''' phass :: corr_thresh=0.3, good_corr=0.7, min_region=200, \n + snaphu :: nlooks=1, tiles=[1, 1, 0, 0], min_region=200, nproc=1, single_tile_opt=False ''' + + parser = argparse.ArgumentParser(description='Phass Unwrapper', formatter_class=argparse.RawTextHelpFormatter) + parser.add_argument('-i', '--interferogram_file', type=str, dest='igram', + required=True, help='Input interferogram file in complex format, example: file.int') + parser.add_argument('-c', '--correlation_file', type=str, dest='correlation', + required=True, help='Input correlation file, example: tcorr_ds_ps.bin') + parser.add_argument('-m', '--mask', type=str, dest='mask', default=None, + help='Input mask byte file') + parser.add_argument('-o', '--output_dir', type=str, dest='outputDir', + default='./', help='Output directory path') + parser.add_argument('-u', '--unwrapper', type=str, dest='unwrapper', + default='phass', help='Unwrapper method: phass/snaphu') + parser.add_argument('-k', '--kwargs', type=str, dest='kwargs', nargs='*', + help=unw_options, action=ParseKwargs) + + return parser.parse_args() + +############################################################################################################ +############################ UTILITS FOR UNWRAPPING OUTPUT METDADATA ####################################### +############################################################################################################ + +GDAL_DATATYPE = { + 1 : "Byte", + 2 : "UInt16", + 3 : "Int16", + 4 : "UInt32", + 5 : "UInt64", + 6 : "Float32", + 7 : "Float64", + 8 : "CInt16", + 9 : "CInt32", + 10 : "CFloat32", + 11 : "CFloat64" + } + +GDAL2ISCE = { + 1 : 'BYTE', + 2 : 'BYTE', + 3 : 'SHORT', + 5 : 'INT', + 6 : 'FLOAT', + 7 : 'DOUBLE', + 10 : 'CFLOAT', + 11 : 'CFLOAT' + } + +def write_vrt(infile, width, length, dtype): + ## infile :: str path to inputfile + + infile_dir = os.path.dirname(os.path.abspath(infile)) + + vrttmpl=''' + + {PATH} + 0 + 4 + {linewidth} + LSB + +''' + + outfile = os.path.join(infile_dir, '{0}.vrt'.format(infile)) + + with open(outfile, 'w') as fid: + fid.write(vrttmpl.format(width = width, + length = length, + DTYPE = GDAL_DATATYPE[dtype], + PATH = infile, + linewidth = 4 * width)) + +def property_xml(root, name, value, doc): + property = ET.SubElement(root, 'property', {'name':name}) + ET.SubElement(property, 'value').text = value + ET.SubElement(property, 'doc').text=doc + +def component_xml(root, name, size, doc): + coord = ET.SubElement(root, 'component', {'name':name}) + ET.SubElement(coord, 'factorymodule').text = 'isceobj.Image' + ET.SubElement(coord, 'factoryname').text = 'createCoordinate' + ET.SubElement(coord, 'doc').text = doc + property_xml(coord, 'delta', '1.0', 'Coordinate quantization.') + property_xml(coord, 'endingvalues', str(float(size)), 'Ending value of the coordinate.') + property_xml(coord, 'family', 'imagecoordinate', 'Instance family name.') + property_xml(coord, 'name', 'imagecoordinate_name', 'Instance name.') + property_xml(coord, 'size', str(width), 'Coordinate size.') + property_xml(coord, 'startingvalue', '0.0', 'Starting value of the coordinate.') + +def write_xml(infile, width, length, bands, datatype, scheme): + + root = ET.Element('imageFile') + property_xml(root, 'access_mode', 'READ', 'Image access mode.') + property_xml(root, 'byte_order', 'l', 'Endianness of the image.') + component_xml(root, 'coordinate1', width, 'First coordinate of a 2D image (width).') + component_xml(root, 'coordinate2', length, 'Second coordinate of a 2D image (length).') + property_xml(root, 'data_type', GDAL2ISCE[dtype], 'Image data type.') + property_xml(root, 'family', 'image', 'Instance family name.') + property_xml(root, 'file_name', os.path.abspath(infile), 'Name of the image file.') + property_xml(root, 'length', str(length), 'Image length.') + property_xml(root, 'name', 'image_name', 'Instance name.') + property_xml(root, 'number_bands', str(bands), 'Number of image bands.') + property_xml(root, 'scheme', scheme, 'Interleaving scheme of the image.') + property_xml(root, 'width', str(width), 'Image width.') + property_xml(root, 'x_max', str(width), 'Maximum range value.') + property_xml(root, 'x_min', '0.0', 'Minimum range value.') + + #### Add ident to xml file + xmlstr = minidom.parseString(ET.tostring(root)).toprettyxml(indent=" ") + with open(infile+'.xml', "w") as f: + f.write(xmlstr) + +def unw_vrt_xml(infile, scheme): + output = Raster(infile) + width = output.width + length = output.length + dtype = output.datatype + bands = output.band + + write_vrt(infile, width, length, dtype) + write_xml(infile, width, length, bands, dtype, scheme) + + +############################################################################################################ +################################### UNWRAPPING PHASS/SNAPHU ######################################### +############################################################################################################ + +def vrt_cpx2rad(infile, width, length, stype): + ''' + Function to extract interferogram phase from complex file + needed when unwrapping interferogram in radians. + + Find the input interferogram dtype, if complex32 or complex64 + reate vrt file that extract phase values onthefly + + Input: + infile :: str wrapped interferogram filename + Output + outfile :: str output .vrt file + ''' + + # Get path to directory and name of inputfile + infile_dir = os.path.dirname(os.path.abspath(infile)) + infile_name = os.path.splitext(os.path.basename(infile))[0] + + #Input file is COMPLEX, create VRT with default pixelfunction : phase to extract phase in radian + if stype == 10 or stype == 11: + + vrttmpl=''' + + Wrapped interferogram in radians on-the-fly + + {PATH} + 0 + 4 + {linewidth} + LSB + + phase + {STYPE} + +''' + + if stype == 10: + dtype = 6 + elif stype == 11: + dtype = 7 + + outfile = os.path.join(infile_dir, '{0}.phase.vrt'.format(infile_name)) + + with open(outfile, 'w') as fid: + fid.write(vrttmpl.format(width = width, + length = length, + DTYPE = GDAL_DATATYPE[dtype], + STYPE = GDAL_DATATYPE[stype], + PATH = infile, + linewidth = 4 * width)) + + else: + outfile = os.path.abspath(infile) + + # Return the path to VRT file as input to Phass + return outfile + +def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, + correlation_threshold=0.3, good_correlation=0.7, min_region=200, print_msg=True): + ''' + Unwrap wrapped interferogram [in radians] using Phass + + Input: + infile :: str Wrapped phase in radians [interferogram *.int in PS_DS folder for fringe], input format .vrt + corr_file :: str Correlation file [tcorr_ds_ps.bin for fringe], input format .vrt + unw_filename :: str Output filename of unwrapped interferogram + conncomp_filename :: str Output filename of connected components + correlation_threshold :: float Threshold for minimum coherence + good_correlation :: float Good coherence + min_region :: int Minimum number of pixels per region (connected component) + + Output: + unw_igram :: Raster, Unwrapped interferogram, type: isce3.io.Raster + conn_comp :: Raster, Connected component labels, type: isce3.io.Raster + ''' + + from isce3.io import Raster # Phass module takes isce3.io.Raster inputs + from isce3.unwrap import Phass + + phass = Phass() + + # If input interferogram is in COMPLEX format, create .vrt file + # to get interferogram phase in radians on-the-fly + infile_name = infile + infile = Raster(infile) + phass_input = vrt_cpx2rad(infile_name, infile.width, infile.length, infile.datatype()) + + # Input files + igram = Raster(phass_input) + correlation = Raster(corr_file) + + # Output files + unw_igram = Raster(unw_filename, igram.width, igram.length,1, 6, "ENVI") #works also with ISCE format + conn_comp = Raster(conncomp_filename, igram.width, igram.length, 1, 4, "ENVI") + + # Phass unwrapping parameters + phass.correlation_threshold = correlation_threshold + phass.good_correlation = good_correlation + phass.min_pixels_region = min_region + + if print_msg: + print('\nPhass Unwrapping:') + print('Input:') + print(' Phase: {}'.format(os.path.abspath(phase_input))) + print(' width: {}, length {}, dtype: {}'.format(igram.width, igram.length, GDAL_DATATYPE[igram.datatype()])) + print(' Correlation: {}'.format(os.path.abspath(corr_file))) + print(' width: {}, length {}, dtype: {}'.format(correlation.width, correlation.length, GDAL_DATATYPE[correlation.datatype()])) + print(' Parameters: correlation_threshold={}, good_correlation={}, min_region={}'.format(correlation_threshold, good_correlation, min_region)) + + # Unwrap + phass.unwrap(igram, correlation, unw_igram, conn_comp) + + if print_msg: + print('#############################################') + print('Finished unwrapping, output files:') + print(' Unwrapped interferogram: {}'.format(os.path.abspath(unw_filename))) + print(' width: {}, length {}, dtype: {}'.format(unw_igram.width, unw_igram.length, GDAL_DATATYPE[unw_igram.datatype()])) + print(' Connected components: {}'.format(os.path.abspath(conncomp_filename))) + print(' width: {}, length {}, dtype: {}'.format(conn_comp.width, conn_comp.length, GDAL_DATATYPE[conn_comp.datatype()])) + +def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, + cost_mode='defo', nlooks=1, tiles=[1,1,0,0], tile_cost_thresh=500, min_region=100, + nproc=1, single_tile_opt=False, defo_thresh=1.25, print_msg=True): + + ''' + Unwrap wrapped interferogram using Snaphu + + Input: + igram :: str Input interferogram. Must be GDT_CFloat32 datatype + corr_file :: str Correlation magnitude, normalized to the interval [0, 1], GDT_Float32 datatype + unw_filename :: str Output filename of unwrapped interferogram + conncomp_filename :: str Output filename of connected components + mask :: str Binary mask of valid pixels. Zeros in this raster indicate interferogram + pixels that should be masked out, GDT_Byte datatype + nlooks :: float Effective number of looks used to form the input correlation data + tiles :: list [tile_nrows, tile_ncols, row_overlap, col_overlap] + tile_nrows, tile_ncols - Number of tiles along the row/column directions + row_overlap, col_overlap - verlap, in number of rows/columns, between neighboring tiles + tile_cost_thresh :: int Cost threshold to use for determining boundaries of reliable regions. + Larger cost threshold implies smaller regions (safer, but more expensive computationally + min_region :: int Minimum number of pixels per region (connected component) + nproc :: int Number of threads for parallel processing (use with tiles) + single_tile_optv :: bool Single Tile Reoptimization + defo_thresh :: float Factor applied to rho0 to get threshold for whether or not phase discontinuity is possible + + Output: + unw_igram :: Raster, Unwrapped interferogram, type: isce3.io.Raster + conn_comp :: Raster, Connected component labels, type: isce3.io.Raster + ''' + + from isce3.io.gdal import Raster # Snaphu module takes isce3.io.gdal.Raster inputs + from isce3.unwrap import snaphu + + # Input file + igram = Raster(infile) + correlation = Raster(corr_file) + + if mask: + print('Using mask file {}'.format(mask)) + mask = Raster(mask) + + # Output file + unw_igram = Raster(unw_filename, igram.width, igram.length, isce3.io.gdal.GDT_Float32, "ENVI") #works also with ISCE format + conn_comp = Raster(conncomp_filename, igram.width, igram.length, isce3.io.gdal.GDT_UInt32, "ENVI") + + # Snaphu parameters + tiling_params = snaphu.TilingParams(nproc=nproc, tile_nrows=tiles[0], tile_ncols=tiles[1], + row_overlap=tiles[2], col_overlap=tiles[3], tile_cost_thresh=tile_cost_thresh, + min_region_size=min_region, single_tile_reoptimize=single_tile_opt) + + conncomp_params = snaphu.ConnCompParams() + corr_bias_model_params = snaphu.CorrBiasModelParams(min_corr_factor = defo_thresh) + + cost_params = None + solver_params = None + init_method = 'mcf' #not yet included in the isce3.unwrap.snaphu + pwr = None + unwest = None + + phase_stddev_model_params = None + scratchdir = ./ + verbose = print_msg + debug = False + + if print_msg: + print('\nSnaphu Unwrapping:') + print('Input:') + print(' Phase: {}'.format(os.path.abspath(infile))) + print(' width: {}, length {}, dtype: {}'.format(igram.width, igram.length, GDAL_DATATYPE[igram.datatype.value])) + print(' Correlation: {}'.format(os.path.abspath(corr_file))) + print(' width: {}, length {}, dtype: {}'.format(correlation.width, correlation.length, GDAL_DATATYPE[correlation.datatype.value])) + print(' Parameters: nlooks={}, tiles={}, tile_cost_thresh={}, min_region={}, nproc={}, single_tile_opt={}, defo_thresh={}'.format(nlooks, tiles, tile_cost_thresh, min_region, nproc, single_tile_opt, defo_thresh)) + + # Unwrap + snaphu.unwrap(unw_igram, conn_comp, igram, correlation, nlooks, + "defo", cost_params, pwr, mask, unwest, tiling_params, solver_params,conncomp_params, + corr_bias_model_params, phase_stddev_model_params, scratchdir, verbose, debug) + + if print_msg: + print('#############################################') + print('Finished unwrapping, output files:') + print(' Unwrapped interferogram: {}'.format(os.path.abspath(unw_filename))) + print(' width: {}, length {}, dtype: {}'.format(unw_igram.width, unw_igram.length, GDAL_DATATYPE[unw_igram.datatype.value])) + print(' Connected components: {}'.format(os.path.abspath(conncomp_filename))) + print(' width: {}, length {}, dtype: {}'.format(conn_comp.width, conn_comp.length, GDAL_DATATYPE[conn_comp.datatype.value])) + + +#################################################################################### +#################################################################################### +if __name__ == '__main__': + + ''' + Main driver + ''' + + inps = cmdLineParser() + + ######################### Prepare input and output files ####################### + input_name = os.path.splitext(os.path.basename(inps.igram))[0] + unw_output = os.path.join(inps.outputDir, input_name + '.unw') + conncomp_output = os.path.join(inps.outputDir, input_name + '_int.conn') + + # Need gdal vrt file to read correctly the tcorr_ds_ps.bin format + corr = Raster(inps.correlation) + write_vrt(inps.correlation, corr.width, corr.length, corr.datatype.value) + inps.correlation = os.path.abspath(inps.correlation + '.vrt') + + ################################ Set-up and Unwrap ############################# + if inps.unwrapper == 'phass': + corr_thresh = 0.3 + good_corr = 0.7 + min_region = 200 + + if inps.kwargs: + if 'corr_thresh' in inps.kwargs: + corr_thresh = float(inps.kwargs['corr_thresh']) + if 'good_corr' in inps.kwargs: + good_corr = float(inps.kwargs['good_corr']) + if 'min_region' in inps.kwargs: + min_region = int(inps.kwargs['min_region']) + + unwrap_phass(inps.igram + '.vrt', inps.correlation, unw_output, conncomp_output, + correlation_threshold=corr_thresh, good_correlation=good_corr, min_region=min_region) + + elif inps.unwrapper == 'snaphu': + nlooks = 27 + tiles = [1,1,0,0] + min_region = 200 + nproc = 1 + single_tile_opt = False + + if inps.kwargs: + if 'nlooks' in inps.kwargs: + nlooks = int(inps.kwargs['nlooks']) + if 'tiles' in inps.kwargs: + tile = inps.kwargs['tiles'].strip('][').split(',') + tiles = [int(i) for i in tile] + if 'min_region' in inps.kwargs: + min_region = int(inps.kwargs['min_region']) + if 'nproc' in inps.kwargs: + nproc = int(inps.kwargs['nproc']) + if 'single_tile_opt' in inps.kwargs: + single_tile_opt = int(inps.kwargs['single_tile_opt']) + + unwrap_snaphu(inps.igram + '.vrt', inps.correlation, unw_output, conncomp_output, mask = inps.mask, + nlooks=nlooks, tiles=tiles, min_region=min_region, nproc=nproc, single_tile_opt=single_tile_opt) + + ################################ Write output vrt and xml file ############################# + + unw_vrt_xml(unw_output, "BIL") + unw_vrt_xml(conncomp_output, "BIP") + + + + + + From 79502fe9adc9e3098052c07403902b5049d98d2b Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 11 Apr 2022 17:33:20 -0700 Subject: [PATCH 2/7] Update unwrap_isce3.py --- python/unwrap_isce3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/unwrap_isce3.py b/python/unwrap_isce3.py index ac60778..c1399e5 100644 --- a/python/unwrap_isce3.py +++ b/python/unwrap_isce3.py @@ -338,7 +338,7 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, unwest = None phase_stddev_model_params = None - scratchdir = ./ + scratchdir = './' verbose = print_msg debug = False From 211b2cf724e38ab56e3d5750c042641a6cd39ef6 Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 11 Apr 2022 17:58:28 -0700 Subject: [PATCH 3/7] Fixed most of codacy issues left issues with xml.Elementtree unsolved --- python/unwrap_isce3.py | 41 ++++++++++++++++++----------------------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/python/unwrap_isce3.py b/python/unwrap_isce3.py index c1399e5..28bb58f 100644 --- a/python/unwrap_isce3.py +++ b/python/unwrap_isce3.py @@ -6,7 +6,6 @@ import os import argparse import isce3 -import numpy as np from isce3.io.gdal import Raster import xml.etree.ElementTree as ET from xml.dom import minidom @@ -19,8 +18,7 @@ def __call__(self, parser, namespace, values, option_string=None): getattr(namespace, self.dest)[key] = value def cmdLineParser(): - ''' - Command Line Parser + ''' Command Line Parser Script to use Phass/Snaphu in the isce3 environment USAGE: @@ -29,10 +27,10 @@ def cmdLineParser(): unwrap_isce3.py -i 20190610_20190716.int -c tcorr_ds_ps.bin -k corr_thresh=0.3, good_corr=0.7 min_region=200 Snaphu: - unwrap_isce3.py -i 20190610_20190716.int -c tcorr_ds_ps.bin -u snaphu -k nlooks=27 tiles=[4,4,500,500] nproc=8 min_region=200 + unwrap_isce3.py -i 20190610_20190716.int -c tcorr_ds_ps.bin -u snaphu -k nlooks=27 tiles=[4,4,500,500] nproc=8 min_region=200 ''' - unw_options = ''' phass :: corr_thresh=0.3, good_corr=0.7, min_region=200, \n - snaphu :: nlooks=1, tiles=[1, 1, 0, 0], min_region=200, nproc=1, single_tile_opt=False ''' + unw_options = ''' phass :: corr_thresh=0.3, good_corr=0.7, min_region=200, \n + snaphu :: nlooks=1, tiles=[1, 1, 0, 0], min_region=200, nproc=1, single_tile_opt=False ''' parser = argparse.ArgumentParser(description='Phass Unwrapper', formatter_class=argparse.RawTextHelpFormatter) parser.add_argument('-i', '--interferogram_file', type=str, dest='igram', @@ -41,9 +39,9 @@ def cmdLineParser(): required=True, help='Input correlation file, example: tcorr_ds_ps.bin') parser.add_argument('-m', '--mask', type=str, dest='mask', default=None, help='Input mask byte file') - parser.add_argument('-o', '--output_dir', type=str, dest='outputDir', + parser.add_argument('-o', '--output_dir', type=str, dest='outputDir', default='./', help='Output directory path') - parser.add_argument('-u', '--unwrapper', type=str, dest='unwrapper', + parser.add_argument('-u', '--unwrapper', type=str, dest='unwrapper', default='phass', help='Unwrapper method: phass/snaphu') parser.add_argument('-k', '--kwargs', type=str, dest='kwargs', nargs='*', help=unw_options, action=ParseKwargs) @@ -121,7 +119,6 @@ def component_xml(root, name, size, doc): property_xml(coord, 'startingvalue', '0.0', 'Starting value of the coordinate.') def write_xml(infile, width, length, bands, datatype, scheme): - root = ET.Element('imageFile') property_xml(root, 'access_mode', 'READ', 'Image access mode.') property_xml(root, 'byte_order', 'l', 'Endianness of the image.') @@ -163,7 +160,7 @@ def vrt_cpx2rad(infile, width, length, stype): Function to extract interferogram phase from complex file needed when unwrapping interferogram in radians. - Find the input interferogram dtype, if complex32 or complex64 + Find the input interferogram dtype, if complex32 or complex64 reate vrt file that extract phase values onthefly Input: @@ -215,9 +212,9 @@ def vrt_cpx2rad(infile, width, length, stype): # Return the path to VRT file as input to Phass return outfile -def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, +def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, correlation_threshold=0.3, good_correlation=0.7, min_region=200, print_msg=True): - ''' + ''' Unwrap wrapped interferogram [in radians] using Phass Input: @@ -239,14 +236,14 @@ def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, phass = Phass() - # If input interferogram is in COMPLEX format, create .vrt file + # If input interferogram is in COMPLEX format, create .vrt file # to get interferogram phase in radians on-the-fly infile_name = infile infile = Raster(infile) phass_input = vrt_cpx2rad(infile_name, infile.width, infile.length, infile.datatype()) # Input files - igram = Raster(phass_input) + igram = Raster(phass_input) correlation = Raster(corr_file) # Output files @@ -261,7 +258,7 @@ def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, if print_msg: print('\nPhass Unwrapping:') print('Input:') - print(' Phase: {}'.format(os.path.abspath(phase_input))) + print(' Phase: {}'.format(os.path.abspath(phass_input))) print(' width: {}, length {}, dtype: {}'.format(igram.width, igram.length, GDAL_DATATYPE[igram.datatype()])) print(' Correlation: {}'.format(os.path.abspath(corr_file))) print(' width: {}, length {}, dtype: {}'.format(correlation.width, correlation.length, GDAL_DATATYPE[correlation.datatype()])) @@ -296,7 +293,7 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, tiles :: list [tile_nrows, tile_ncols, row_overlap, col_overlap] tile_nrows, tile_ncols - Number of tiles along the row/column directions row_overlap, col_overlap - verlap, in number of rows/columns, between neighboring tiles - tile_cost_thresh :: int Cost threshold to use for determining boundaries of reliable regions. + tile_cost_thresh :: int Cost threshold to use for determining boundaries of reliable regions. Larger cost threshold implies smaller regions (safer, but more expensive computationally min_region :: int Minimum number of pixels per region (connected component) nproc :: int Number of threads for parallel processing (use with tiles) @@ -312,7 +309,7 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, from isce3.unwrap import snaphu # Input file - igram = Raster(infile) + igram = Raster(infile) correlation = Raster(corr_file) if mask: @@ -333,7 +330,7 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, cost_params = None solver_params = None - init_method = 'mcf' #not yet included in the isce3.unwrap.snaphu + #init_method = 'mcf' #not yet included in the isce3.unwrap.snaphu pwr = None unwest = None @@ -369,9 +366,7 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, #################################################################################### if __name__ == '__main__': - ''' - Main driver - ''' + # Main driver inps = cmdLineParser() @@ -403,10 +398,10 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, correlation_threshold=corr_thresh, good_correlation=good_corr, min_region=min_region) elif inps.unwrapper == 'snaphu': - nlooks = 27 + nlooks = 27 tiles = [1,1,0,0] min_region = 200 - nproc = 1 + nproc = 1 single_tile_opt = False if inps.kwargs: From ca93c29b7befc32e6faf669ce46168de0162e406 Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 11 Apr 2022 18:05:41 -0700 Subject: [PATCH 4/7] Update unwrap_isce3.py --- python/unwrap_isce3.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/python/unwrap_isce3.py b/python/unwrap_isce3.py index 28bb58f..3647094 100644 --- a/python/unwrap_isce3.py +++ b/python/unwrap_isce3.py @@ -18,7 +18,8 @@ def __call__(self, parser, namespace, values, option_string=None): getattr(namespace, self.dest)[key] = value def cmdLineParser(): - ''' Command Line Parser + ''' + Command Line Parser Script to use Phass/Snaphu in the isce3 environment USAGE: @@ -78,10 +79,10 @@ def cmdLineParser(): } def write_vrt(infile, width, length, dtype): - ## infile :: str path to inputfile + ## infile :: str path to inputfile infile_dir = os.path.dirname(os.path.abspath(infile)) - + vrttmpl=''' {PATH} From e839b0266c5f01a19274948c1ffdf074c225c982 Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 11 Apr 2022 18:28:44 -0700 Subject: [PATCH 5/7] Fix Codacy issues ver2 --- python/unwrap_isce3.py | 67 +++++++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 34 deletions(-) diff --git a/python/unwrap_isce3.py b/python/unwrap_isce3.py index 3647094..e073eee 100644 --- a/python/unwrap_isce3.py +++ b/python/unwrap_isce3.py @@ -5,10 +5,13 @@ import os import argparse -import isce3 -from isce3.io.gdal import Raster import xml.etree.ElementTree as ET from xml.dom import minidom +import isce3 +from isce3.unwrap import snaphu +from isce3.unwrap import Phass +from isce3.io.gdal import Raster as gRaster # Snaphu module takes isce3.io.gdal.Raster inputs +from isce3.io import Raster # Phass module takes isce3.io.Raster inputs class ParseKwargs(argparse.Action): def __call__(self, parser, namespace, values, option_string=None): @@ -18,6 +21,7 @@ def __call__(self, parser, namespace, values, option_string=None): getattr(namespace, self.dest)[key] = value def cmdLineParser(): + ''' Command Line Parser Script to use Phass/Snaphu in the isce3 environment @@ -80,7 +84,6 @@ def cmdLineParser(): def write_vrt(infile, width, length, dtype): ## infile :: str path to inputfile - infile_dir = os.path.dirname(os.path.abspath(infile)) vrttmpl=''' @@ -116,7 +119,7 @@ def component_xml(root, name, size, doc): property_xml(coord, 'endingvalues', str(float(size)), 'Ending value of the coordinate.') property_xml(coord, 'family', 'imagecoordinate', 'Instance family name.') property_xml(coord, 'name', 'imagecoordinate_name', 'Instance name.') - property_xml(coord, 'size', str(width), 'Coordinate size.') + property_xml(coord, 'size', str(size), 'Coordinate size.') property_xml(coord, 'startingvalue', '0.0', 'Starting value of the coordinate.') def write_xml(infile, width, length, bands, datatype, scheme): @@ -125,7 +128,7 @@ def write_xml(infile, width, length, bands, datatype, scheme): property_xml(root, 'byte_order', 'l', 'Endianness of the image.') component_xml(root, 'coordinate1', width, 'First coordinate of a 2D image (width).') component_xml(root, 'coordinate2', length, 'Second coordinate of a 2D image (length).') - property_xml(root, 'data_type', GDAL2ISCE[dtype], 'Image data type.') + property_xml(root, 'data_type', GDAL2ISCE[datatype], 'Image data type.') property_xml(root, 'family', 'image', 'Instance family name.') property_xml(root, 'file_name', os.path.abspath(infile), 'Name of the image file.') property_xml(root, 'length', str(length), 'Image length.') @@ -231,10 +234,7 @@ def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, unw_igram :: Raster, Unwrapped interferogram, type: isce3.io.Raster conn_comp :: Raster, Connected component labels, type: isce3.io.Raster ''' - - from isce3.io import Raster # Phass module takes isce3.io.Raster inputs - from isce3.unwrap import Phass - + phass = Phass() # If input interferogram is in COMPLEX format, create .vrt file @@ -276,9 +276,9 @@ def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, print(' Connected components: {}'.format(os.path.abspath(conncomp_filename))) print(' width: {}, length {}, dtype: {}'.format(conn_comp.width, conn_comp.length, GDAL_DATATYPE[conn_comp.datatype()])) -def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, - cost_mode='defo', nlooks=1, tiles=[1,1,0,0], tile_cost_thresh=500, min_region=100, - nproc=1, single_tile_opt=False, defo_thresh=1.25, print_msg=True): +def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, cost_mode='defo', + nlooks=1, tile_nrows=1, tile_ncols=1, row_overlap=0, col_overlap=0, tile_cost_thresh=500, + min_region=100,nproc=1, single_tile_opt=False, defo_thresh=1.25, print_msg=True): ''' Unwrap wrapped interferogram using Snaphu @@ -291,9 +291,10 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, mask :: str Binary mask of valid pixels. Zeros in this raster indicate interferogram pixels that should be masked out, GDT_Byte datatype nlooks :: float Effective number of looks used to form the input correlation data - tiles :: list [tile_nrows, tile_ncols, row_overlap, col_overlap] - tile_nrows, tile_ncols - Number of tiles along the row/column directions - row_overlap, col_overlap - verlap, in number of rows/columns, between neighboring tiles + tile_nrows :: int Number of tiles along the row directions + tile_ncols :: int Number of tiles along the column directions + row_overlap :: int Overlap in number of rows between neighboring tiles + col_overlap :: int Overlap in number of columns between neighboring tiles tile_cost_thresh :: int Cost threshold to use for determining boundaries of reliable regions. Larger cost threshold implies smaller regions (safer, but more expensive computationally min_region :: int Minimum number of pixels per region (connected component) @@ -306,24 +307,21 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, conn_comp :: Raster, Connected component labels, type: isce3.io.Raster ''' - from isce3.io.gdal import Raster # Snaphu module takes isce3.io.gdal.Raster inputs - from isce3.unwrap import snaphu - # Input file - igram = Raster(infile) - correlation = Raster(corr_file) + igram = gRaster(infile) + correlation = gRaster(corr_file) if mask: print('Using mask file {}'.format(mask)) - mask = Raster(mask) + mask = gRaster(mask) # Output file - unw_igram = Raster(unw_filename, igram.width, igram.length, isce3.io.gdal.GDT_Float32, "ENVI") #works also with ISCE format - conn_comp = Raster(conncomp_filename, igram.width, igram.length, isce3.io.gdal.GDT_UInt32, "ENVI") + unw_igram = gRaster(unw_filename, igram.width, igram.length, isce3.io.gdal.GDT_Float32, "ENVI") #works also with ISCE format + conn_comp = gRaster(conncomp_filename, igram.width, igram.length, isce3.io.gdal.GDT_UInt32, "ENVI") # Snaphu parameters - tiling_params = snaphu.TilingParams(nproc=nproc, tile_nrows=tiles[0], tile_ncols=tiles[1], - row_overlap=tiles[2], col_overlap=tiles[3], tile_cost_thresh=tile_cost_thresh, + tiling_params = snaphu.TilingParams(nproc=nproc, tile_nrows=tile_nrows, tile_ncols=tile_ncols, + row_overlap=row_overlap, col_overlap=col_overlap, tile_cost_thresh=tile_cost_thresh, min_region_size=min_region, single_tile_reoptimize=single_tile_opt) conncomp_params = snaphu.ConnCompParams() @@ -377,7 +375,7 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, conncomp_output = os.path.join(inps.outputDir, input_name + '_int.conn') # Need gdal vrt file to read correctly the tcorr_ds_ps.bin format - corr = Raster(inps.correlation) + corr = gRaster(inps.correlation) write_vrt(inps.correlation, corr.width, corr.length, corr.datatype.value) inps.correlation = os.path.abspath(inps.correlation + '.vrt') @@ -404,6 +402,8 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, min_region = 200 nproc = 1 single_tile_opt = False + defo_thresh = 1.25 + tile_cost_thresh = 500 if inps.kwargs: if 'nlooks' in inps.kwargs: @@ -417,17 +417,16 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, nproc = int(inps.kwargs['nproc']) if 'single_tile_opt' in inps.kwargs: single_tile_opt = int(inps.kwargs['single_tile_opt']) + if 'defo_thresh' in inps.kwargs: + defo_thresh = float(inps.kwargs['defo_thresh']) + if 'tile_cost_thresh' in inps.kwargs: + tile_cost_thresh = float(inps.kwargs['tile_cost_thresh']) - unwrap_snaphu(inps.igram + '.vrt', inps.correlation, unw_output, conncomp_output, mask = inps.mask, - nlooks=nlooks, tiles=tiles, min_region=min_region, nproc=nproc, single_tile_opt=single_tile_opt) + unwrap_snaphu(inps.igram + '.vrt', inps.correlation, unw_output, conncomp_output, mask = inps.mask, nlooks=nlooks, + tile_nrows=tiles[0], tile_ncols=tiles[1], row_overlap=tiles[2], col_overlap=tiles[3], min_region=min_region, + nproc=nproc, single_tile_opt=single_tile_opt, defo_thresh=defo_thresh, tile_cost_thresh=tile_cost_thresh)) ################################ Write output vrt and xml file ############################# unw_vrt_xml(unw_output, "BIL") unw_vrt_xml(conncomp_output, "BIP") - - - - - - From 72027eda346ce07fdecf0dc7f21c9ae32223df1a Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 11 Apr 2022 18:32:01 -0700 Subject: [PATCH 6/7] Update unwrap_isce3.py --- python/unwrap_isce3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/unwrap_isce3.py b/python/unwrap_isce3.py index e073eee..04ec992 100644 --- a/python/unwrap_isce3.py +++ b/python/unwrap_isce3.py @@ -424,7 +424,7 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, unwrap_snaphu(inps.igram + '.vrt', inps.correlation, unw_output, conncomp_output, mask = inps.mask, nlooks=nlooks, tile_nrows=tiles[0], tile_ncols=tiles[1], row_overlap=tiles[2], col_overlap=tiles[3], min_region=min_region, - nproc=nproc, single_tile_opt=single_tile_opt, defo_thresh=defo_thresh, tile_cost_thresh=tile_cost_thresh)) + nproc=nproc, single_tile_opt=single_tile_opt, defo_thresh=defo_thresh, tile_cost_thresh=tile_cost_thresh) ################################ Write output vrt and xml file ############################# From 0db5ec4d5ebe7c9b29e1f6f773f1237fb85f35ff Mon Sep 17 00:00:00 2001 From: Marin Govorcin Date: Mon, 11 Apr 2022 20:41:35 -0700 Subject: [PATCH 7/7] Update unwrap_isce3.py --- python/unwrap_isce3.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/python/unwrap_isce3.py b/python/unwrap_isce3.py index 04ec992..8a4b700 100644 --- a/python/unwrap_isce3.py +++ b/python/unwrap_isce3.py @@ -8,10 +8,9 @@ import xml.etree.ElementTree as ET from xml.dom import minidom import isce3 -from isce3.unwrap import snaphu -from isce3.unwrap import Phass -from isce3.io.gdal import Raster as gRaster # Snaphu module takes isce3.io.gdal.Raster inputs +from isce3.unwrap import snaphu, Phass from isce3.io import Raster # Phass module takes isce3.io.Raster inputs +from isce3.io.gdal import Raster as gRaster # Snaphu module takes isce3.io.gdal.Raster inputs class ParseKwargs(argparse.Action): def __call__(self, parser, namespace, values, option_string=None): @@ -21,8 +20,8 @@ def __call__(self, parser, namespace, values, option_string=None): getattr(namespace, self.dest)[key] = value def cmdLineParser(): - ''' + Command Line Parser Script to use Phass/Snaphu in the isce3 environment @@ -278,7 +277,7 @@ def unwrap_phass(infile, corr_file, unw_filename, conncomp_filename, def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, cost_mode='defo', nlooks=1, tile_nrows=1, tile_ncols=1, row_overlap=0, col_overlap=0, tile_cost_thresh=500, - min_region=100,nproc=1, single_tile_opt=False, defo_thresh=1.25, print_msg=True): + min_region=100, nproc=1, single_tile_opt=False, defo_thresh=1.25, print_msg=True): ''' Unwrap wrapped interferogram using Snaphu @@ -379,6 +378,10 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, write_vrt(inps.correlation, corr.width, corr.length, corr.datatype.value) inps.correlation = os.path.abspath(inps.correlation + '.vrt') + # Check if output folder exists + if not os.path.exists(inps.outputDir): + os.makedirs(inps.outputDir) + ################################ Set-up and Unwrap ############################# if inps.unwrapper == 'phass': corr_thresh = 0.3 @@ -422,8 +425,8 @@ def unwrap_snaphu(infile, corr_file, unw_filename, conncomp_filename, mask=None, if 'tile_cost_thresh' in inps.kwargs: tile_cost_thresh = float(inps.kwargs['tile_cost_thresh']) - unwrap_snaphu(inps.igram + '.vrt', inps.correlation, unw_output, conncomp_output, mask = inps.mask, nlooks=nlooks, - tile_nrows=tiles[0], tile_ncols=tiles[1], row_overlap=tiles[2], col_overlap=tiles[3], min_region=min_region, + unwrap_snaphu(inps.igram + '.vrt', inps.correlation, unw_output, conncomp_output, mask = inps.mask, nlooks=nlooks, + tile_nrows=tiles[0], tile_ncols=tiles[1], row_overlap=tiles[2], col_overlap=tiles[3], min_region=min_region, nproc=nproc, single_tile_opt=single_tile_opt, defo_thresh=defo_thresh, tile_cost_thresh=tile_cost_thresh) ################################ Write output vrt and xml file #############################