diff --git a/contrib/stack/topsStack/unwrap.py b/contrib/stack/topsStack/unwrap.py index c297dc8..fd3931e 100755 --- a/contrib/stack/topsStack/unwrap.py +++ b/contrib/stack/topsStack/unwrap.py @@ -38,9 +38,11 @@ from isceobj.Constants import SPEED_OF_LIGHT import argparse import os import pickle -import sys +import gdal +import numpy as np #import shelve import s1a_isce_utils as ut +from isceobj.Util.ImageUtil import ImageLib as IML def createParser(): ''' @@ -70,6 +72,9 @@ def createParser(): parser.add_argument('-m', '--method', dest='method', type=str, default='icu', help='unwrapping method') + parser.add_argument('--rmfilter', action='store_true', default=False, + help='remove the effect of filtering from final unwrapped interferograms') + return parser @@ -276,6 +281,53 @@ def runUnwrapIcu(infile, outfile): #unwImage.finalizeImage() unwImage.renderHdr() + +def remove_filter(intfile, filtfile, unwfile): + + outunw = os.path.abspath(unwfile).split('filt_') + outunw = outunw[0] + outunw[1] + + ds_unw = gdal.Open(unwfile + ".vrt", gdal.GA_ReadOnly) + unw_phase = ds_unw.ReadAsArray() + + bands, length, width = unw_phase.shape + print(bands, length, width) + outfile = IML.memmap(outunw, mode='write', nchannels=2, + nxx=width, nyy=length, scheme='BIL', dataType='f') + + outfile.bands[0][:, :] = unw_phase[0, :, :] + + ds_filt = gdal.Open(filtfile + ".vrt", gdal.GA_ReadOnly) + filt_phase = np.angle(ds_filt.ReadAsArray()) + + integer_jumps = unw_phase[1, :, :] - filt_phase + + del unw_phase, filt_phase + + ds_int = gdal.Open(intfile + ".vrt", gdal.GA_ReadOnly) + int_phase = np.angle(ds_int.ReadAsArray()) + + out_phase = integer_jumps + int_phase + + del int_phase, integer_jumps + + outfile.bands[1][:, :] = out_phase[:, :] + + unwImage = isceobj.Image.createImage() + unwImage.setFilename(outunw) + unwImage.setWidth(width) + unwImage.setLength(length) + unwImage.imageType = 'unw' + unwImage.bands = 2 + unwImage.scheme = 'BIL' + unwImage.dataType = 'FLOAT' + unwImage.setAccessMode('read') + #unwImage.createImage() + unwImage.renderHdr() + #unwImage.finalizeImage() + + return + def main(iargs=None): ''' The main driver. @@ -284,6 +336,7 @@ def main(iargs=None): inps = cmdLineParse(iargs) interferogramDir = os.path.dirname(inps.intfile) print ('unwrapping method : ' , inps.method) + ''' if inps.method == 'snaphu': if inps.nomcf: fncall = runUnwrap @@ -297,6 +350,14 @@ def main(iargs=None): elif inps.method == 'icu': runUnwrapIcu(inps.intfile, inps.unwfile) + ''' + + if inps.rmfilter: + filtfile = os.path.abspath(inps.intfile) + intfile = filtfile.split('filt_') + intfile = intfile[0] + intfile[1] + + remove_filter(intfile, filtfile, inps.unwfile) if __name__ == '__main__':