From d22bf1048ff0f7a077981ea4fbe8e0e6bc563961 Mon Sep 17 00:00:00 2001 From: CunrenLiang <56097947+CunrenLiang@users.noreply.github.com> Date: Fri, 24 Apr 2020 21:44:18 -0700 Subject: [PATCH] enhanced merging for ALOS-2 ionospheric correction --- components/isceobj/Alos2Proc/runIonFilt.py | 65 ++++++++--- .../isceobj/Alos2Proc/runSwathMosaic.py | 109 ++++++++++++++---- 2 files changed, 131 insertions(+), 43 deletions(-) diff --git a/components/isceobj/Alos2Proc/runIonFilt.py b/components/isceobj/Alos2Proc/runIonFilt.py index 330faf2..9c2ed5e 100644 --- a/components/isceobj/Alos2Proc/runIonFilt.py +++ b/components/isceobj/Alos2Proc/runIonFilt.py @@ -44,7 +44,8 @@ def runIonFilt(self): ################################### #SET PARAMETERS HERE #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self) - corThresholdAdj = 0.85 + corThresholdAdj = 0.97 + corOrderAdj = 20 ################################### print('\ncomputing ionosphere') @@ -75,11 +76,16 @@ def runIonFilt(self): upperUnw[area[0]:area[1], area[2]:area[3]] = 0 cor[area[0]:area[1], area[2]:area[3]] = 0 + #remove possible wired values in coherence + cor[np.nonzero(cor<0)] = 0.0 + cor[np.nonzero(cor>1)] = 0.0 + #compute ionosphere fl = SPEED_OF_LIGHT / self._insar.subbandRadarWavelength[0] fu = SPEED_OF_LIGHT / self._insar.subbandRadarWavelength[1] adjFlag = 1 - ionos = computeIonosphere(lowerUnw, upperUnw, cor, fl, fu, adjFlag, corThresholdAdj, 0) + cor[np.nonzero(cor= 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) - corThresholdIon = 0.85 + #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 = 20 + else: + #parameters for using lower/upper band coherence + corfile = subbandPrefix[0]+ml2+'.cor' + corThresholdFit = 0.4 + corOrderFit = 10 + corOrderFilt = 10 + ################################################# print('\nfiltering ionosphere') ionfile = 'ion'+ml2+'.ion' - corfile = 'diff'+ml2+'.cor' + #corfile = 'diff'+ml2+'.cor' ionfiltfile = 'filt_ion'+ml2+'.ion' img = isceobj.createImage() @@ -145,14 +174,17 @@ def runIonFilt(self): # cor[np.nonzero(wbd!=0)] = 0.00001 if fit: - ion_fit = weight_fitting(ion, cor, width, length, 1, 1, 1, 1, 2, corThresholdIon) + import copy + wgt = copy.deepcopy(cor) + wgt[np.nonzero(wgt=corThresholdAdj) + flag = (lowerUnw!=0)*(wgt!=0) index = np.nonzero(flag!=0) mv = np.mean((lowerUnw - upperUnw)[index], dtype=np.float64) print('mean value of phase difference: {}'.format(mv)) diff = mv #adjust phase using a surface else: - diff = weight_fitting(lowerUnw - upperUnw, cor, width, length, 1, 1, 1, 1, 2, corThresholdAdj) + diff = weight_fitting(lowerUnw - upperUnw, wgt, width, length, 1, 1, 1, 1, 2) flag2 = (lowerUnw!=0) index2 = np.nonzero(flag2) @@ -435,10 +466,10 @@ def cal_surface(x, y, c, order): return z -def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, coth): +def weight_fitting(ionos, weight, width, length, nrli, nali, nrlo, nalo, order): ''' ionos: input ionospheric phase - cor: coherence of the interferogram + weight: weight width: file width length: file length nrli: number of range looks of the input interferograms @@ -446,7 +477,6 @@ def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, cot nrlo: number of range looks of the output ionosphere phase nalo: number of azimuth looks of the ioutput ionosphere phase order: the order of the polynomial for fitting ionosphere phase estimates - coth: coherence threshhold for ionosphere phase estimation ''' from isceobj.Alos2Proc.Alos2ProcPublic import create_multi_index2 @@ -460,11 +490,8 @@ def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, cot rgindex = create_multi_index2(widtho, nrli, nrlo) azindex = create_multi_index2(lengtho, nali, nalo) - #convert coherence to weight - cor = cor**2/(1.009-cor**2) - #look for data to use - flag = (cor>coth)*(ionos!=0) + flag = (weight!=0)*(ionos!=0) point_index = np.nonzero(flag) m = point_index[0].shape[0] @@ -475,7 +502,7 @@ def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, cot x = x0[point_index].reshape(m, 1) y = y0[point_index].reshape(m, 1) z = ionos[point_index].reshape(m, 1) - w = cor[point_index].reshape(m, 1) + w = weight[point_index].reshape(m, 1) #convert to higher precision type before use x=np.asfarray(x,np.float64) diff --git a/components/isceobj/Alos2Proc/runSwathMosaic.py b/components/isceobj/Alos2Proc/runSwathMosaic.py index 8dad041..9d6b175 100644 --- a/components/isceobj/Alos2Proc/runSwathMosaic.py +++ b/components/isceobj/Alos2Proc/runSwathMosaic.py @@ -378,33 +378,94 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num #get difference - corth = 0.87 - if filt: - corth = 0.90 - diffMean0 = 0.0 - for k in range(5): - dataDiff = data1 * np.conj(data2) - cor = cal_coherence_1(dataDiff, win=3) - index = np.nonzero(np.logical_and(cor>corth, dataDiff!=0)) + dataDiff = data1 * np.conj(data2) + cor = cal_coherence_1(dataDiff, win=5) + index = np.nonzero(np.logical_and(cor>0.85, dataDiff!=0)) + + DEBUG=False + if DEBUG: + from isceobj.Alos2Proc.Alos2ProcPublic import create_xml + (length7, width7)=dataDiff.shape + filename = 'diff_ori_s{}-s{}.int'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber) + dataDiff.astype(np.complex64).tofile(filename) + create_xml(filename, width7, length7, 'int') + filename = 'cor_ori_s{}-s{}.cor'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber) + cor.astype(np.float32).tofile(filename) + create_xml(filename, width7, length7, 'float') + + print('\ncompute phase difference between subswaths {} and {}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber)) + print('number of pixels with coherence > 0.85: {}'.format(index[0].size)) + + #if already filtered the subswath overlap interferograms (MAI), do not filtered differential interferograms + if (filt == False) and (index[0].size < 4000): + #coherence too low, filter subswath overlap differential interferogram + diffMean0 = 0.0 + breakFlag = False + for (filterStrength, filterWinSize) in zip([3.0, 9.0], [64, 128]): + dataDiff = data1 * np.conj(data2) + dataDiff /= (np.absolute(dataDiff)+(dataDiff==0)) + dataDiff = filterInterferogram(dataDiff, filterStrength, filterWinSize, 1) + cor = cal_coherence_1(dataDiff, win=7) + + DEBUG=False + if DEBUG: + from isceobj.Alos2Proc.Alos2ProcPublic import create_xml + (length7, width7)=dataDiff.shape + filename = 'diff_filt_s{}-s{}_strength_{}_winsize_{}.int'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, filterStrength, filterWinSize) + dataDiff.astype(np.complex64).tofile(filename) + create_xml(filename, width7, length7, 'int') + filename = 'cor_filt_s{}-s{}_strength_{}_winsize_{}.cor'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, filterStrength, filterWinSize) + cor.astype(np.float32).tofile(filename) + create_xml(filename, width7, length7, 'float') + + for corth in [0.99999, 0.9999]: + index = np.nonzero(np.logical_and(cor>corth, dataDiff!=0)) + if index[0].size > 30000: + breakFlag = True + break + if breakFlag: + break + if index[0].size < 100: diffMean0 = 0.0 - print('\n\nWARNING: too few high coherence pixels for swath phase difference estimation between swath {} and {}'.format(i-1, i)) - print(' : first swath swath number: 0\n\n') - 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)) + 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)) + else: + print('filtered coherence threshold used: {}, number of pixels used: {}'.format(corth, index[0].size)) + 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 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') - DEBUG=False - if DEBUG and (k==0): - from isceobj.Alos2Proc.Alos2ProcPublic import create_xml - (lengthxx, widthxx)=dataDiff.shape - filtnamePrefix = 'subswath{}_subswath{}_loop{}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, k) - cor.astype(np.float32).tofile(filtnamePrefix+'.cor') - create_xml(filtnamePrefix+'.cor', widthxx, lengthxx, 'float') - dataDiff.astype(np.complex64).tofile(filtnamePrefix+'.int') - create_xml(filtnamePrefix+'.int', widthxx, lengthxx, 'int') diffMean.append(diffMean0) print('phase offset: subswath{} - subswath{}: {}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, diffMean0))