From 90d578fa5d8d27c3ed6e9ffd4bb93753c66defa3 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Fri, 29 May 2020 13:37:07 -0700 Subject: [PATCH] cuDenseOffsets: expose --raw-osf and --corr-win-size option contrib/PyCuAmpcor/cuDenseOffsets.py: + expose --raw-over-samp-factor option to change the raw data oversampling factor with default of 2 + expose --corr-win-size option, which by default, is calculated automatically based on the search window size and raw oversampling factor + change the namespace key for the oversampling factor of correlation surface from 'oversample' to 'corr_oversample' to avoid ambiguity, the cmd option name is still the same, thus no different in the user side. + use vrt file for master / slave slc data loading folllowing https://github.com/isce-framework/isce2/pull/78 + re-organize createParse() for readability + add example to argparse help stripmapStack/MaskAndFilter.py: + fix a bug introduced in #https://github.com/isce-framework/isce2/pull/112 + add example usage to help message + mask pixels with SNR == 0, for plotting only without touching data files + show input parameter in the figure title + show percentage after SNR thresholding + add interpolation for matplotlib while plotting topsStack/grossOffsets: move obsolete basemap import from the top to the related commented code block --- contrib/PyCuAmpcor/examples/cuDenseOffsets.py | 204 ++++++++++-------- contrib/stack/stripmapStack/MaskAndFilter.py | 76 +++++-- contrib/stack/topsStack/grossOffsets.py | 3 +- 3 files changed, 174 insertions(+), 109 deletions(-) diff --git a/contrib/PyCuAmpcor/examples/cuDenseOffsets.py b/contrib/PyCuAmpcor/examples/cuDenseOffsets.py index 1495047..a77edb1 100755 --- a/contrib/PyCuAmpcor/examples/cuDenseOffsets.py +++ b/contrib/PyCuAmpcor/examples/cuDenseOffsets.py @@ -1,109 +1,127 @@ #!/usr/bin/env python3 -# Author: Minyan Zhong, Lijun Zhu +# Author: Minyan Zhong, Lijun Zhu + -import argparse import os +import argparse +import numpy as np + import isce import isceobj from isceobj.Util.decorators import use_api - -import numpy as np from contrib.PyCuAmpcor.PyCuAmpcor import PyCuAmpcor + +EXAMPLE = '''example + cuDenseOffsets.py -m ./merged/SLC/20151120/20151120.slc.full -s ./merged/SLC/20151214/20151214.slc.full + --masterxml ./master/IW1.xml --outprefix ./merged/offsets/20151120_20151214/offset + --ww 256 --wh 256 --oo 32 --kw 300 --kh 100 --nwac 100 --nwdc 1 --sw 8 --sh 8 --gpuid 2 +''' + + def createParser(): ''' Command line parser. ''' - parser = argparse.ArgumentParser( description='Generate offset field between two Sentinel slc') + parser = argparse.ArgumentParser(description='Generate offset field between two Sentinel slc', + formatter_class=argparse.RawTextHelpFormatter, + epilog=EXAMPLE) parser.add_argument('-m','--master', type=str, dest='master', required=True, - help='Master image') + help='Master image') parser.add_argument('-s', '--slave',type=str, dest='slave', required=True, - help='Slave image') + help='Slave image') parser.add_argument('-l', '--lat',type=str, dest='lat', required=False, - help='Latitude') + help='Latitude') parser.add_argument('-L', '--lon',type=str, dest='lon', required=False, - help='Longitude') + help='Longitude') parser.add_argument('--los',type=str, dest='los', required=False, - help='Line of Sight') - parser.add_argument('--masterxml',type=str, dest='masterxml', required=False, - help='Master Image Xml File') + help='Line of Sight') + parser.add_argument('-x', '--masterxml',type=str, dest='masterxml', required=False, + help='Master Image XML File') + + parser.add_argument('--op','--outprefix','--output-prefix', type=str, dest='outprefix', + default='offset', required=True, + help='Output prefix, default: offset.') + parser.add_argument('--os','--outsuffix', type=str, dest='outsuffix', default='', + help='Output suffix, default:.') + parser.add_argument('--ww', type=int, dest='winwidth', default=64, - help='Window Width') + help='Window width (default: %(default)s).') parser.add_argument('--wh', type=int, dest='winhgt', default=64, - help='Window height') - parser.add_argument('--sw', type=int, dest='srcwidth', default=20, - help='Search window width') - parser.add_argument('--sh', type=int, dest='srchgt', default=20, - help='Search window height') + help='Window height (default: %(default)s).') + + parser.add_argument('--sw', type=int, dest='srcwidth', default=20, choices=range(8, 33), + help='Search window width (default: %(default)s).') + parser.add_argument('--sh', type=int, dest='srchgt', default=20, choices=range(8, 33), + help='Search window height (default: %(default)s).') parser.add_argument('--mm', type=int, dest='margin', default=50, - help='Margin') + help='Margin (default: %(default)s).') + parser.add_argument('--kw', type=int, dest='skipwidth', default=64, - help='Skip across') + help='Skip across (default: %(default)s).') parser.add_argument('--kh', type=int, dest='skiphgt', default=64, - help='Skip down') + help='Skip down (default: %(default)s).') + + parser.add_argument('--raw-osf','--raw-over-samp-factor', type=int, dest='raw_oversample', + default=2, choices=range(2,5), + help='raw data oversampling factor (default: %(default)s).') + + gross = parser.add_argument_group('Initial gross offset') + gross.add_argument('-g','--gross', type=int, dest='gross', default=0, + help='Use gross offset or not') + gross.add_argument('--aa', type=int, dest='azshift', default=0, + help='Gross azimuth offset (default: %(default)s).') + gross.add_argument('--rr', type=int, dest='rgshift', default=0, + help='Gross range offset (default: %(default)s).') + + corr = parser.add_argument_group('Correlation surface') + corr.add_argument('--corr-win-size', type=int, dest='corr_win_size', default=-1, + help='Zoom-in window size of the correlation surface for oversampling (default: %(default)s).') + corr.add_argument('--corr-osf', '--oo', '--corr-over-samp-factor', type=int, dest='corr_oversample', default=32, + help = 'Oversampling factor of the zoom-in correlation surface (default: %(default)s).') parser.add_argument('--nwa', type=int, dest='numWinAcross', default=-1, - help='Number of Window Across') + help='Number of window across (default: %(default)s).') parser.add_argument('--nwd', type=int, dest='numWinDown', default=-1, - help='Number of Window Down') + help='Number of window down (default: %(default)s).') parser.add_argument('--nwac', type=int, dest='numWinAcrossInChunk', default=1, - help='Number of Window Across in Chunk') + help='Number of window across in chunk (default: %(default)s).') parser.add_argument('--nwdc', type=int, dest='numWinDownInChunk', default=1, - help='Number of Window Down in Chunk') + help='Number of window down in chunk (default: %(default)s).') + parser.add_argument('-r', '--redo', dest='redo', action='store_true', + help='To redo by force (ignore the existing offset fields).') - parser.add_argument('-op','--outprefix', type=str, dest='outprefix', default='dense_ampcor', required=True, - help='Output prefix') + parser.add_argument('--drmp', '--deramp', dest='deramp', type=int, default=0, + help='deramp method (0: mag, 1: complex) (default: %(default)s).') - parser.add_argument('-os','--outsuffix', type=str, dest='outsuffix',default='dense_ampcor', - help='Output suffix') - - parser.add_argument('-g','--gross', type=int, dest='gross', default=0, - help='Use gross offset or not') - - parser.add_argument('--aa', type=int, dest='azshift', default=0, - help='Gross azimuth offset') - - parser.add_argument('--rr', type=int, dest='rgshift', default=0, - help='Gross range offset') - - parser.add_argument('--oo', type=int, dest='oversample', default=32, - help = 'Oversampling factor') - - parser.add_argument('-r', '--redo', dest='redo', type=int, default=0 - , help='To redo or not') - - parser.add_argument('-drmp', '--deramp', dest='deramp', type=int, default=0 - , help='deramp method (0: mag, 1: complex)') - - parser.add_argument('-gid', '--gpuid', dest='gpuid', type=int, default=-1 - , help='GPU ID') + parser.add_argument('--gpuid', '--gid', '--gpu-id', dest='gpuid', type=int, default=-1, + help='GPU ID (default: %(default)s).') return parser + def cmdLineParse(iargs = None): parser = createParser() inps = parser.parse_args(args=iargs) return inps + @use_api def estimateOffsetField(master, slave, inps=None): - import pathlib - ###Loading the slave image object sim = isceobj.createSlcImage() - sim.load(pathlib.Path(slave).with_suffix('.xml')) + sim.load(slave+'.xml') sim.setAccessMode('READ') sim.createImage() - ###Loading the master image object sar = isceobj.createSlcImage() - sar.load(pathlib.Path(master).with_suffix('.xml')) + sar.load(master+'.xml') sar.setAccessMode('READ') sar.createImage() @@ -111,19 +129,19 @@ def estimateOffsetField(master, slave, inps=None): length = sar.getLength() objOffset = PyCuAmpcor() - + objOffset.algorithm = 0 objOffset.deviceID = inps.gpuid # -1:let system find the best GPU - objOffset.nStreams = 1 #cudaStreams + objOffset.nStreams = 2 #cudaStreams objOffset.derampMethod = inps.deramp + print('deramp method (0 for magnitude, 1 for complex): ', objOffset.derampMethod) - objOffset.masterImageName = master + objOffset.masterImageName = master+'.vrt' objOffset.masterImageHeight = length objOffset.masterImageWidth = width - objOffset.slaveImageName = slave + objOffset.slaveImageName = slave+'.vrt' objOffset.slaveImageHeight = length objOffset.slaveImageWidth = width - print("image length:",length) print("image width:",width) @@ -132,39 +150,50 @@ def estimateOffsetField(master, slave, inps=None): if (inps.numWinDown != -1): objOffset.numberWindowDown = inps.numWinDown - if (inps.numWinAcross != -1): objOffset.numberWindowAcross = inps.numWinAcross - print("offset field length: ",objOffset.numberWindowDown) print("offset field width: ",objOffset.numberWindowAcross) # window size objOffset.windowSizeHeight = inps.winhgt objOffset.windowSizeWidth = inps.winwidth - + print('cross correlation window size: {} by {}'.format(objOffset.windowSizeHeight, objOffset.windowSizeWidth)) + # search range objOffset.halfSearchRangeDown = inps.srchgt objOffset.halfSearchRangeAcross = inps.srcwidth + print('half search range: {} by {}'.format(inps.srchgt, inps.srcwidth)) # starting pixel objOffset.masterStartPixelDownStatic = inps.margin objOffset.masterStartPixelAcrossStatic = inps.margin - - # skip size + + # skip/step size objOffset.skipSampleDown = inps.skiphgt objOffset.skipSampleAcross = inps.skipwidth + print('search step: {} by {}'.format(inps.skiphgt, inps.skipwidth)) + + # oversample raw data (SLC) + objOffset.rawDataOversamplingFactor = inps.raw_oversample + print('raw data oversampling factor:', inps.raw_oversample) + + # correlation surface + if inps.corr_win_size == -1: + corr_win_size_orig = min(inps.srchgt, inps.srcwidth) * inps.raw_oversample + 1 + inps.corr_win_size = np.power(2, int(np.log2(corr_win_size_orig))) + objOffset.corrSurfaceZoomInWindow = inps.corr_win_size + print('correlation surface zoom-in window size:', inps.corr_win_size) - # oversampling objOffset.corrSufaceOverSamplingMethod = 0 - objOffset.corrSurfaceOverSamplingFactor = inps.oversample + objOffset.corrSurfaceOverSamplingFactor = inps.corr_oversample + print('correlation surface oversampling factor:', inps.corr_oversample) # output filenames objOffset.offsetImageName = str(inps.outprefix) + str(inps.outsuffix) + '.bip' objOffset.grossOffsetImageName = str(inps.outprefix) + str(inps.outsuffix) + '_gross.bip' objOffset.snrImageName = str(inps.outprefix) + str(inps.outsuffix) + '_snr.bip' objOffset.covImageName = str(inps.outprefix) + str(inps.outsuffix) + '_cov.bip' - print("offsetfield: ",objOffset.offsetImageName) print("gross offsetfield: ",objOffset.grossOffsetImageName) print("snr: ",objOffset.snrImageName) @@ -175,44 +204,45 @@ def estimateOffsetField(master, slave, inps=None): snrImageName = objOffset.snrImageName.decode('utf8') covImageName = objOffset.covImageName.decode('utf8') - if os.path.exists(offsetImageName) and inps.redo==0: - + print(offsetImageName) + print(inps.redo) + if os.path.exists(offsetImageName) and not inps.redo: print('offsetfield file exists') - exit() + return 0 # generic control objOffset.numberWindowDownInChunk = inps.numWinDownInChunk objOffset.numberWindowAcrossInChunk = inps.numWinAcrossInChunk objOffset.useMmap = 0 objOffset.mmapSize = 8 - objOffset.setupParams() - + ## Set Gross Offset ### if inps.gross == 0: print("Set constant grossOffset") print("By default, the gross offsets are zero") print("You can override the default values here") objOffset.setConstantGrossOffset(0, 0) + else: print("Set varying grossOffset") print("By default, the gross offsets are zero") print("You can override the default grossDown and grossAcross arrays here") - objOffset.setVaryingGrossOffset(np.zeros(shape=grossDown.shape,dtype=np.int32), np.zeros(shape=grossAcross.shape,dtype=np.int32)) - - # check + objOffset.setVaryingGrossOffset(np.zeros(shape=grossDown.shape,dtype=np.int32), + np.zeros(shape=grossAcross.shape,dtype=np.int32)) + + # check objOffset.checkPixelInImageRange() # Run the code print('Running PyCuAmpcor') - - objOffset.runAmpcor() + objOffset.runAmpcor() print('Finished') sar.finalizeImage() sim.finalizeImage() - + # Finalize the results # offsetfield outImg = isceobj.createImage() @@ -257,18 +287,20 @@ def estimateOffsetField(master, slave, inps=None): covImg.setAccessMode('read') covImg.renderHdr() - return 0 - -def main(iargs=None): + return + + +def main(iargs=None): inps = cmdLineParse(iargs) outDir = os.path.dirname(inps.outprefix) print(inps.outprefix) - if not os.path.exists(outDir): - os.makedirs(outDir) - + os.makedirs(outDir, exist_ok=True) + estimateOffsetField(inps.master, inps.slave, inps) + return + if __name__ == '__main__': - + main() diff --git a/contrib/stack/stripmapStack/MaskAndFilter.py b/contrib/stack/stripmapStack/MaskAndFilter.py index 06d4417..d9277a7 100755 --- a/contrib/stack/stripmapStack/MaskAndFilter.py +++ b/contrib/stack/stripmapStack/MaskAndFilter.py @@ -3,8 +3,8 @@ # Author: Heresh Fattahi # Copyright 2016 # - -import os + +import os import argparse import numpy as np from scipy import ndimage @@ -32,7 +32,14 @@ GDAL2NUMPY_DATATYPE = { 7 : np.float64, 10: np.complex64, 11: np.complex128, -} +} + + +EXAMPLE = '''example: + MaskAndFilter.py -d offset.bip -s offset_snr.bip + MaskAndFilter.py -d offset.bip -s offset_snr.bip --plot +''' + EXAMPLE = '''example: @@ -42,10 +49,10 @@ EXAMPLE = '''example: def createParser(): - ''' + ''' Command line parser. ''' - + parser = argparse.ArgumentParser(description='Mask and filter the densOffset', formatter_class=argparse.RawTextHelpFormatter, epilog=EXAMPLE) @@ -88,7 +95,7 @@ def cmdLineParse(iargs = None): def read(file, processor='ISCE', bands=None, dataType=None): ''' raeder based on GDAL. - + Args: * file -> File name to be read @@ -96,7 +103,7 @@ def read(file, processor='ISCE', bands=None, dataType=None): Kwargs: * processor -> the processor used for the InSAR processing. default: ISCE - * bands -> a list of bands to be extracted. If not specified all bands will be extracted. + * bands -> a list of bands to be extracted. If not specified all bands will be extracted. * dataType -> if not specified, it will be extracted from the data itself Returns: * data : A numpy array with dimensions : number_of_bands * length * width @@ -116,7 +123,7 @@ def read(file, processor='ISCE', bands=None, dataType=None): if bands is None: bands = range(1,dataset.RasterCount+1) ###################################### - # if dataType is not known let's get it from the data: + # if dataType is not known let's get it from the data: if dataType is None: band = dataset.GetRasterBand(1) dataType = GDAL2NUMPY_DATATYPE[band.DataType] @@ -149,13 +156,13 @@ def fill(data, invalid=None): """ Replace the value of invalid 'data' cells (indicated by 'invalid') by the value of the nearest valid data cell - + Input: data: numpy array of any dimension invalid: a binary array of same shape as 'data'. data value are replaced where invalid is True If None (default), use: invalid = np.isnan(data) - + Output: Return a filled array. """ @@ -171,7 +178,7 @@ def mask_filter(inps, band, outName): """masking and Filtering""" # read offset - offset = read(inps.denseOffset, bands=band) + offset = read(inps.denseOffset, bands=band) offset = offset[0,:,:] # read SNR @@ -184,17 +191,21 @@ def mask_filter(inps, band, outName): offset1 = np.array(offset) offset1[snr < inps.snrThreshold] = np.nan + # percentage of masked out pixels among all non-zero SNR pixels + perc = np.sum(snr >= inps.snrThreshold) / np.sum(snr > 0) + print('percentage of pixels with SNR >= {} among pixels with SNR > 0: {:.0%}'.format(inps.snrThreshold, perc)) + # fill the hole in offset with nearest data print('fill the masked out region with nearest data') offset2 = fill(offset1) # median filtering - print('filtering with median filter with size : ', inps.filterSize) + print('filtering with median filter with size: {}'.format(inps.filterSize)) offset3 = ndimage.median_filter(offset2, size=inps.filterSize) length, width = offset3.shape # write data to file - print('writing masked and filtered offsets to: ', outName) + print('writing masked and filtered offsets to: {}'.format(outName)) write(offset3, outName, 1, 6) # write the xml/vrt/hdr file @@ -216,33 +227,55 @@ def mask_filter(inps, band, outName): def plot_mask_and_filtering(az_list, rg_list, inps=None): + print('-'*30) + print('plotting mask and filtering result ...') + print('mask pixels with SNR == 0 (for plotting ONLY; data files are untouched)') + snr = az_list[0] + for i in range(1, len(az_list)): + az_list[i][snr == 0] = np.nan + rg_list[i][snr == 0] = np.nan + + # percentage of masked out pixels among all non-zero SNR pixels + perc = np.sum(snr >= inps.snrThreshold) / np.sum(snr > 0) + print('percentage of pixels with SNR >= {} among pixels with SNR > 0: {:.0%}'.format(inps.snrThreshold, perc)) + fig, axs = plt.subplots(nrows=2, ncols=5, figsize=inps.figsize, sharex=True, sharey=True) - titles = ['SNR', 'offset', 'offset (mask)', 'offset (mask/fill)', 'offset (mask/fill/filter)'] + titles = ['SNR', + 'offset', + 'offset (mask {} - {:.0%} remain)'.format(inps.snrThreshold, perc), + 'offset (mask {} / fill)'.format(inps.snrThreshold), + 'offset (mask {} / fill / filter {})'.format(inps.snrThreshold, inps.filterSize)] # plot SNR - im0 = axs[0,0].imshow(az_list[0], vmin=inps.vlim_snr[0], vmax=inps.vlim_snr[1], cmap='RdBu') - im0 = axs[1,0].imshow(rg_list[0], vmin=inps.vlim_snr[0], vmax=inps.vlim_snr[1], cmap='RdBu') + kwargs = dict(vmin=inps.vlim_snr[0], vmax=inps.vlim_snr[1], cmap='RdBu', interpolation='nearest') + im0 = axs[0,0].imshow(snr, **kwargs) + im0 = axs[1,0].imshow(snr, **kwargs) axs[0,0].set_title('SNR', fontsize=12) + print('SNR data range: [{}, {}]'.format(np.nanmin(snr), np.nanmax(snr))) # label axs[0,0].set_ylabel('azimuth', fontsize=12) axs[1,0].set_ylabel('range', fontsize=12) # plot offset + kwargs = dict(vmin=inps.vlim[0], vmax=inps.vlim[1], cmap='jet', interpolation='nearest') for i in range(1,len(az_list)): - im1 = axs[0,i].imshow(az_list[i], vmin=inps.vlim[0], vmax=inps.vlim[1], cmap='jet') - im1 = axs[1,i].imshow(rg_list[i], vmin=inps.vlim[0], vmax=inps.vlim[1], cmap='jet') + im1 = axs[0,i].imshow(az_list[i], **kwargs) + im1 = axs[1,i].imshow(rg_list[i], **kwargs) axs[0,i].set_title(titles[i], fontsize=12) + print('{} data range'.format(titles[i])) + print('azimuth offset: [{:.3f}, {:.3f}]'.format(np.nanmin(az_list[i]), np.nanmax(az_list[i]))) + print('range offset: [{:.3f}, {:.3f}]'.format(np.nanmin(rg_list[i]), np.nanmax(rg_list[i]))) fig.tight_layout() # colorbar fig.subplots_adjust(bottom=0.15) - cax0 = fig.add_axes([0.09, 0.1, 0.08, 0.015]) + cax0 = fig.add_axes([0.08, 0.1, 0.08, 0.015]) cbar0 = plt.colorbar(im0, cax=cax0, orientation='horizontal') cax0.yaxis.set_ticks_position('left') #fig.subplots_adjust(right=0.93) - cax1 = fig.add_axes([0.57, 0.1, 0.08, 0.015]) + cax1 = fig.add_axes([0.60, 0.1, 0.15, 0.015]) cbar1 = plt.colorbar(im1, cax=cax1, orientation='horizontal') cbar1.set_label('pixel', fontsize=12) @@ -259,7 +292,7 @@ def main(iargs=None): inps = cmdLineParse(iargs) - os.makedirs(inps.outD, exist_ok=True) + os.makedirs(inps.outDir, exist_ok=True) ####################### # masking the dense offsets based on SNR and median filter the masked offs @@ -283,4 +316,3 @@ if __name__ == '__main__': Main driver. ''' main() - diff --git a/contrib/stack/topsStack/grossOffsets.py b/contrib/stack/topsStack/grossOffsets.py index 180e82e..0247049 100755 --- a/contrib/stack/topsStack/grossOffsets.py +++ b/contrib/stack/topsStack/grossOffsets.py @@ -11,13 +11,13 @@ import isceobj from iscesys.Component.ProductManager import ProductManager as PM import numpy as np from netCDF4 import Dataset -from mpl_toolkits.basemap import Basemap import gdal from scipy.interpolate import interp2d, griddata import matplotlib.pyplot as plt + class grossOffsets: def __init__(self): @@ -116,6 +116,7 @@ class grossOffsets: #x,y = np.meshgrid(self.x0,self.y0) + #from mpl_toolkits.basemap import Basemap #self.AntVeloDataMap = Basemap(width=5600000,height=5600000,\ # resolution='l',projection='stere',\ # lat_ts=-71,lat_0=-90,lon_0=0)