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/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/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