diff --git a/applications/alos2App.py b/applications/alos2App.py index fec726f..c56b666 100755 --- a/applications/alos2App.py +++ b/applications/alos2App.py @@ -323,6 +323,14 @@ MASKED_AREAS_ION = Application.Parameter('maskedAreasIon', container = list, doc = 'areas masked out in ionospheric phase estimation') +SWATH_PHASE_DIFF_SNAP_ION = Application.Parameter('swathPhaseDiffSnapIon', + public_name = 'swath phase difference snap to fixed values', + default = None, + type = bool, + mandatory = False, + container = list, + doc = 'swath phase difference snap to fixed values') + FIT_ION = Application.Parameter('fitIon', public_name = 'apply polynomial fit before filtering ionosphere phase', default = True, @@ -330,16 +338,30 @@ FIT_ION = Application.Parameter('fitIon', mandatory = False, doc = 'apply polynomial fit before filtering ionosphere phase') +FILT_ION = Application.Parameter('filtIon', + public_name = 'whether filtering ionosphere phase', + default = True, + type = bool, + mandatory = False, + doc = 'whether filtering ionosphere phase') + +FIT_ADAPTIVE_ION = Application.Parameter('fitAdaptiveIon', + public_name = 'apply polynomial fit in adaptive filtering window', + default = True, + type = bool, + mandatory = False, + doc = 'apply polynomial fit in adaptive filtering window') + FILTERING_WINSIZE_MAX_ION = Application.Parameter('filteringWinsizeMaxIon', public_name='maximum window size for filtering ionosphere phase', - default=151, + default=301, type=int, mandatory=False, doc='maximum window size for filtering ionosphere phase') FILTERING_WINSIZE_MIN_ION = Application.Parameter('filteringWinsizeMinIon', public_name='minimum window size for filtering ionosphere phase', - default=41, + default=11, type=int, mandatory=False, doc='minimum window size for filtering ionosphere phase') @@ -608,7 +630,10 @@ class Alos2InSAR(Application): NUMBER_RANGE_LOOKS_ION, NUMBER_AZIMUTH_LOOKS_ION, MASKED_AREAS_ION, + SWATH_PHASE_DIFF_SNAP_ION, FIT_ION, + FILT_ION, + FIT_ADAPTIVE_ION, FILTERING_WINSIZE_MAX_ION, FILTERING_WINSIZE_MIN_ION, FILTER_SUBBAND_INT, diff --git a/applications/alos2burstApp.py b/applications/alos2burstApp.py index e078678..8470051 100755 --- a/applications/alos2burstApp.py +++ b/applications/alos2burstApp.py @@ -313,6 +313,14 @@ MASKED_AREAS_ION = Application.Parameter('maskedAreasIon', container = list, doc = 'areas masked out in ionospheric phase estimation') +SWATH_PHASE_DIFF_SNAP_ION = Application.Parameter('swathPhaseDiffSnapIon', + public_name = 'swath phase difference snap to fixed values', + default = None, + type = bool, + mandatory = False, + container = list, + doc = 'swath phase difference snap to fixed values') + FIT_ION = Application.Parameter('fitIon', public_name = 'apply polynomial fit before filtering ionosphere phase', default = True, @@ -320,16 +328,30 @@ FIT_ION = Application.Parameter('fitIon', mandatory = False, doc = 'apply polynomial fit before filtering ionosphere phase') +FILT_ION = Application.Parameter('filtIon', + public_name = 'whether filtering ionosphere phase', + default = True, + type = bool, + mandatory = False, + doc = 'whether filtering ionosphere phase') + +FIT_ADAPTIVE_ION = Application.Parameter('fitAdaptiveIon', + public_name = 'apply polynomial fit in adaptive filtering window', + default = True, + type = bool, + mandatory = False, + doc = 'apply polynomial fit in adaptive filtering window') + FILTERING_WINSIZE_MAX_ION = Application.Parameter('filteringWinsizeMaxIon', public_name='maximum window size for filtering ionosphere phase', - default=151, + default=301, type=int, mandatory=False, doc='maximum window size for filtering ionosphere phase') FILTERING_WINSIZE_MIN_ION = Application.Parameter('filteringWinsizeMinIon', public_name='minimum window size for filtering ionosphere phase', - default=41, + default=11, type=int, mandatory=False, doc='minimum window size for filtering ionosphere phase') @@ -543,7 +565,10 @@ class Alos2burstInSAR(Application): NUMBER_RANGE_LOOKS_ION, NUMBER_AZIMUTH_LOOKS_ION, MASKED_AREAS_ION, + SWATH_PHASE_DIFF_SNAP_ION, FIT_ION, + FILT_ION, + FIT_ADAPTIVE_ION, FILTERING_WINSIZE_MAX_ION, FILTERING_WINSIZE_MIN_ION, FILTER_SUBBAND_INT, diff --git a/components/isceobj/Alos2Proc/Alos2ProcPublic.py b/components/isceobj/Alos2Proc/Alos2ProcPublic.py index ccf124f..edd4fcf 100644 --- a/components/isceobj/Alos2Proc/Alos2ProcPublic.py +++ b/components/isceobj/Alos2Proc/Alos2ProcPublic.py @@ -536,6 +536,7 @@ def waterBodyRadar(latFile, lonFile, wbdFile, wbdOutFile): latFp = open(latFile, 'rb') lonFp = open(lonFile, 'rb') wbdOutFp = open(wbdOutFile, 'wb') + wbdOutIndex = np.arange(width, dtype=np.int32) print("create water body in radar coordinates...") for i in range(length): if (((i+1)%200) == 0): @@ -551,7 +552,7 @@ def waterBodyRadar(latFile, lonFile, wbdFile, wbdOutFile): np.logical_and(sampleIndex>=0, sampleIndex<=demImage.width-1) ) #keep SRTM convention. water body. (0) --- land; (-1) --- water; (-2 or other value) --- no data. - wbdOut = wbd[(lineIndex[inboundIndex], sampleIndex[inboundIndex])] + wbdOut[(wbdOutIndex[inboundIndex],)] = wbd[(lineIndex[inboundIndex], sampleIndex[inboundIndex])] wbdOut.astype(np.int8).tofile(wbdOutFp) print("processing line %6d of %6d" % (length, length)) #create_xml(wbdOutFile, width, length, 'byte') @@ -1183,7 +1184,74 @@ def create_multi_index2(width2, l1, l2): return ((l2 - l1) / 2.0 + np.arange(width2) * l2) / l1 - +def computePhaseDiff(data1, data22, coherenceWindowSize=5, coherenceThreshold=0.85): + import copy + import numpy as np + from isceobj.Alos2Proc.Alos2ProcPublic import cal_coherence_1 + + #data22 will be changed in the processing, so make a copy here + data2 = copy.deepcopy(data22) + + dataDiff = data1 * np.conj(data2) + cor = cal_coherence_1(dataDiff, win=coherenceWindowSize) + index = np.nonzero(np.logical_and(cor>coherenceThreshold, dataDiff!=0)) + + #check if there are valid pixels + if index[0].size == 0: + phaseDiff = 0.0 + numberOfValidSamples = 0 + return (phaseDiff, numberOfValidSamples) + else: + numberOfValidSamples = index[0].size + + #in case phase difference is around PI, sum of +PI and -PI is zero, which affects the following + #mean phase difference computation. + #remove magnitude before doing sum? + dataDiff = dataDiff / (np.absolute(dataDiff)+(dataDiff==0)) + phaseDiff0 = np.angle(np.sum(dataDiff[index], dtype=np.complex128)) + #now the phase difference values are mostly centered at 0 + data2 *= np.exp(np.complex64(1j) * phaseDiff0) + phaseDiff = phaseDiff0 + + #compute phase difference + numberOfIterations = 1000000 + threshold = 0.000001 + for k in range(numberOfIterations): + dataDiff = data1 * np.conj(data2) + angle = np.mean(np.angle(dataDiff[index]), dtype=np.float64) + phaseDiff += angle + data2 *= np.exp(np.complex64(1j) * angle) + print('phase offset: %15.12f rad after iteration: %3d'%(phaseDiff, k+1)) + if (k+1 >= 5) and (angle <= threshold): + break + + #only take the value within -pi--pi + if phaseDiff > np.pi: + phaseDiff -= 2.0 * np.pi + if phaseDiff < -np.pi: + phaseDiff += 2.0 * np.pi + + # mean phase difference + # number of valid samples to compute the phase difference + return (phaseDiff, numberOfValidSamples) + + +def snap(inputValue, fixedValues, snapThreshold): + ''' + fixedValues can be a list or numpy array + ''' + import numpy as np + + diff = np.absolute(np.absolute(np.array(fixedValues)) - np.absolute(inputValue)) + indexMin = np.argmin(diff) + if diff[indexMin] < snapThreshold: + outputValue = np.sign(inputValue) * np.absolute(fixedValues[indexMin]) + snapped = True + else: + outputValue = inputValue + snapped = False + + return (outputValue, snapped) diff --git a/components/isceobj/Alos2Proc/runFrameMosaic.py b/components/isceobj/Alos2Proc/runFrameMosaic.py index 45f2464..152b6c0 100644 --- a/components/isceobj/Alos2Proc/runFrameMosaic.py +++ b/components/isceobj/Alos2Proc/runFrameMosaic.py @@ -103,13 +103,18 @@ def runFrameMosaic(self): rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=False, phaseCompensation=False, resamplingMethod=0) #mosaic interferograms - frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, + (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=True, phaseCompensation=True, resamplingMethod=1) create_xml(self._insar.amplitude, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'amp') create_xml(self._insar.interferogram, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'int') + catalog.addItem('frame phase diff estimated', phaseDiffEst[1:], 'runFrameMosaic') + catalog.addItem('frame phase diff used', phaseDiffUsed[1:], 'runFrameMosaic') + catalog.addItem('frame phase diff used source', phaseDiffSource[1:], 'runFrameMosaic') + catalog.addItem('frame phase diff samples used', numberOfValidSamples[1:], 'runFrameMosaic') + #update secondary parameters here #do not match for secondary, always use geometrical rangeOffsets = self._insar.frameRangeOffsetGeometricalSecondary @@ -125,7 +130,7 @@ def runFrameMosaic(self): self._insar.procDoc.addAllFromCatalog(catalog) -def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks, updateTrack=False, phaseCompensation=False, resamplingMethod=0): +def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks, updateTrack=False, phaseCompensation=False, phaseDiffFixed=None, snapThreshold=None, resamplingMethod=0): ''' mosaic frames @@ -138,6 +143,8 @@ def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num numberOfAzimuthLooks: number of azimuth looks of the input files updateTrack: whether update track parameters phaseCompensation: whether do phase compensation for each frame + phaseDiffFixed: if provided, the estimated value will snap to one of these values, which is nearest to the estimated one. + snapThreshold: this is used with phaseDiffFixed resamplingMethod: 0: amp resampling. 1: int resampling. 2: slc resampling ''' import numpy as np @@ -149,6 +156,8 @@ def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num from isceobj.Alos2Proc.Alos2ProcPublic import create_xml from isceobj.Alos2Proc.Alos2ProcPublic import find_vrt_file from isceobj.Alos2Proc.Alos2ProcPublic import find_vrt_keyword + from isceobj.Alos2Proc.Alos2ProcPublic import computePhaseDiff + from isceobj.Alos2Proc.Alos2ProcPublic import snap numberOfFrames = len(track.frames) frames = track.frames @@ -305,6 +314,15 @@ def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num #compute phase offset if phaseCompensation: + + phaseDiffEst = [0.0 for i in range(numberOfFrames)] + phaseDiffUsed = [0.0 for i in range(numberOfFrames)] + phaseDiffSource = ['estimated' for i in range(numberOfFrames)] + numberOfValidSamples = [0 for i in range(numberOfFrames)] + #phaseDiffEst = [0.0] + #phaseDiffUsed = [0.0] + #phaseDiffSource = ['estimated'] + phaseOffsetPolynomials = [np.array([0.0])] for i in range(1, numberOfFrames): upperframe = np.zeros((ye[i-1]-ys[i]+1, outWidth), dtype=np.complex128) @@ -323,8 +341,29 @@ def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num diff = np.sum(upperframe * np.conj(lowerframe), axis=0) (firstLine, lastLine, firstSample, lastSample) = findNonzero(np.reshape(diff, (1, outWidth))) #here i use mean value(deg=0) in case difference is around -pi or pi. + #!!!!!there have been updates, now deg must be 0 deg = 0 p = np.polyfit(np.arange(firstSample, lastSample+1), np.angle(diff[firstSample:lastSample+1]), deg) + + #need to use a more sophisticated method to compute the mean phase difference + (phaseDiffEst[i], numberOfValidSamples[i]) = computePhaseDiff(upperframe, lowerframe, coherenceWindowSize=9, coherenceThreshold=0.80) + + #snap phase difference to fixed values + if phaseDiffFixed is not None: + (outputValue, snapped) = snap(phaseDiffEst[i], phaseDiffFixed, snapThreshold) + if snapped == True: + phaseDiffUsed[i] = outputValue + phaseDiffSource[i] = 'estimated+snap' + else: + phaseDiffUsed[i] = phaseDiffEst[i] + phaseDiffSource[i] = 'estimated' + else: + phaseDiffUsed[i] = phaseDiffEst[i] + phaseDiffSource[i] = 'estimated' + + #use new phase constant value + p[-1] = phaseDiffUsed[i] + phaseOffsetPolynomials.append(p) @@ -435,6 +474,10 @@ def frameMosaic(track, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num track.azimuthPixelSize = frames[0].azimuthPixelSize track.azimuthLineInterval = frames[0].azimuthLineInterval + if phaseCompensation: + # estimated phase diff, used phase diff, used phase diff source + return (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) + def frameMosaicParameters(track, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks): ''' diff --git a/components/isceobj/Alos2Proc/runIonFilt.py b/components/isceobj/Alos2Proc/runIonFilt.py index ffe3d3c..53ccf50 100644 --- a/components/isceobj/Alos2Proc/runIonFilt.py +++ b/components/isceobj/Alos2Proc/runIonFilt.py @@ -113,53 +113,53 @@ def runIonFilt(self): ################################################# #SET PARAMETERS HERE - #if applying polynomial fitting - #False: no fitting, True: with fitting + #fit and filter ionosphere fit = self.fitIon - #gaussian filtering window size + filt = self.filtIon + fitAdaptive = self.fitAdaptiveIon + if (fit == False) and (filt == False): + raise Exception('either fit ionosphere or filt ionosphere should be True when doing ionospheric correction\n') + + #filtering window size size_max = self.filteringWinsizeMaxIon size_min = self.filteringWinsizeMinIon + if size_min > size_max: + print('\n\nWARNING: minimum window size for filtering ionosphere phase {} > maximum window size {}'.format(size_min, size_max)) + print(' re-setting maximum window size to {}\n\n'.format(size_min)) + size_max = size_min - if size_min >= size_max: - print('\n\nWARNING: minimum window size for filtering ionosphere phase {} >= maximum window size {}'.format(size_min, size_max)) - print(' resetting maximum window size to {}\n\n'.format(size_min+5)) - size_max = size_min + 5 - - #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self) - #corThresholdFit = 0.85 - - #Now changed to use lower band coherence. crl, 23-apr-2020. - useDiffCoherence = False - if useDiffCoherence: - #parameters for using diff coherence - corfile = 'diff'+ml2+'.cor' - corThresholdFit = 0.95 - # 1 is not good for low coherence case, changed to 20 - #corOrderFit = 1 - corOrderFit = 20 - corOrderFilt = 14 - else: - #parameters for using lower/upper band coherence - corfile = subbandPrefix[0]+ml2+'.cor' - corThresholdFit = 0.4 - corOrderFit = 10 - corOrderFilt = 4 + #coherence threshold for fitting a polynomial + corThresholdFit = 0.25 + #ionospheric phase standard deviation after filtering + std_out0 = 0.1 ################################################# print('\nfiltering ionosphere') + + #input files ionfile = 'ion'+ml2+'.ion' #corfile = 'diff'+ml2+'.cor' + corLowerfile = subbandPrefix[0]+ml2+'.cor' + corUpperfile = subbandPrefix[1]+ml2+'.cor' + #output files ionfiltfile = 'filt_ion'+ml2+'.ion' + stdfiltfile = 'filt_ion'+ml2+'.std' + windowsizefiltfile = 'filt_ion'+ml2+'.win' + #read data img = isceobj.createImage() img.load(ionfile + '.xml') width = img.width length = img.length - #ion = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :] + ion = np.fromfile(ionfile, dtype=np.float32).reshape(length, width) - cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :] - #amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :] + corLower = (np.fromfile(corLowerfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :] + corUpper = (np.fromfile(corUpperfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :] + cor = (corLower + corUpper) / 2.0 + index = np.nonzero(np.logical_or(corLower==0, corUpper==0)) + cor[index] = 0 + del corLower, corUpper #masked out user-specified areas if self.maskedAreasIon != None: @@ -172,7 +172,7 @@ def runIonFilt(self): cor[np.nonzero(cor<0)] = 0.0 cor[np.nonzero(cor>1)] = 0.0 - #remove water body + #remove water body. Not helpful, just leave it here wbd = np.fromfile('wbd'+ml2+'.wbd', dtype=np.int8).reshape(length, width) cor[np.nonzero(wbd==-1)] = 0.0 @@ -183,32 +183,113 @@ def runIonFilt(self): # wbd = np.fromfile(waterBodyFile, dtype=np.int8).reshape(length, width) # cor[np.nonzero(wbd!=0)] = 0.00001 - if fit: - import copy - wgt = copy.deepcopy(cor) - wgt[np.nonzero(wgt= 1 + + zero samples in data and weight are OK. + ''' + #import numpy as np + + if order < 1: + raise Exception('order must >= 1!\n') + + if data.shape != weight.shape: + raise Exception('data and weight must be of same size!\n') + + (length, width) = data.shape + #length*width, but below is better since no need to convert to int + n = data.size + + #number of coefficients + ncoeff = 1 + for i in range(1, order+1): + for j in range(i+1): + ncoeff += 1 + + #row, column + y, x = np.indices((length, width)) + x = x.flatten() + y = y.flatten() + z = data.flatten() + weight = np.sqrt(weight.flatten()) + + #linear functions: H theta = s + #compute observation matrix H (n*ncoeff) + H = np.zeros((n, ncoeff)) + H[:,0] += 1 + k = 1 + for i in range(1, order+1): + for j in range(i+1): + #x and y do not need to be column vector here + H[:, k] = x**(i-j)*y**(j) + k += 1 + + #least squares + #this is robust to singular cases + coeff = np.linalg.lstsq(H*weight[:,None], z*weight, rcond=-1)[0] + #this uses multiple threads, should be faster + #coeff = least_sqares(H*weight[:,None], z*weight, W=None) + + #fit surface + data_fit = (np.dot(H, coeff)).reshape(length, width) + + return (data_fit, coeff) + + +def adaptive_gaussian(data, std, size_min, size_max, std_out0, fit=True): + ''' + This program performs Gaussian filtering with adaptive window size. + Cunren Liang, 11-JUN-2020 + + data: input raw data, numpy array + std: standard deviation of raw data, numpy array + size_min: minimum filter window size + size_max: maximum filter window size (size_min <= size_max, size_min == size_max is allowed) + std_out0: standard deviation of output data + fit: whether do fitting before gaussian filtering + ''' + import scipy.signal as ss + + + (length, width) = data.shape + + #assume zero-value samples are invalid + index = np.nonzero(np.logical_or(data==0, std==0)) + data[index] = 0 + std[index] = 0 + #compute weight using standard deviation + wgt = 1.0 / (std**2 + (std==0)) + wgt[index] = 0 + + #compute number of gaussian filters + if size_min > size_max: + raise Exception('size_min: {} > size_max: {}\n'.format(size_min, size_max)) + + if size_min % 2 == 0: + size_min += 1 + if size_max % 2 == 0: + size_max += 1 + + size_num = int((size_max - size_min) / 2 + 1) + #'size_num == 1' is checked to be OK starting from here + + + #create gaussian filters + print('compute Gaussian filters\n') + gaussian_filters = [] + for i in range(size_num): + size = int(size_min + i * 2) + gaussian_filters.append(gaussian(size, size/2.0, scale=1.0)) + + + #compute standard deviation after filtering coresponding to each of gaussian_filters + #if value is 0, there is no valid sample in the gaussian window + print('compute standard deviation after filtering for each filtering window size') + std_filt = np.zeros((length, width, size_num)) + for i in range(size_num): + size = int(size_min + i * 2) + print('current window size: %4d, min window size: %4d, max window size: %4d' % (size, size_min, size_max), end='\r', flush=True) + #robust zero value detector. non-zero convolution result at least >= 1, so can use 0.5 + #as threshold to detect zero-value result + index = np.nonzero(ss.fftconvolve(wgt!=0, gaussian_filters[i]!=0, mode='same') < 0.5) + scale = ss.fftconvolve(wgt, gaussian_filters[i], mode='same') + scale[index] = 0 + #variance of resulting filtered sample + var_filt = ss.fftconvolve(wgt, gaussian_filters[i]**2, mode='same') / (scale**2 + (scale==0)) + var_filt[index] = 0 + std_filt[:, :, i] = np.sqrt(var_filt) + print('\n') + + + #find gaussian window size (3rd-dimension index of the window size in gaussian_filters) + #if value is -1, there is no valid sample in any of the gaussian windows + #and therefore no filtering in the next step is needed + print('find Gaussian window size to use') + gaussian_index = np.zeros((length, width), dtype=np.int32) + std_filt2 = np.zeros((length, width)) + for i in range(length): + if (((i+1)%50) == 0): + print('processing line %6d of %6d' % (i+1, length), end='\r', flush=True) + for j in range(width): + if np.sum(std_filt[i, j, :]) == 0: + gaussian_index[i, j] = -1 + else: + gaussian_index[i, j] = size_num - 1 + for k in range(size_num): + if (std_filt[i, j, k] != 0) and (std_filt[i, j, k] <= std_out0): + gaussian_index[i, j] = k + break + if gaussian_index[i, j] != -1: + std_filt2[i, j] = std_filt[i, j, gaussian_index[i, j]] + del std_filt + print("processing line %6d of %6d\n" % (length, length)) + + + #adaptive gaussian filtering + print('filter image') + data_out = np.zeros((length, width)) + std_out = np.zeros((length, width)) + window_size_out = np.zeros((length, width), dtype=np.int16) + for i in range(length): + #if (((i+1)%5) == 0): + print('processing line %6d of %6d' % (i+1, length), end='\r', flush=True) + for j in range(width): + #if value is -1, there is no valid sample in any of the gaussian windows + #and therefore no filtering in the next step is needed + if gaussian_index[i, j] == -1: + continue + + #1. extract data + size = int(size_min + gaussian_index[i, j] * 2) + size_half = int((size - 1) / 2) + window_size_out[i, j] = size + + #index in original data + first_line = max(i-size_half, 0) + last_line = min(i+size_half, length-1) + first_column = max(j-size_half, 0) + last_column = min(j+size_half, width-1) + length_valid = last_line - first_line + 1 + width_valid = last_column - first_column + 1 + + #index in filter window + if first_line == 0: + last_line2 = size - 1 + first_line2 = last_line2 - (length_valid - 1) + else: + first_line2 = 0 + last_line2 = first_line2 + (length_valid - 1) + if first_column == 0: + last_column2 = size - 1 + first_column2 = last_column2 - (width_valid - 1) + else: + first_column2 = 0 + last_column2 = first_column2 + (width_valid - 1) + + #prepare data and weight within the window + data_window = np.zeros((size, size)) + wgt_window = np.zeros((size, size)) + data_window[first_line2:last_line2+1, first_column2:last_column2+1] = data[first_line:last_line+1, first_column:last_column+1] + wgt_window[first_line2:last_line2+1, first_column2:last_column2+1] = wgt[first_line:last_line+1, first_column:last_column+1] + #number of valid samples in the filtering window + n_valid = np.sum(data_window!=0) + + #2. fit + #order, n_coeff = (1, 3) + order, n_coeff = (2, 6) + if fit: + #must have enough samples to do fitting + #even if order is 2, n_coeff * 3 is much smaller than size_min*size_min in most cases. + if n_valid > n_coeff * 3: + #data_fit = weight_fitting(data_window, wgt_window, size, size, 1, 1, 1, 1, order) + data_fit, coeff = polyfit_2d(data_window, wgt_window, order) + index = np.nonzero(data_window!=0) + data_window[index] -= data_fit[index] + + #3. filter + wgt_window_2 = wgt_window * gaussian_filters[gaussian_index[i, j]] + scale = 1.0/np.sum(wgt_window_2) + wgt_window_2 *= scale + data_out[i, j] = np.sum(wgt_window_2 * data_window) + #std_out[i, j] = scale * np.sqrt(np.sum(wgt_window*(gaussian_filters[gaussian_index[i, j]]**2))) + #already computed + std_out[i, j] = std_filt2[i, j] + #print('std_out[i, j], std_filt2[i, j]', std_out[i, j], std_filt2[i, j]) + + #4. add back filtered value + if fit: + if n_valid > n_coeff * 3: + data_out[i, j] += data_fit[size_half, size_half] + print('\n') + + return (data_out, std_out, window_size_out) + + +def reformatMaskedAreas(maskedAreas, length, width): + ''' + reformat masked areas coordinates that are ready to use + 'maskedAreas' is a 2-D list. Each element in the 2-D list is a four-element list: [firstLine, + lastLine, firstColumn, lastColumn], with line/column numbers starting with 1. If one of the + four elements is specified with -1, the program will use firstLine/lastLine/firstColumn/ + lastColumn instead. + + output is a 2-D list containing the corresponding python-list/array-format indexes. + ''' + numberOfAreas = len(maskedAreas) + maskedAreasReformated = [[0, length, 0, width] for i in range(numberOfAreas)] + + for i in range(numberOfAreas): + if maskedAreas[i][0] != -1: + maskedAreasReformated[i][0] = maskedAreas[i][0] - 1 + if maskedAreas[i][1] != -1: + maskedAreasReformated[i][1] = maskedAreas[i][1] + if maskedAreas[i][2] != -1: + maskedAreasReformated[i][2] = maskedAreas[i][2] - 1 + if maskedAreas[i][3] != -1: + maskedAreasReformated[i][3] = maskedAreas[i][3] + if (not (0 <= maskedAreasReformated[i][0] <= length-1)) or \ + (not (1 <= maskedAreasReformated[i][1] <= length)) or \ + (not (0 <= maskedAreasReformated[i][2] <= width-1)) or \ + (not (1 <= maskedAreasReformated[i][3] <= width)) or \ + (not (maskedAreasReformated[i][1]-maskedAreasReformated[i][0]>=1)) or \ + (not (maskedAreasReformated[i][3]-maskedAreasReformated[i][2]>=1)): + raise Exception('area {} masked out in ionospheric phase estimation not correct'.format(i+1)) + + return maskedAreasReformated + + +######################################################################################################## +# The following functions are not used anywhere, and can be deleted. +######################################################################################################## + def fit_surface(x, y, z, wgt, order): # x: x coordinate, a column vector # y: y coordinate, a column vector @@ -530,106 +968,15 @@ def weight_fitting(ionos, weight, width, length, nrli, nali, nrlo, nalo, order): return phase_fit -def gaussian(size, sigma, scale = 1.0): - - if size % 2 != 1: - raise Exception('size must be odd') - hsize = (size - 1) / 2 - x = np.arange(-hsize, hsize + 1) * scale - f = np.exp(-x**2/(2.0*sigma**2)) / (sigma * np.sqrt(2.0*np.pi)) - f2d=np.matlib.repmat(f, size, 1) * np.matlib.repmat(f.reshape(size, 1), 1, size) - - return f2d/np.sum(f2d) -def adaptive_gaussian(ionos, wgt, size_max, size_min): - ''' - This program performs Gaussian filtering with adaptive window size. - ionos: ionosphere - wgt: weight - size_max: maximum window size - size_min: minimum window size - ''' - import scipy.signal as ss - - length = (ionos.shape)[0] - width = (ionos.shape)[1] - flag = (ionos!=0) * (wgt!=0) - ionos *= flag - wgt *= flag - - size_num = 100 - size = np.linspace(size_min, size_max, num=size_num, endpoint=True) - std = np.zeros((length, width, size_num)) - flt = np.zeros((length, width, size_num)) - out = np.zeros((length, width, 1)) - - #calculate filterd image and standard deviation - #sigma of window size: size_max - sigma = size_max / 2.0 - for i in range(size_num): - size2 = np.int(np.around(size[i])) - if size2 % 2 == 0: - size2 += 1 - if (i+1) % 10 == 0: - print('min win: %4d, max win: %4d, current win: %4d'%(np.int(np.around(size_min)), np.int(np.around(size_max)), size2)) - g2d = gaussian(size2, sigma*size2/size_max, scale=1.0) - scale = ss.fftconvolve(wgt, g2d, mode='same') - flt[:, :, i] = ss.fftconvolve(ionos*wgt, g2d, mode='same') / (scale + (scale==0)) - #variance of resulting filtered sample - scale = scale**2 - var = ss.fftconvolve(wgt, g2d**2, mode='same') / (scale + (scale==0)) - #in case there is a large area without data where scale is very small, which leads to wired values in variance - var[np.nonzero(var<0)] = 0 - std[:, :, i] = np.sqrt(var) - - std_mv = np.mean(std[np.nonzero(std!=0)], dtype=np.float64) - diff_max = np.amax(np.absolute(std - std_mv)) + std_mv + 1 - std[np.nonzero(std==0)] = diff_max - - index = np.nonzero(np.ones((length, width))) + ((np.argmin(np.absolute(std - std_mv), axis=2)).reshape(length*width), ) - out = flt[index] - out = out.reshape((length, width)) - - #remove artifacts due to varying wgt - size_smt = size_min - if size_smt % 2 == 0: - size_smt += 1 - g2d = gaussian(size_smt, size_smt/2.0, scale=1.0) - scale = ss.fftconvolve((out!=0), g2d, mode='same') - out2 = ss.fftconvolve(out, g2d, mode='same') / (scale + (scale==0)) - - return out2 -def reformatMaskedAreas(maskedAreas, length, width): - ''' - reformat masked areas coordinates that are ready to use - 'maskedAreas' is a 2-D list. Each element in the 2-D list is a four-element list: [firstLine, - lastLine, firstColumn, lastColumn], with line/column numbers starting with 1. If one of the - four elements is specified with -1, the program will use firstLine/lastLine/firstColumn/ - lastColumn instead. - output is a 2-D list containing the corresponding python-list/array-format indexes. - ''' - numberOfAreas = len(maskedAreas) - maskedAreasReformated = [[0, length, 0, width] for i in range(numberOfAreas)] - for i in range(numberOfAreas): - if maskedAreas[i][0] != -1: - maskedAreasReformated[i][0] = maskedAreas[i][0] - 1 - if maskedAreas[i][1] != -1: - maskedAreasReformated[i][1] = maskedAreas[i][1] - if maskedAreas[i][2] != -1: - maskedAreasReformated[i][2] = maskedAreas[i][2] - 1 - if maskedAreas[i][3] != -1: - maskedAreasReformated[i][3] = maskedAreas[i][3] - if (not (0 <= maskedAreasReformated[i][0] <= length-1)) or \ - (not (1 <= maskedAreasReformated[i][1] <= length)) or \ - (not (0 <= maskedAreasReformated[i][2] <= width-1)) or \ - (not (1 <= maskedAreasReformated[i][3] <= width)) or \ - (not (maskedAreasReformated[i][1]-maskedAreasReformated[i][0]>=1)) or \ - (not (maskedAreasReformated[i][3]-maskedAreasReformated[i][2]>=1)): - raise Exception('area {} masked out in ionospheric phase estimation not correct'.format(i+1)) - return maskedAreasReformated + + + + + diff --git a/components/isceobj/Alos2Proc/runIonSubband.py b/components/isceobj/Alos2Proc/runIonSubband.py index ea1b058..1991f4f 100644 --- a/components/isceobj/Alos2Proc/runIonSubband.py +++ b/components/isceobj/Alos2Proc/runIonSubband.py @@ -329,6 +329,16 @@ def runIonSubband(self): #These are for ALOS-2, may need to change for ALOS-4! phaseDiffFixed = [0.0, 0.4754024578084084, 0.9509913179406437, 1.4261648478671614, 2.179664007520499, 2.6766909968024932, 3.130810857] + #if (referenceTrack.frames[i].processingSoftwareVersion == '2.025' and secondaryTrack.frames[i].processingSoftwareVersion == '2.023') or \ + # (referenceTrack.frames[i].processingSoftwareVersion == '2.023' and secondaryTrack.frames[i].processingSoftwareVersion == '2.025'): + + # # changed value number of samples to estimate new value new values estimate area + # ########################################################################################################################### + # # 2.6766909968024932-->2.6581660335779866 1808694 d169-f2850, north CA + # # 2.179664007520499 -->2.204125866652153 131120 d169-f2850, north CA + + # phaseDiffFixed = [0.0, 0.4754024578084084, 0.9509913179406437, 1.4261648478671614, 2.204125866652153, 2.6581660335779866, 3.130810857] + snapThreshold = 0.2 #the above preparetions only applies to 'self._insar.modeCombination == 21' @@ -338,24 +348,36 @@ def runIonSubband(self): phaseDiffFixed = None snapThreshold = None - (phaseDiffEst, phaseDiffUsed, phaseDiffSource) = swathMosaic(referenceTrack.frames[i], inputInterferograms, self._insar.interferogram, + #whether snap for each swath + if self.swathPhaseDiffSnapIon == None: + snapSwath = [[True for jjj in range(numberOfSwaths-1)] for iii in range(numberOfFrames)] + else: + snapSwath = self.swathPhaseDiffSnapIon + if len(snapSwath) != numberOfFrames: + raise Exception('please specify each frame for parameter: swath phase difference snap to fixed values') + for iii in range(numberOfFrames): + if len(snapSwath[iii]) != (numberOfSwaths-1): + raise Exception('please specify correct number of swaths for parameter: swath phase difference snap to fixed values') + + (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = swathMosaic(referenceTrack.frames[i], inputInterferograms, self._insar.interferogram, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateFrame=False, - phaseCompensation=True, phaseDiff=phaseDiff, phaseDiffFixed=phaseDiffFixed, snapThreshold=snapThreshold, pcRangeLooks=1, pcAzimuthLooks=4, + phaseCompensation=True, phaseDiff=phaseDiff, phaseDiffFixed=phaseDiffFixed, snapThreshold=snapThreshold, snapSwath=snapSwath[i], pcRangeLooks=1, pcAzimuthLooks=4, filt=False, resamplingMethod=1) #the first item is meaningless for all the following list, so only record the following items if phaseDiff == None: phaseDiff = [None for iii in range(self._insar.startingSwath, self._insar.endingSwath + 1)] - catalog.addItem('{} subswath phase difference input'.format(ionDir['subband'][k]), phaseDiff[1:], 'runIonSubband') - catalog.addItem('{} subswath phase difference estimated'.format(ionDir['subband'][k]), phaseDiffEst[1:], 'runIonSubband') - catalog.addItem('{} subswath phase difference used'.format(ionDir['subband'][k]), phaseDiffUsed[1:], 'runIonSubband') - catalog.addItem('{} subswath phase difference used source'.format(ionDir['subband'][k]), phaseDiffSource[1:], 'runIonSubband') + catalog.addItem('frame {} {} band swath phase diff input'.format(frameNumber, ionDir['subband'][k]), phaseDiff[1:], 'runIonSubband') + catalog.addItem('frame {} {} band swath phase diff estimated'.format(frameNumber, ionDir['subband'][k]), phaseDiffEst[1:], 'runIonSubband') + catalog.addItem('frame {} {} band swath phase diff used'.format(frameNumber, ionDir['subband'][k]), phaseDiffUsed[1:], 'runIonSubband') + catalog.addItem('frame {} {} band swath phase diff used source'.format(frameNumber, ionDir['subband'][k]), phaseDiffSource[1:], 'runIonSubband') + catalog.addItem('frame {} {} band swath phase diff samples used'.format(frameNumber, ionDir['subband'][k]), numberOfValidSamples[1:], 'runIonSubband') #check if there is value around 3.130810857, which may not be stable phaseDiffUnstableExist = False for xxx in phaseDiffUsed: if abs(abs(xxx) - 3.130810857) < 0.2: phaseDiffUnstableExist = True - catalog.addItem('{} subswath phase difference unstable exists'.format(ionDir['subband'][k]), phaseDiffUnstableExist, 'runIonSubband') + catalog.addItem('frame {} {} band swath phase diff unstable exists'.format(frameNumber, ionDir['subband'][k]), phaseDiffUnstableExist, 'runIonSubband') create_xml(self._insar.amplitude, referenceTrack.frames[i].numberOfSamples, referenceTrack.frames[i].numberOfLines, 'amp') create_xml(self._insar.interferogram, referenceTrack.frames[i].numberOfSamples, referenceTrack.frames[i].numberOfLines, 'int') @@ -426,13 +448,18 @@ def runIonSubband(self): rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=False, phaseCompensation=False, resamplingMethod=0) #mosaic interferograms - frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, + (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=False, phaseCompensation=True, resamplingMethod=1) create_xml(self._insar.amplitude, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'amp') create_xml(self._insar.interferogram, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'int') + catalog.addItem('{} band frame phase diff estimated'.format(ionDir['subband'][k]), phaseDiffEst[1:], 'runIonSubband') + catalog.addItem('{} band frame phase diff used'.format(ionDir['subband'][k]), phaseDiffUsed[1:], 'runIonSubband') + catalog.addItem('{} band frame phase diff used source'.format(ionDir['subband'][k]), phaseDiffSource[1:], 'runIonSubband') + catalog.addItem('{} band frame phase diff samples used'.format(ionDir['subband'][k]), numberOfValidSamples[1:], 'runIonSubband') + #update secondary parameters here, no need to update secondary parameters here os.chdir('../') diff --git a/components/isceobj/Alos2Proc/runSwathMosaic.py b/components/isceobj/Alos2Proc/runSwathMosaic.py index 91c8719..3af99ff 100644 --- a/components/isceobj/Alos2Proc/runSwathMosaic.py +++ b/components/isceobj/Alos2Proc/runSwathMosaic.py @@ -162,7 +162,7 @@ def runSwathMosaic(self): self._insar.procDoc.addAllFromCatalog(catalog) -def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks, updateFrame=False, phaseCompensation=False, phaseDiff=None, phaseDiffFixed=None, snapThreshold=None, pcRangeLooks=1, pcAzimuthLooks=4, filt=False, resamplingMethod=0): +def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks, updateFrame=False, phaseCompensation=False, phaseDiff=None, phaseDiffFixed=None, snapThreshold=None, snapSwath=None, pcRangeLooks=1, pcAzimuthLooks=4, filt=False, resamplingMethod=0): ''' mosaic swaths @@ -181,6 +181,7 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num phaseDiff: pre-computed compensation phase for each swath phaseDiffFixed: if provided, the estimated value will snap to one of these values, which is nearest to the estimated one. snapThreshold: this is used with phaseDiffFixed + snapSwath: indicate whether snap to fixed values for each swath phase diff, must be specified if phaseDiffFixed!=None pcRangeLooks: number of range looks to take when compute swath phase difference pcAzimuthLooks: number of azimuth looks to take when compute swath phase difference filt: whether do filtering when compute swath phase difference @@ -193,6 +194,8 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num from isceobj.Alos2Proc.Alos2ProcPublic import multilook from isceobj.Alos2Proc.Alos2ProcPublic import cal_coherence_1 from isceobj.Alos2Proc.Alos2ProcPublic import filterInterferogram + from isceobj.Alos2Proc.Alos2ProcPublic import computePhaseDiff + from isceobj.Alos2Proc.Alos2ProcPublic import snap numberOfSwaths = len(frame.swaths) swaths = frame.swaths @@ -355,6 +358,8 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num # 2. 'estimated+snap': estimated from subswath overlap and snap to a fixed value # 3. 'input': pre-computed # confidence level: 3 > 2 > 1 + numberOfValidSamples = [None for i in range(numberOfSwaths)] + # only record when (filt == False) and (index[0].size >= 4000) if phaseCompensation: #compute swath phase offset diffMean = [0.0] @@ -469,48 +474,30 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num data2 *= np.exp(np.complex64(1j) * angle) print('phase offset: %15.12f rad with filter strength: %f, window size: %3d'%(diffMean0, filterStrength, filterWinSize)) else: - diffMean0 = 0.0 - for k in range(30): - dataDiff = data1 * np.conj(data2) - cor = cal_coherence_1(dataDiff, win=5) - if filt: - index = np.nonzero(np.logical_and(cor>0.95, dataDiff!=0)) - else: - index = np.nonzero(np.logical_and(cor>0.85, dataDiff!=0)) - if index[0].size < 100: - diffMean0 = 0.0 - print('\n\nWARNING: too few high coherence pixels for swath phase difference estimation') - print(' number of high coherence pixels: {}\n\n'.format(index[0].size)) - break - angle = np.mean(np.angle(dataDiff[index]), dtype=np.float64) - diffMean0 += angle - data2 *= np.exp(np.complex64(1j) * angle) - print('phase offset: %15.12f rad after loop: %3d'%(diffMean0, k)) - - DEBUG=False - if DEBUG and (k==0): - from isceobj.Alos2Proc.Alos2ProcPublic import create_xml - (length7, width7)=dataDiff.shape - filename = 'diff_ori_s{}-s{}_loop_{}.int'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, k) - dataDiff.astype(np.complex64).tofile(filename) - create_xml(filename, width7, length7, 'int') - filename = 'cor_ori_s{}-s{}_loop_{}.cor'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, k) - cor.astype(np.float32).tofile(filename) - create_xml(filename, width7, length7, 'float') + if filt: + (diffMean0, numberOfValidSamples[i]) = computePhaseDiff(data1, data2, coherenceWindowSize=5, coherenceThreshold=0.95) + else: + (diffMean0, numberOfValidSamples[i]) = computePhaseDiff(data1, data2, coherenceWindowSize=5, coherenceThreshold=0.85) + if numberOfValidSamples[i] < 100: + diffMean0 = 0.0 + print('\n\nWARNING: too few high coherence pixels for swath phase difference estimation') + print(' number of high coherence pixels: {}\n\n'.format(numberOfValidSamples[i])) + #do not record when filt + if filt: + numberOfValidSamples[i] = None + #save purely estimated diff phase phaseDiffEst[i] = diffMean0 #if fixed diff phase provided and the estimated diff phase is close enough to a fixed value, snap to it - ############################################################################################################ if phaseDiffFixed != None: - phaseDiffTmp = np.absolute(np.absolute(np.array(phaseDiffFixed)) - np.absolute(diffMean0)) - phaseDiffTmpMinIndex = np.argmin(phaseDiffTmp) - if phaseDiffTmp[phaseDiffTmpMinIndex] < snapThreshold: - diffMean0 = np.sign(diffMean0) * np.absolute(phaseDiffFixed[phaseDiffTmpMinIndex]) - phaseDiffSource[i] = 'estimated+snap' - ############################################################################################################ + if snapSwath[i-1] == True: + (outputValue, snapped) = snap(diffMean0, phaseDiffFixed, snapThreshold) + if snapped == True: + diffMean0 = outputValue + phaseDiffSource[i] = 'estimated+snap' diffMean.append(diffMean0) print('phase offset: subswath{} - subswath{}: {}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, diffMean0)) @@ -550,7 +537,7 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num if phaseCompensation: # estimated phase diff, used phase diff, used phase diff source - return (phaseDiffEst, diffMean, phaseDiffSource) + return (phaseDiffEst, diffMean, phaseDiffSource, numberOfValidSamples) def swathMosaicParameters(frame, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks): ''' diff --git a/components/isceobj/Alos2burstProc/runFrameMosaic.py b/components/isceobj/Alos2burstProc/runFrameMosaic.py index 482206f..d545ad4 100644 --- a/components/isceobj/Alos2burstProc/runFrameMosaic.py +++ b/components/isceobj/Alos2burstProc/runFrameMosaic.py @@ -102,13 +102,18 @@ def runFrameMosaic(self): rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=False, phaseCompensation=False, resamplingMethod=0) #mosaic interferograms - frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, + (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=True, phaseCompensation=True, resamplingMethod=1) create_xml(self._insar.amplitude, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'amp') create_xml(self._insar.interferogram, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'int') + catalog.addItem('frame phase diff estimated', phaseDiffEst[1:], 'runFrameMosaic') + catalog.addItem('frame phase diff used', phaseDiffUsed[1:], 'runFrameMosaic') + catalog.addItem('frame phase diff used source', phaseDiffSource[1:], 'runFrameMosaic') + catalog.addItem('frame phase diff samples used', numberOfValidSamples[1:], 'runFrameMosaic') + #update secondary parameters here #do not match for secondary, always use geometrical rangeOffsets = self._insar.frameRangeOffsetGeometricalSecondary @@ -153,11 +158,17 @@ def runFrameMosaic(self): inputSd[k].append(os.path.join('../', frameDir, 'mosaic', sdFile)) #mosaic spectral diversity interferograms - for inputSdList, outputSdFile in zip(inputSd, self._insar.interferogramSd): - frameMosaic(referenceTrack, inputSdList, outputSdFile, + for i, (inputSdList, outputSdFile) in enumerate(zip(inputSd, self._insar.interferogramSd)): + (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = frameMosaic(referenceTrack, inputSdList, outputSdFile, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=False, phaseCompensation=True, resamplingMethod=1) + catalog.addItem('sd {} frame phase diff estimated'.format(i+1), phaseDiffEst[1:], 'runFrameMosaic') + catalog.addItem('sd {} frame phase diff used'.format(i+1), phaseDiffUsed[1:], 'runFrameMosaic') + catalog.addItem('sd {} frame phase diff used source'.format(i+1), phaseDiffSource[1:], 'runFrameMosaic') + catalog.addItem('sd {} frame phase diff samples used'.format(i+1), numberOfValidSamples[1:], 'runFrameMosaic') + + for sdFile in self._insar.interferogramSd: create_xml(sdFile, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'int') diff --git a/components/isceobj/Alos2burstProc/runIonSubband.py b/components/isceobj/Alos2burstProc/runIonSubband.py index 7ae85f6..b4ecb9a 100644 --- a/components/isceobj/Alos2burstProc/runIonSubband.py +++ b/components/isceobj/Alos2burstProc/runIonSubband.py @@ -294,24 +294,36 @@ def runIonSubband(self): phaseDiffFixed = None snapThreshold = None - (phaseDiffEst, phaseDiffUsed, phaseDiffSource) = swathMosaic(referenceTrack.frames[i], inputInterferograms, self._insar.interferogram, + #whether snap for each swath + if self.swathPhaseDiffSnapIon == None: + snapSwath = [[True for jjj in range(numberOfSwaths-1)] for iii in range(numberOfFrames)] + else: + snapSwath = self.swathPhaseDiffSnapIon + if len(snapSwath) != numberOfFrames: + raise Exception('please specify each frame for parameter: swath phase difference snap to fixed values') + for iii in range(numberOfFrames): + if len(snapSwath[iii]) != (numberOfSwaths-1): + raise Exception('please specify correct number of swaths for parameter: swath phase difference snap to fixed values') + + (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = swathMosaic(referenceTrack.frames[i], inputInterferograms, self._insar.interferogram, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateFrame=False, - phaseCompensation=True, phaseDiff=phaseDiff, phaseDiffFixed=phaseDiffFixed, snapThreshold=snapThreshold, pcRangeLooks=1, pcAzimuthLooks=3, + phaseCompensation=True, phaseDiff=phaseDiff, phaseDiffFixed=phaseDiffFixed, snapThreshold=snapThreshold, snapSwath=snapSwath[i], pcRangeLooks=1, pcAzimuthLooks=3, filt=False, resamplingMethod=1) #the first item is meaningless for all the following list, so only record the following items if phaseDiff == None: phaseDiff = [None for iii in range(self._insar.startingSwath, self._insar.endingSwath + 1)] - catalog.addItem('{} subswath phase difference input'.format(ionDir['subband'][k]), phaseDiff[1:], 'runIonSubband') - catalog.addItem('{} subswath phase difference estimated'.format(ionDir['subband'][k]), phaseDiffEst[1:], 'runIonSubband') - catalog.addItem('{} subswath phase difference used'.format(ionDir['subband'][k]), phaseDiffUsed[1:], 'runIonSubband') - catalog.addItem('{} subswath phase difference used source'.format(ionDir['subband'][k]), phaseDiffSource[1:], 'runIonSubband') + catalog.addItem('frame {} {} band subswath phase diff input'.format(frameNumber, ionDir['subband'][k]), phaseDiff[1:], 'runIonSubband') + catalog.addItem('frame {} {} band subswath phase diff estimated'.format(frameNumber, ionDir['subband'][k]), phaseDiffEst[1:], 'runIonSubband') + catalog.addItem('frame {} {} band subswath phase diff used'.format(frameNumber, ionDir['subband'][k]), phaseDiffUsed[1:], 'runIonSubband') + catalog.addItem('frame {} {} band subswath phase diff used source'.format(frameNumber, ionDir['subband'][k]), phaseDiffSource[1:], 'runIonSubband') + catalog.addItem('frame {} {} band subswath phase diff samples used'.format(frameNumber, ionDir['subband'][k]), numberOfValidSamples[1:], 'runIonSubband') #check if there is value around 3.130810857, which may not be stable phaseDiffUnstableExist = False for xxx in phaseDiffUsed: if abs(abs(xxx) - 3.130810857) < 0.2: phaseDiffUnstableExist = True - catalog.addItem('{} subswath phase difference unstable exists'.format(ionDir['subband'][k]), phaseDiffUnstableExist, 'runIonSubband') + catalog.addItem('frame {} {} band subswath phase diff unstable exists'.format(frameNumber, ionDir['subband'][k]), phaseDiffUnstableExist, 'runIonSubband') create_xml(self._insar.amplitude, referenceTrack.frames[i].numberOfSamples, referenceTrack.frames[i].numberOfLines, 'amp') create_xml(self._insar.interferogram, referenceTrack.frames[i].numberOfSamples, referenceTrack.frames[i].numberOfLines, 'int') @@ -378,13 +390,18 @@ def runIonSubband(self): rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=False, phaseCompensation=False, resamplingMethod=0) #mosaic interferograms - frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, + (phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = frameMosaic(referenceTrack, inputInterferograms, self._insar.interferogram, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateTrack=False, phaseCompensation=True, resamplingMethod=1) create_xml(self._insar.amplitude, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'amp') create_xml(self._insar.interferogram, referenceTrack.numberOfSamples, referenceTrack.numberOfLines, 'int') + catalog.addItem('{} band frame phase diff estimated'.format(ionDir['subband'][k]), phaseDiffEst[1:], 'runIonSubband') + catalog.addItem('{} band frame phase diff used'.format(ionDir['subband'][k]), phaseDiffUsed[1:], 'runIonSubband') + catalog.addItem('{} band frame phase diff used source'.format(ionDir['subband'][k]), phaseDiffSource[1:], 'runIonSubband') + catalog.addItem('{} band frame phase diff samples used'.format(ionDir['subband'][k]), numberOfValidSamples[1:], 'runIonSubband') + os.chdir('../') os.chdir('../') diff --git a/components/isceobj/Sensor/MultiMode/ALOS2.py b/components/isceobj/Sensor/MultiMode/ALOS2.py index 629a0da..c302c4a 100644 --- a/components/isceobj/Sensor/MultiMode/ALOS2.py +++ b/components/isceobj/Sensor/MultiMode/ALOS2.py @@ -333,6 +333,8 @@ class ALOS2(Component): swath.rangePixelSize = Const.c/(2.0*swath.rangeSamplingRate) swath.rangeBandwidth =abs((sceneHeaderRecord.metadata['Nominal range pulse (chirp) amplitude coefficient linear term']) * (sceneHeaderRecord.metadata['Range pulse length in microsec']*1.0e-6)) + #this value is also correct + #swath.rangeBandwidth = sceneHeaderRecord.metadata['Total processor bandwidth in range'] * 1000.0 #sensingStart yr = imageData.metadata['Sensor acquisition year'] @@ -357,9 +359,16 @@ class ALOS2(Component): # '64': 'Manual observation' #print('ScanSAR mode, using PRF from the line header') swath.prf = imageData.metadata['PRF'] * 1.0e-3 + #entire azimuth spectrum is processed for ScanSAR. Here we 0.85 * minimum PRF of '08': 'ScanSAR nominal mode' (subswath 4) + swath.azimuthBandwidth = 2270.575 * 0.85 + #if operationMode == '08': + # swath.azimuthBandwidth = 2270.575 * 0.85 / 5.0 + #else: + # swath.azimuthBandwidth = 2270.575 * 0.85 / 7.0 else: #print('not ScanSAR mode, using PRF from leader file') swath.prf = sceneHeaderRecord.metadata['Pulse Repetition Frequency in mHz']*1.0e-3 + swath.azimuthBandwidth = sceneHeaderRecord.metadata['Total processor bandwidth in azimuth'] #azimuth pixel size at swath center on ground azimuthTime = swath.sensingStart + datetime.timedelta(seconds=swath.numberOfLines/swath.prf/2.0) diff --git a/components/isceobj/Sensor/MultiMode/Swath.py b/components/isceobj/Sensor/MultiMode/Swath.py index 63b5d80..2e8134f 100644 --- a/components/isceobj/Sensor/MultiMode/Swath.py +++ b/components/isceobj/Sensor/MultiMode/Swath.py @@ -91,6 +91,13 @@ AZIMUTH_PIXEL_SIZE = Component.Parameter('azimuthPixelSize', mandatory = True, doc = 'azimuth pixel size on ground in m') +AZIMUTH_BANDWIDTH = Component.Parameter('azimuthBandwidth', + public_name = 'azimuth bandwidth', + default = None, + type=float, + mandatory = True, + doc = 'azimuth bandwidth in Hz') + AZIMUTH_LINE_INTERVAL = Component.Parameter('azimuthLineInterval', public_name = 'azimuth line interval', default = None, @@ -206,6 +213,7 @@ class Swath(Component): SENSING_START, PRF, AZIMUTH_PIXEL_SIZE, + AZIMUTH_BANDWIDTH, AZIMUTH_LINE_INTERVAL, DOPPLER_VS_PIXEL, AZIMUTH_FMRATE_VS_PIXEL, diff --git a/examples/input_files/alos2/alos2App.xml b/examples/input_files/alos2/alos2App.xml index 83fc135..b5398e1 100644 --- a/examples/input_files/alos2/alos2App.xml +++ b/examples/input_files/alos2/alos2App.xml @@ -251,9 +251,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + + + + + - - + + + + - @@ -238,9 +237,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -238,9 +237,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -238,9 +237,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -238,9 +237,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -238,9 +237,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + + + + + - - + + + + + + + + - - + + + + + + + + - - + + + + + + + + - - + + + + - @@ -236,9 +235,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -236,9 +235,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -236,9 +235,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -238,9 +237,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -236,9 +235,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + + - @@ -236,9 +235,20 @@ IEEE Transactions on Geoscience and Remote Sensing, vol. 56, no. 8, pp. 4492-450 [[100, 200, 100, 200],[1000, 1200, 500, 600]] ==========================================================================================--> + + + + - - + + + +