diff --git a/applications/stripmapApp.py b/applications/stripmapApp.py index fbbeed5..938cafb 100755 --- a/applications/stripmapApp.py +++ b/applications/stripmapApp.py @@ -242,14 +242,20 @@ FILTER_STRENGTH = Application.Parameter('filterStrength', mandatory=False, doc='') - -DO_RUBBERSHEETING = Application.Parameter('doRubbersheeting', - public_name='do rubbersheeting', +############################################## Modified by V.Brancato 10.07.2019 +DO_RUBBERSHEETINGAZIMUTH = Application.Parameter('doRubbersheetingAzimuth', + public_name='do rubbersheetingAzimuth', default=False, type=bool, mandatory=False, doc='') - +DO_RUBBERSHEETINGRANGE = Application.Parameter('doRubbersheetingRange', + public_name='do rubbersheetingRange', + default=False, + type=bool, + mandatory=False, + doc='') +################################################################################# RUBBERSHEET_SNR_THRESHOLD = Application.Parameter('rubberSheetSNRThreshold', public_name='rubber sheet SNR Threshold', default = 5.0, @@ -533,7 +539,8 @@ class _RoiBase(Application, FrameMixin): GEOCODE_BOX, REGION_OF_INTEREST, HEIGHT_RANGE, - DO_RUBBERSHEETING, + DO_RUBBERSHEETINGRANGE, #Modified by V. Brancato 10.07.2019 + DO_RUBBERSHEETINGAZIMUTH, #Modified by V. Brancato 10.07.2019 RUBBERSHEET_SNR_THRESHOLD, RUBBERSHEET_FILTER_SIZE, DO_DENSEOFFSETS, @@ -724,7 +731,8 @@ class _RoiBase(Application, FrameMixin): self.runResampleSlc = StripmapProc.createResampleSlc(self) self.runRefineSlaveTiming = StripmapProc.createRefineSlaveTiming(self) self.runDenseOffsets = StripmapProc.createDenseOffsets(self) - self.runRubbersheet = StripmapProc.createRubbersheet(self) + self.runRubbersheetRange = StripmapProc.createRubbersheetRange(self) #Modified by V. Brancato 10.07.2019 + self.runRubbersheetAzimuth =StripmapProc.createRubbersheetAzimuth(self) #Modified by V. Brancato 10.07.2019 self.runResampleSubbandSlc = StripmapProc.createResampleSubbandSlc(self) self.runInterferogram = StripmapProc.createInterferogram(self) self.runFilter = StripmapProc.createFilter(self) @@ -774,8 +782,11 @@ class _RoiBase(Application, FrameMixin): args=('refined',)) self.step('dense_offsets', func=self.runDenseOffsets) - - self.step('rubber_sheet', func=self.runRubbersheet) +######################################################################## Modified by V. Brancato 10.07.2019 + self.step('rubber_sheet_range', func=self.runRubbersheetRange) + + self.step('rubber_sheet_azimuth',func=self.runRubbersheetAzimuth) +######################################################################### self.step('fine_resample', func=self.runResampleSlc, args=('fine',)) @@ -852,10 +863,14 @@ class _RoiBase(Application, FrameMixin): # run dense offsets self.runDenseOffsets() - - # adding the azimuth offsets computed from cross correlation to geometry offsets - self.runRubbersheet() - + +############ Modified by V. Brancato 10.07.2019 + # adding the azimuth offsets computed from cross correlation to geometry offsets + self.runRubbersheetAzimuth() + + # adding the range offsets computed from cross correlation to geometry offsets + self.runRubbersheetRange() +#################################################################################### # resampling using rubbersheeted offsets # which include geometry + constant range + constant azimuth # + dense azimuth offsets diff --git a/components/isceobj/StripmapProc/StripmapProc.py b/components/isceobj/StripmapProc/StripmapProc.py index ea9df09..c0f5344 100755 --- a/components/isceobj/StripmapProc/StripmapProc.py +++ b/components/isceobj/StripmapProc/StripmapProc.py @@ -325,14 +325,21 @@ AZIMUTH_OFFSET_FILENAME = Component.Parameter('azimuthOffsetFilename', doc='') - +# Modified by V. Brancato 10.07.2019 AZIMUTH_RUBBERSHEET_FILENAME = Component.Parameter('azimuthRubbersheetFilename', public_name='azimuth Rubbersheet Image Name', default = 'azimuth_sheet.off', type=str, mandatory=False, doc='') - + +RANGE_RUBBERSHEET_FILENAME = Component.Parameter('rangeRubbersheetFilename', + public_name='range Rubbersheet Image Name', + default = 'range_sheet.off', + type=str, + mandatory=False, + doc='') +# End of modification MISREG_FILENAME = Component.Parameter('misregFilename', public_name='misreg file name', default='misreg', @@ -346,14 +353,21 @@ DENSE_OFFSET_FILENAME = Component.Parameter('denseOffsetFilename', type=str, mandatory=False, doc='file name of dense offsets computed from cross correlating two SLC images') - +# Modified by V. Brancato 10.07.2019 FILT_AZIMUTH_OFFSET_FILENAME = Component.Parameter('filtAzimuthOffsetFilename', public_name='filtered azimuth offset filename', default='filtAzimuth.off', type=str, mandatory=False, doc='Filtered azimuth dense offsets') - + +FILT_RANGE_OFFSET_FILENAME = Component.Parameter('filtRangeOffsetFilename', + public_name='filtered range offset filename', + default='filtRange.off', + type=str, + mandatory=False, + doc='Filtered range dense offsets') +# End of modification DISPERSIVE_FILENAME = Component.Parameter('dispersiveFilename', public_name = 'dispersive phase filename', default='dispersive.bil', @@ -470,8 +484,10 @@ class StripmapProc(Component, FrameMixin): LOS_FILENAME, RANGE_OFFSET_FILENAME, AZIMUTH_OFFSET_FILENAME, - AZIMUTH_RUBBERSHEET_FILENAME, - FILT_AZIMUTH_OFFSET_FILENAME, + AZIMUTH_RUBBERSHEET_FILENAME, # Added by V. Brancato 10.07.2019 + RANGE_RUBBERSHEET_FILENAME, # Added by V. Brancato 10.07.2019 + FILT_AZIMUTH_OFFSET_FILENAME, # Added by V. Brancato 10.07.2019 + FILT_RANGE_OFFSET_FILENAME, # Added by V. Brancato 10.07.2019 DENSE_OFFSET_FILENAME, MISREG_FILENAME, DISPERSIVE_FILENAME, diff --git a/components/isceobj/StripmapProc/runInterferogram.py b/components/isceobj/StripmapProc/runInterferogram.py index 08bc98a..71e3ddd 100644 --- a/components/isceobj/StripmapProc/runInterferogram.py +++ b/components/isceobj/StripmapProc/runInterferogram.py @@ -1,14 +1,73 @@ + # # Author: Heresh Fattahi, 2017 -# +# Modified by V. Brancato (10.2019) +# (Included flattening when rubbersheeting in range is turned on import isceobj import logging from components.stdproc.stdproc import crossmul from iscesys.ImageUtil.ImageUtil import ImageUtil as IU import os +import gdal +import numpy as np + logger = logging.getLogger('isce.insar.runInterferogram') +# Added by V. Brancato 10.09.2019 +def write_xml(fileName,width,length,bands,dataType,scheme): + + img = isceobj.createImage() + img.setFilename(fileName) + img.setWidth(width) + img.setLength(length) + img.setAccessMode('READ') + img.bands = bands + img.dataType = dataType + img.scheme = scheme + img.renderHdr() + img.renderVRT() + + return None + + +def compute_FlatEarth(self,width,length): + from imageMath import IML + import logging + + # If rubbersheeting has been performed add back the range sheet offsets + + info = self._insar.loadProduct(self._insar.slaveSlcCropProduct) + radarWavelength = info.getInstrument().getRadarWavelength() + rangePixelSize = info.getInstrument().getRangePixelSize() + fact = 4 * np.pi* rangePixelSize / radarWavelength + + cJ = np.complex64(-1j) + + # Open the range sheet offset + rngOff = os.path.join(self.insar.offsetsDirname, self.insar.rangeOffsetFilename ) + + print(rngOff) + if os.path.exists(rngOff): + rng2 = np.memmap(rngOff, dtype=np.float64, mode='r', shape=(length,width)) + else: + print('No range offsets provided') + rng2 = np.zeros((length,width)) + + # Open the interferogram + ifgFilename= os.path.join(self.insar.ifgDirname, self.insar.ifgFilename) + intf = np.memmap(ifgFilename+'.full',dtype=np.complex64,mode='r+',shape=(length,width)) + + for ll in range(length): + intf[ll,:] *= np.exp(cJ*fact*rng2[ll,:]) + + del rng2 + del intf + + return + + + def multilook(infile, outname=None, alks=5, rlks=15): ''' Take looks. @@ -66,8 +125,8 @@ def computeCoherence(slc1name, slc2name, corname, virtual=True): slc2.finalizeImage() return - -def generateIgram(imageSlc1, imageSlc2, resampName, azLooks, rgLooks): +# Modified by V. Brancato on 10.09.2019 (added self) +def generateIgram(self,imageSlc1, imageSlc2, resampName, azLooks, rgLooks): objSlc1 = isceobj.createSlcImage() IU.copyAttributes(imageSlc1, objSlc1) objSlc1.setAccessMode('read') @@ -79,8 +138,13 @@ def generateIgram(imageSlc1, imageSlc2, resampName, azLooks, rgLooks): objSlc2.createImage() slcWidth = imageSlc1.getWidth() - intWidth = int(slcWidth / rgLooks) - + + + if not self.doRubbersheetingRange: + intWidth = int(slcWidth/rgLooks) # Modified by V. Brancato intWidth = int(slcWidth / rgLooks) + else: + intWidth = int(slcWidth) + lines = min(imageSlc1.getLength(), imageSlc2.getLength()) if '.flat' in resampName: @@ -93,7 +157,7 @@ def generateIgram(imageSlc1, imageSlc2, resampName, azLooks, rgLooks): resampInt = resampName objInt = isceobj.createIntImage() - objInt.setFilename(resampInt) + objInt.setFilename(resampInt+'.full') objInt.setWidth(intWidth) imageInt = isceobj.createIntImage() IU.copyAttributes(objInt, imageInt) @@ -101,21 +165,41 @@ def generateIgram(imageSlc1, imageSlc2, resampName, azLooks, rgLooks): objInt.createImage() objAmp = isceobj.createAmpImage() - objAmp.setFilename(resampAmp) + objAmp.setFilename(resampAmp+'.full') objAmp.setWidth(intWidth) imageAmp = isceobj.createAmpImage() IU.copyAttributes(objAmp, imageAmp) objAmp.setAccessMode('write') objAmp.createImage() + + if not self.doRubbersheetingRange: + print('Rubbersheeting in range is off, interferogram is already flattened') + objCrossmul = crossmul.createcrossmul() + objCrossmul.width = slcWidth + objCrossmul.length = lines + objCrossmul.LooksDown = azLooks + objCrossmul.LooksAcross = rgLooks - objCrossmul = crossmul.createcrossmul() - objCrossmul.width = slcWidth - objCrossmul.length = lines - objCrossmul.LooksDown = azLooks - objCrossmul.LooksAcross = rgLooks - - objCrossmul.crossmul(objSlc1, objSlc2, objInt, objAmp) - + objCrossmul.crossmul(objSlc1, objSlc2, objInt, objAmp) + else: + # Modified by V. Brancato 10.09.2019 (added option to add Range Rubber sheet Flat-earth back) + print('Rubbersheeting in range is on, removing flat-Earth phase') + objCrossmul = crossmul.createcrossmul() + objCrossmul.width = slcWidth + objCrossmul.length = lines + objCrossmul.LooksDown = 1 + objCrossmul.LooksAcross = 1 + objCrossmul.crossmul(objSlc1, objSlc2, objInt, objAmp) + + # Remove Flat-Earth component + compute_FlatEarth(self,intWidth,lines) + + # Perform Multilook + multilook(resampInt+'.full', outname=resampInt, alks=azLooks, rlks=rgLooks) #takeLooks(objAmp,azLooks,rgLooks) + multilook(resampAmp+'.full', outname=resampAmp, alks=azLooks, rlks=rgLooks) #takeLooks(objInt,azLooks,rgLooks) + + #os.system('rm ' + resampInt+'.full* ' + resampAmp + '.full* ') + # End of modification for obj in [objInt, objAmp, objSlc1, objSlc2]: obj.finalizeImage() @@ -142,7 +226,7 @@ def subBandIgram(self, masterSlc, slaveSlc, subBandDir): interferogramName = os.path.join(ifgDir , self.insar.ifgFilename) - generateIgram(img1, img2, interferogramName, azLooks, rgLooks) + generateIgram(self,img1, img2, interferogramName, azLooks, rgLooks) return interferogramName @@ -185,7 +269,7 @@ def runFullBandInterferogram(self): masterFrame = self._insar.loadProduct( self._insar.masterSlcCropProduct) masterSlc = masterFrame.getImage().filename - if self.doRubbersheeting: + if (self.doRubbersheetingRange | self.doRubbersheetingAzimuth): slaveSlc = os.path.join(self._insar.coregDirname, self._insar.fineCoregFilename) else: slaveSlc = os.path.join(self._insar.coregDirname, self._insar.refinedCoregFilename) @@ -212,7 +296,7 @@ def runFullBandInterferogram(self): interferogramName = os.path.join(ifgDir , self.insar.ifgFilename) - generateIgram(img1, img2, interferogramName, azLooks, rgLooks) + generateIgram(self,img1, img2, interferogramName, azLooks, rgLooks) ###Compute coherence @@ -221,7 +305,7 @@ def runFullBandInterferogram(self): multilook(cohname+'.full', outname=cohname, alks=azLooks, rlks=rgLooks) - ###Multilook relevant geometry products + ##Multilook relevant geometry products for fname in [self.insar.latFilename, self.insar.lonFilename, self.insar.losFilename]: inname = os.path.join(self.insar.geometryDirname, fname) multilook(inname + '.full', outname= inname, alks=azLooks, rlks=rgLooks) diff --git a/components/isceobj/StripmapProc/runResampleSlc.py b/components/isceobj/StripmapProc/runResampleSlc.py index cdc4a34..c108c98 100644 --- a/components/isceobj/StripmapProc/runResampleSlc.py +++ b/components/isceobj/StripmapProc/runResampleSlc.py @@ -23,7 +23,7 @@ def runResampleSlc(self, kind='coarse'): raise Exception('Unknown operation type {0} in runResampleSlc'.format(kind)) if kind == 'fine': - if not self.doRubbersheeting: + if not (self.doRubbersheetingRange | self.doRubbersheetingAzimuth): # Modified by V. Brancato 10.10.2019 print('Rubber sheeting not requested, skipping resampling ....') return @@ -68,12 +68,25 @@ def runResampleSlc(self, kind='coarse'): #Since the app is based on geometry module we expect pixel-by-pixel offset #field offsetsDir = self.insar.offsetsDirname - rgname = os.path.join(offsetsDir, self.insar.rangeOffsetFilename) + + # Modified by V. Brancato 10.10.2019 + #rgname = os.path.join(offsetsDir, self.insar.rangeOffsetFilename) + if kind in ['coarse', 'refined']: azname = os.path.join(offsetsDir, self.insar.azimuthOffsetFilename) + rgname = os.path.join(offsetsDir, self.insar.rangeOffsetFilename) else: azname = os.path.join(offsetsDir, self.insar.azimuthRubbersheetFilename) - + if self.doRubbersheetingRange: + print('Rubbersheeting in range is turned on, taking the cross-correlation offsets') + print('Setting Flattening to False') + rgname = os.path.join(offsetsDir, self.insar.rangeRubbersheetFilename) + flatten=False + else: + print('Rubbersheeting in range is turned off, taking range geometric offsets') + rgname = os.path.join(offsetsDir, self.insar.rangeOffsetFilename) + flatten=True + rngImg = isceobj.createImage() rngImg.load(rgname + '.xml') rngImg.setAccessMode('READ') @@ -85,8 +98,8 @@ def runResampleSlc(self, kind='coarse'): width = rngImg.getWidth() length = rngImg.getLength() - - flatten = True +# Modified by V. Brancato 10.10.2019 + #flatten = True rObj.flatten = flatten rObj.outputWidth = width rObj.outputLines = length diff --git a/components/isceobj/StripmapProc/runResampleSubbandSlc.py b/components/isceobj/StripmapProc/runResampleSubbandSlc.py index e3f7ff7..2cd7a12 100644 --- a/components/isceobj/StripmapProc/runResampleSubbandSlc.py +++ b/components/isceobj/StripmapProc/runResampleSubbandSlc.py @@ -14,7 +14,8 @@ import shelve logger = logging.getLogger('isce.insar.runResampleSubbandSlc') -def resampleSlc(masterFrame, slaveFrame, imageSlc2, radarWavelength, coregDir, +# Modified by V. Brancato 10.14.2019 added "self" as input parameter of resampleSLC +def resampleSlc(self,masterFrame, slaveFrame, imageSlc2, radarWavelength, coregDir, azoffname, rgoffname, azpoly = None, rgpoly = None, misreg=False): logger.info("Resampling slave SLC") @@ -56,8 +57,15 @@ def resampleSlc(masterFrame, slaveFrame, imageSlc2, radarWavelength, coregDir, width = rngImg.getWidth() length = rngImg.getLength() - - flatten = True +# Modified by V. Brancato on 10.14.2019 (if Rubbersheeting in range is turned on, flatten the interferogram during cross-correlation) + if not self.doRubbersheetingRange: + print('Rubber sheeting in range is turned off, flattening the interferogram during resampling') + flatten = True + else: + print('Rubber sheeting in range is turned on, flattening the interferogram during interferogram formation') + flatten=False +# end of Modification + rObj.flatten = flatten rObj.outputWidth = width rObj.outputLines = length @@ -105,15 +113,25 @@ def runResampleSubbandSlc(self, misreg=False): masterFrame = self._insar.loadProduct( self._insar.masterSlcCropProduct) slaveFrame = self._insar.loadProduct( self._insar.slaveSlcCropProduct) - if self.doRubbersheeting: - print('Using rubber sheeted offsets for resampling sub-bands') +# Modified by V. Brancato 10.14.2019 + + if self.doRubbersheetingAzimuth: + print('Using rubber in azimuth sheeted offsets for resampling sub-bands') azoffname = os.path.join( self.insar.offsetsDirname, self.insar.azimuthRubbersheetFilename) else: print('Using refined offsets for resampling sub-bands') azoffname = os.path.join( self.insar.offsetsDirname, self.insar.azimuthOffsetFilename) - rgoffname = os.path.join( self.insar.offsetsDirname, self.insar.rangeOffsetFilename) + if self.doRubbersheetingRange: + print('Using rubber in range sheeted offsets for resampling sub-bands') + rgoffname = os.path.join( self.insar.offsetsDirname, self.insar.rangeRubbersheetFilename) + else: + print('Using refined offsets for resampling sub-bands') + rgoffname = os.path.join( self.insar.offsetsDirname, self.insar.rangeOffsetFilename) +# ****************** End of Modification + + # rgoffname = os.path.join( self.insar.offsetsDirname, self.insar.rangeOffsetFilename) azpoly = self.insar.loadProduct( os.path.join(self.insar.misregDirname, self.insar.misregFilename) + '_az.xml') rgpoly = self.insar.loadProduct( os.path.join(self.insar.misregDirname, self.insar.misregFilename) + '_rg.xml') @@ -124,7 +142,7 @@ def runResampleSubbandSlc(self, misreg=False): wvlL = self.insar.lowBandRadarWavelength coregDir = os.path.join(self.insar.coregDirname, self.insar.lowBandSlcDirname) - lowbandCoregFilename = resampleSlc(masterFrame, slaveFrame, imageSlc2, wvlL, coregDir, + lowbandCoregFilename = resampleSlc(self,masterFrame, slaveFrame, imageSlc2, wvlL, coregDir, azoffname, rgoffname, azpoly=azpoly, rgpoly=rgpoly,misreg=False) imageSlc2 = os.path.join(self.insar.splitSpectrumDirname, self.insar.highBandSlcDirname, @@ -132,7 +150,7 @@ def runResampleSubbandSlc(self, misreg=False): wvlH = self.insar.highBandRadarWavelength coregDir = os.path.join(self.insar.coregDirname, self.insar.highBandSlcDirname) - highbandCoregFilename = resampleSlc(masterFrame, slaveFrame, imageSlc2, wvlH, coregDir, + highbandCoregFilename = resampleSlc(self,masterFrame, slaveFrame, imageSlc2, wvlH, coregDir, azoffname, rgoffname, azpoly=azpoly, rgpoly=rgpoly, misreg=False) self.insar.lowBandSlc2 = lowbandCoregFilename diff --git a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py new file mode 100755 index 0000000..5035de3 --- /dev/null +++ b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py @@ -0,0 +1,255 @@ +# +# Author: Heresh Fattahi +# Copyright 2017 +# + +import isce +import isceobj +from osgeo import gdal +from scipy import ndimage +from astropy.convolution import convolve +import numpy as np +import os + +def mask_filterNoSNR(denseOffsetFile,filterSize,outName): + # Masking the offsets with a data-based approach + + # Open the offsets + ds = gdal.Open(denseOffsetFile+'.vrt',gdal.GA_ReadOnly) + off_az = ds.GetRasterBand(1).ReadAsArray() + off_rg = ds.GetRasterBand(2).ReadAsArray() + ds = None + + off_az_copy = np.copy(off_az) + off_rg_copy = np.copy(off_rg) + + # Compute MAD + off_azm = ndimage.median_filter(off_az,filterSize) + off_rgm = ndimage.median_filter(off_rg,filterSize) + + metric_rg = np.abs(off_rgm-off_rg) + metric_az = np.abs(off_azm-off_az) + + idx = np.where((metric_rg > 3) | (metric_az > 3)) + + + # Remove missing values in the offset map + off_az_copy[np.where(off_az_copy<-9999)]=np.nan + off_az_copy[idx]=np.nan + + + # Fill the offsets with smoothed values + off_az_filled = fill_with_smoothed(off_az_copy,filterSize) + + # Median filter the offsets + off_az_filled = ndimage.median_filter(off_az_filled,filterSize) + + # Save the filtered offsets + length, width = off_az_filled.shape + + # writing the masked and filtered offsets to a file + print ('writing masked and filtered offsets to: ', outName) + + ##Write array to offsetfile + off_az_filled.tofile(outName) + + # write the xml file + img = isceobj.createImage() + img.setFilename(outName) + img.setWidth(width) + img.setAccessMode('READ') + img.bands = 1 + img.dataType = 'FLOAT' + img.scheme = 'BIP' + img.renderHdr() + + return None + + +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. + """ + if invalid is None: invalid = np.isnan(data) + + ind = ndimage.distance_transform_edt(invalid, + return_distances=False, + return_indices=True) + return data[tuple(ind)] + + +def mask_filter(denseOffsetFile, snrFile, band, snrThreshold, filterSize, outName): + #masking and Filtering + + ##Read in the offset file + ds = gdal.Open(denseOffsetFile + '.vrt', gdal.GA_ReadOnly) + Offset = ds.GetRasterBand(band).ReadAsArray() + ds = None + + ##Read in the SNR file + ds = gdal.Open(snrFile + '.vrt', gdal.GA_ReadOnly) + snr = ds.GetRasterBand(1).ReadAsArray() + ds = None + + # Masking the dense offsets based on SNR + print ('masking the dense offsets with SNR threshold: ', snrThreshold) + Offset[snr 3) | (metric_az > 3)) + + # Remove missing data in the offset map + off_rg_copy[np.where(off_rg_copy<-9999)]=np.nan + off_rg_copy[idx]=np.nan + + # Fill the offsets with smoothed values + off_rg_filled = fill_with_smoothed(off_rg_copy,filterSize) + + # Median filter the offsets + off_rg_filled = ndimage.median_filter(off_rg_filled,filterSize) + + # Save the filtered offsets + length, width = off_rg_filled.shape + + # writing the masked and filtered offsets to a file + print ('writing masked and filtered offsets to: ', outName) + + ##Write array to offsetfile + off_rg_filled.tofile(outName) + + # write the xml file + img = isceobj.createImage() + img.setFilename(outName) + img.setWidth(width) + img.setAccessMode('READ') + img.bands = 1 + img.dataType = 'FLOAT' + img.scheme = 'BIP' + img.renderHdr() + + return None + +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. + """ + if invalid is None: invalid = np.isnan(data) + + ind = ndimage.distance_transform_edt(invalid, + return_distances=False, + return_indices=True) + return data[tuple(ind)] + +def fill_with_smoothed(off,filterSize): + + off_2filt=np.copy(off) + kernel = np.ones((filterSize,filterSize),np.float32)/(filterSize*filterSize) + loop = 0 + cnt2=1 + + while (loop<100): + loop += 1 + idx2= np.isnan(off_2filt) + cnt2 = np.sum(np.count_nonzero(np.isnan(off_2filt))) + + if cnt2 != 0: + off_filt= convolve(off_2filt,kernel,boundary='extend',nan_treatment='interpolate') + off_2filt[idx2]=off_filt[idx2] + idx3 = np.where(off_filt == 0) + off_2filt[idx3]=np.nan + off_filt=None + + return off_2filt + +def mask_filter(denseOffsetFile, snrFile, band, snrThreshold, filterSize, outName): + #masking and Filtering + + ##Read in the offset file + ds = gdal.Open(denseOffsetFile + '.vrt', gdal.GA_ReadOnly) + Offset = ds.GetRasterBand(band).ReadAsArray() + ds = None + + ##Read in the SNR file + ds = gdal.Open(snrFile + '.vrt', gdal.GA_ReadOnly) + snr = ds.GetRasterBand(1).ReadAsArray() + ds = None + + # Masking the dense offsets based on SNR + print ('masking the dense offsets with SNR threshold: ', snrThreshold) + Offset[snr