diff --git a/components/isceobj/StripmapProc/runDenseOffsets.py b/components/isceobj/StripmapProc/runDenseOffsets.py index eb27366..0071b1e 100644 --- a/components/isceobj/StripmapProc/runDenseOffsets.py +++ b/components/isceobj/StripmapProc/runDenseOffsets.py @@ -62,7 +62,7 @@ def estimateOffsetField(master, slave, denseOffsetFileName, else: objOffset.setImageDataType2('real') - + objOffset.offsetImageName = denseOffsetFileName + '.bil' objOffset.snrImageName = denseOffsetFileName +'_snr.bil' objOffset.covImageName = denseOffsetFileName +'_cov.bil' @@ -75,11 +75,11 @@ def estimateOffsetField(master, slave, denseOffsetFileName, def runDenseOffsets(self): - if self.doDenseOffsets or self.doRubbersheeting: + if self.doDenseOffsets or self.doRubbersheetingAzimuth: if self.doDenseOffsets: print('Dense offsets explicitly requested') - if self.doRubbersheeting: + if self.doRubbersheetingAzimuth: print('Generating offsets as rubber sheeting requested') else: return @@ -96,7 +96,7 @@ def runDenseOffsets(self): os.makedirs(dirname) denseOffsetFilename = os.path.join(dirname , self.insar.denseOffsetFilename) - + field = estimateOffsetField(masterSlc, slaveSlc, denseOffsetFilename, ww = self.denseWindowWidth, wh = self.denseWindowHeight, @@ -107,5 +107,5 @@ def runDenseOffsets(self): self._insar.offset_top = field[0] self._insar.offset_left = field[1] - + return None diff --git a/components/isceobj/StripmapProc/runRubbersheet.py b/components/isceobj/StripmapProc/runRubbersheet.py index ea417bd..6f6afc8 100644 --- a/components/isceobj/StripmapProc/runRubbersheet.py +++ b/components/isceobj/StripmapProc/runRubbersheet.py @@ -13,13 +13,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. """ @@ -97,7 +97,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): ###Currently making the assumption that top left of dense offsets and interfeorgrams are the same. ###This is not true for now. We need to update DenseOffsets to have the ability to have same top left ###As the input images. Once that is implemente, the math here should all be consistent. - ###However, this is not too far off since the skip for doing dense offsets is generally large. + ###However, this is not too far off since the skip for doing dense offsets is generally large. ###The offset is not too large to worry about right now. If the skip is decreased, this could be an issue. print('oversampling the filtered and masked offsets to the width and length:', width, ' ', length ) @@ -125,7 +125,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): for ll in range(length): val = geomoff.bands[0][ll,:] + osoff.bands[0][ll,:] val.tofile(fid) - + fid.close() img = isceobj.createImage() @@ -138,13 +138,13 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): img.scheme = 'BIP' img.renderHdr() - + return None def runRubbersheet(self): - if not self.doRubbersheeting: + if not self.doRubbersheetingAzimuth: print('Rubber sheeting not requested ... skipping') return @@ -168,11 +168,9 @@ def runRubbersheet(self): sheetOffset = os.path.join(offsetsDir, self.insar.azimuthRubbersheetFilename) # oversampling the filtAzOffsetFile to the same size of geometryAzimuthOffset - # and then update the geometryAzimuthOffset by adding the oversampled + # and then update the geometryAzimuthOffset by adding the oversampled # filtAzOffsetFile to it. resampleOffset(filtAzOffsetFile, geometryAzimuthOffset, sheetOffset) print("I'm here") return None - - diff --git a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py index 43ee4a8..fc07e25 100755 --- a/components/isceobj/StripmapProc/runRubbersheetAzimuth.py +++ b/components/isceobj/StripmapProc/runRubbersheetAzimuth.py @@ -22,39 +22,39 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): off_az = ds.GetRasterBand(1).ReadAsArray() off_rg = ds.GetRasterBand(2).ReadAsArray() ds = None - - # Remove missing values from ampcor + + # Remove missing values from ampcor off_rg[np.where(off_rg < -9999)]=0 off_az[np.where(off_az < -9999)]=0 - - + + # Store the offsets in a complex variable off = off_rg + 1j*off_az # Mask the azimuth offsets based on the MAD - mask = off_masking(off,filterSize,thre=3) - + mask = off_masking(off,filterSize,thre=3) + xoff_masked = np.ma.array(off.real,mask=mask) yoff_masked = np.ma.array(off.imag,mask=mask) - + # Delete unused variables mask = None off = None - + # Remove residual noisy spots with a median filter on the azimuth offmap yoff_masked.mask = yoff_masked.mask | \ (ndimage.median_filter(xoff_masked.filled(fill_value=0),3) == 0) | \ (ndimage.median_filter(yoff_masked.filled(fill_value=0),3) == 0) - - # Fill the data by iteratively using smoothed values + + # Fill the data by iteratively using smoothed values data = yoff_masked.data data[yoff_masked.mask]=np.nan - + off_az_filled = fill_with_smoothed(data,filterSize) - + # Apply median filter to smooth the azimuth offset map off_az_filled = ndimage.median_filter(off_az_filled,filterSize) - + # Save the filtered offsets length, width = off_az_filled.shape @@ -74,8 +74,8 @@ def mask_filterNoSNR(denseOffsetFile,filterSize,outName): img.scheme = 'BIP' img.renderHdr() - - return + + return def off_masking(off,filterSize,thre=2): @@ -94,16 +94,18 @@ 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. """ + from scipy import ndimage + if invalid is None: invalid = np.isnan(data) ind = ndimage.distance_transform_edt(invalid, @@ -166,7 +168,7 @@ def fill_with_smoothed(off,filterSize): kernel = np.ones((filterSize,filterSize),np.float32)/(filterSize*filterSize) loop = 0 cnt2=1 - + while (cnt2!=0 & loop<100): loop += 1 idx2= np.isnan(off_2filt) @@ -178,9 +180,9 @@ def fill_with_smoothed(off,filterSize): idx3 = np.where(off_filt == 0) off_2filt[idx3]=np.nan off_filt=None - + return off_2filt - + def resampleOffset(maskedFiltOffset, geometryOffset, outName): ''' Oversample offset and add. @@ -198,7 +200,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): ###Currently making the assumption that top left of dense offsets and interfeorgrams are the same. ###This is not true for now. We need to update DenseOffsets to have the ability to have same top left ###As the input images. Once that is implemente, the math here should all be consistent. - ###However, this is not too far off since the skip for doing dense offsets is generally large. + ###However, this is not too far off since the skip for doing dense offsets is generally large. ###The offset is not too large to worry about right now. If the skip is decreased, this could be an issue. print('oversampling the filtered and masked offsets to the width and length:', width, ' ', length ) @@ -226,7 +228,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): for ll in range(length): val = geomoff.bands[0][ll,:] + osoff.bands[0][ll,:] val.tofile(fid) - + fid.close() img = isceobj.createImage() @@ -239,7 +241,7 @@ def resampleOffset(maskedFiltOffset, geometryOffset, outName): img.scheme = 'BIP' img.renderHdr() - + return None @@ -264,7 +266,7 @@ def runRubbersheetAzimuth(self): if not self.doRubbersheetingRange: print('Rubber sheeting in range is off, filtering the offsets with a SNR-based mask') mask_filter(denseOffsetFile, snrFile, band[0], snrThreshold, filterSize, filtAzOffsetFile) - else: + else: print('Rubber sheeting in range is on, filtering the offsets with data-based mask') mask_filterNoSNR(denseOffsetFile, filterSize, filtAzOffsetFile) @@ -274,10 +276,8 @@ def runRubbersheetAzimuth(self): sheetOffset = os.path.join(offsetsDir, self.insar.azimuthRubbersheetFilename) # oversampling the filtAzOffsetFile to the same size of geometryAzimuthOffset - # and then update the geometryAzimuthOffset by adding the oversampled + # and then update the geometryAzimuthOffset by adding the oversampled # filtAzOffsetFile to it. resampleOffset(filtAzOffsetFile, geometryAzimuthOffset, sheetOffset) return None - -