227 lines
9.2 KiB
Python
227 lines
9.2 KiB
Python
|
#!/usr/bin/env python3
|
||
|
|
|
||
|
|
#
|
||
|
|
# Author: Cunren Liang
|
||
|
|
# Copyright 2015-present, NASA-JPL/Caltech
|
||
|
|
#
|
||
|
|
|
||
|
|
import os
|
||
|
|
import glob
|
||
|
|
import shutil
|
||
|
|
import datetime
|
||
|
|
import numpy as np
|
||
|
|
import xml.etree.ElementTree as ET
|
||
|
|
|
||
|
|
import isce, isceobj
|
||
|
|
from isceobj.Alos2Proc.Alos2ProcPublic import create_xml
|
||
|
|
from isceobj.Alos2Proc.runSwathOffset import swathOffset
|
||
|
|
from isceobj.Alos2Proc.runFrameOffset import frameOffset
|
||
|
|
from isceobj.Alos2Proc.runSwathMosaic import swathMosaic
|
||
|
|
from isceobj.Alos2Proc.runFrameMosaic import frameMosaic
|
||
|
|
|
||
|
|
from StackPulic import acquisitionModesAlos2
|
||
|
|
from StackPulic import loadTrack
|
||
|
|
|
||
|
|
def cmdLineParse():
|
||
|
|
'''
|
||
|
|
command line parser.
|
||
|
|
'''
|
||
|
|
import sys
|
||
|
|
import argparse
|
||
|
|
|
||
|
|
parser = argparse.ArgumentParser(description='form interferogram')
|
||
|
|
parser.add_argument('-ref_date_stack', dest='ref_date_stack', type=str, required=True,
|
||
|
|
help = 'reference date of stack. format: YYMMDD')
|
||
|
|
parser.add_argument('-ref_date', dest='ref_date', type=str, required=True,
|
||
|
|
help = 'reference date of this pair. format: YYMMDD')
|
||
|
|
parser.add_argument('-sec_date', dest='sec_date', type=str, required=True,
|
||
|
|
help = 'reference date of this pair. format: YYMMDD')
|
||
|
|
parser.add_argument('-nrlks1', dest='nrlks1', type=int, default=1,
|
||
|
|
help = 'number of range looks 1. default: 1')
|
||
|
|
parser.add_argument('-nalks1', dest='nalks1', type=int, default=1,
|
||
|
|
help = 'number of azimuth looks 1. default: 1')
|
||
|
|
|
||
|
|
if len(sys.argv) <= 1:
|
||
|
|
print('')
|
||
|
|
parser.print_help()
|
||
|
|
sys.exit(1)
|
||
|
|
else:
|
||
|
|
return parser.parse_args()
|
||
|
|
|
||
|
|
|
||
|
|
if __name__ == '__main__':
|
||
|
|
|
||
|
|
inps = cmdLineParse()
|
||
|
|
|
||
|
|
|
||
|
|
#get user parameters from input
|
||
|
|
dateReferenceStack = inps.ref_date_stack
|
||
|
|
dateReference = inps.ref_date
|
||
|
|
dateSecondary = inps.sec_date
|
||
|
|
numberRangeLooks1 = inps.nrlks1
|
||
|
|
numberAzimuthLooks1 = inps.nalks1
|
||
|
|
#######################################################
|
||
|
|
|
||
|
|
logFile = 'process.log'
|
||
|
|
|
||
|
|
pair = '{}-{}'.format(dateReference, dateSecondary)
|
||
|
|
|
||
|
|
ml1 = '_{}rlks_{}alks'.format(numberRangeLooks1, numberAzimuthLooks1)
|
||
|
|
|
||
|
|
interferogram = pair + ml1 + '.int'
|
||
|
|
amplitude = pair + ml1 + '.amp'
|
||
|
|
|
||
|
|
spotlightModes, stripmapModes, scansarNominalModes, scansarWideModes, scansarModes = acquisitionModesAlos2()
|
||
|
|
|
||
|
|
#use one date to find frames and swaths. any date should work, here we use dateIndexReference
|
||
|
|
frames = sorted([x[-4:] for x in glob.glob(os.path.join('./', 'f*_*'))])
|
||
|
|
swaths = sorted([int(x[-1]) for x in glob.glob(os.path.join('./', 'f1_*', 's*'))])
|
||
|
|
|
||
|
|
nframe = len(frames)
|
||
|
|
nswath = len(swaths)
|
||
|
|
|
||
|
|
trackReferenceStack = loadTrack('./', dateReferenceStack)
|
||
|
|
|
||
|
|
#mosaic swaths
|
||
|
|
for i, frameNumber in enumerate(frames):
|
||
|
|
frameDir = 'f{}_{}'.format(i+1, frameNumber)
|
||
|
|
os.chdir(frameDir)
|
||
|
|
|
||
|
|
mosaicDir = 'mosaic'
|
||
|
|
os.makedirs(mosaicDir, exist_ok=True)
|
||
|
|
os.chdir(mosaicDir)
|
||
|
|
|
||
|
|
|
||
|
|
if not (swaths[-1] - swaths[0] >= 1):
|
||
|
|
|
||
|
|
swathDir = 's{}'.format(swaths[0])
|
||
|
|
if not os.path.isfile(interferogram):
|
||
|
|
os.symlink(os.path.join('../', swathDir, interferogram), interferogram)
|
||
|
|
shutil.copy2(os.path.join('../', swathDir, interferogram+'.vrt'), interferogram+'.vrt')
|
||
|
|
shutil.copy2(os.path.join('../', swathDir, interferogram+'.xml'), interferogram+'.xml')
|
||
|
|
if not os.path.isfile(amplitude):
|
||
|
|
os.symlink(os.path.join('../', swathDir, amplitude), amplitude)
|
||
|
|
shutil.copy2(os.path.join('../', swathDir, amplitude+'.vrt'), amplitude+'.vrt')
|
||
|
|
shutil.copy2(os.path.join('../', swathDir, amplitude+'.xml'), amplitude+'.xml')
|
||
|
|
|
||
|
|
os.chdir('../../')
|
||
|
|
|
||
|
|
else:
|
||
|
|
#compute swath offset using reference stack
|
||
|
|
#geometrical offset is enough now
|
||
|
|
offsetReferenceStack = swathOffset(trackReferenceStack.frames[i], dateReferenceStack+'.slc', 'swath_offset_' + dateReferenceStack + '.txt',
|
||
|
|
crossCorrelation=False, numberOfAzimuthLooks=10)
|
||
|
|
#we can faithfully make it integer.
|
||
|
|
#this can also reduce the error due to floating point computation
|
||
|
|
rangeOffsets = [float(round(x)) for x in offsetReferenceStack[0]]
|
||
|
|
azimuthOffsets = [float(round(x)) for x in offsetReferenceStack[1]]
|
||
|
|
|
||
|
|
#list of input files
|
||
|
|
inputInterferograms = []
|
||
|
|
inputAmplitudes = []
|
||
|
|
for j, swathNumber in enumerate(range(swaths[0], swaths[-1] + 1)):
|
||
|
|
swathDir = 's{}'.format(swathNumber)
|
||
|
|
inputInterferograms.append(os.path.join('../', swathDir, interferogram))
|
||
|
|
inputAmplitudes.append(os.path.join('../', swathDir, amplitude))
|
||
|
|
|
||
|
|
#note that frame parameters do not need to be updated after mosaicking
|
||
|
|
#mosaic amplitudes
|
||
|
|
swathMosaic(trackReferenceStack.frames[i], inputAmplitudes, amplitude,
|
||
|
|
rangeOffsets, azimuthOffsets, numberRangeLooks1, numberAzimuthLooks1, resamplingMethod=0)
|
||
|
|
#mosaic interferograms
|
||
|
|
swathMosaic(trackReferenceStack.frames[i], inputInterferograms, interferogram,
|
||
|
|
rangeOffsets, azimuthOffsets, numberRangeLooks1, numberAzimuthLooks1, resamplingMethod=1)
|
||
|
|
|
||
|
|
create_xml(amplitude, trackReferenceStack.frames[i].numberOfSamples, trackReferenceStack.frames[i].numberOfLines, 'amp')
|
||
|
|
create_xml(interferogram, trackReferenceStack.frames[i].numberOfSamples, trackReferenceStack.frames[i].numberOfLines, 'int')
|
||
|
|
|
||
|
|
os.chdir('../../')
|
||
|
|
|
||
|
|
|
||
|
|
#mosaic frame
|
||
|
|
mosaicDir = 'insar'
|
||
|
|
os.makedirs(mosaicDir, exist_ok=True)
|
||
|
|
os.chdir(mosaicDir)
|
||
|
|
|
||
|
|
if nframe == 1:
|
||
|
|
frameDir = os.path.join('f1_{}/mosaic'.format(frames[0]))
|
||
|
|
if not os.path.isfile(interferogram):
|
||
|
|
os.symlink(os.path.join('../', frameDir, interferogram), interferogram)
|
||
|
|
#shutil.copy2() can overwrite
|
||
|
|
shutil.copy2(os.path.join('../', frameDir, interferogram+'.vrt'), interferogram+'.vrt')
|
||
|
|
shutil.copy2(os.path.join('../', frameDir, interferogram+'.xml'), interferogram+'.xml')
|
||
|
|
if not os.path.isfile(amplitude):
|
||
|
|
os.symlink(os.path.join('../', frameDir, amplitude), amplitude)
|
||
|
|
shutil.copy2(os.path.join('../', frameDir, amplitude+'.vrt'), amplitude+'.vrt')
|
||
|
|
shutil.copy2(os.path.join('../', frameDir, amplitude+'.xml'), amplitude+'.xml')
|
||
|
|
else:
|
||
|
|
if trackReferenceStack.operationMode in scansarModes:
|
||
|
|
matchingMode=0
|
||
|
|
else:
|
||
|
|
matchingMode=1
|
||
|
|
|
||
|
|
#geometrical offset is enough
|
||
|
|
offsetReferenceStack = frameOffset(trackReferenceStack, dateReferenceStack+'.slc', 'frame_offset_' + dateReferenceStack + '.txt',
|
||
|
|
crossCorrelation=False, matchingMode=matchingMode)
|
||
|
|
|
||
|
|
#we can faithfully make it integer.
|
||
|
|
#this can also reduce the error due to floating point computation
|
||
|
|
rangeOffsets = [float(round(x)) for x in offsetReferenceStack[0]]
|
||
|
|
azimuthOffsets = [float(round(x)) for x in offsetReferenceStack[1]]
|
||
|
|
|
||
|
|
#list of input files
|
||
|
|
inputInterferograms = []
|
||
|
|
inputAmplitudes = []
|
||
|
|
for i, frameNumber in enumerate(frames):
|
||
|
|
frameDir = 'f{}_{}'.format(i+1, frameNumber)
|
||
|
|
inputInterferograms.append(os.path.join('../', frameDir, 'mosaic', interferogram))
|
||
|
|
inputAmplitudes.append(os.path.join('../', frameDir, 'mosaic', amplitude))
|
||
|
|
|
||
|
|
#note that track parameters do not need to be updated after mosaicking
|
||
|
|
#mosaic amplitudes
|
||
|
|
frameMosaic(trackReferenceStack, inputAmplitudes, amplitude,
|
||
|
|
rangeOffsets, azimuthOffsets, numberRangeLooks1, numberAzimuthLooks1,
|
||
|
|
updateTrack=False, phaseCompensation=False, resamplingMethod=0)
|
||
|
|
#mosaic interferograms
|
||
|
|
(phaseDiffEst, phaseDiffUsed, phaseDiffSource, numberOfValidSamples) = \
|
||
|
|
frameMosaic(trackReferenceStack, inputInterferograms, interferogram,
|
||
|
|
rangeOffsets, azimuthOffsets, numberRangeLooks1, numberAzimuthLooks1,
|
||
|
|
updateTrack=False, phaseCompensation=True, resamplingMethod=1)
|
||
|
|
|
||
|
|
create_xml(amplitude, trackReferenceStack.numberOfSamples, trackReferenceStack.numberOfLines, 'amp')
|
||
|
|
create_xml(interferogram, trackReferenceStack.numberOfSamples, trackReferenceStack.numberOfLines, 'int')
|
||
|
|
|
||
|
|
#if multiple frames, remove frame amplitudes/inteferograms to save space
|
||
|
|
for x in inputAmplitudes:
|
||
|
|
os.remove(x)
|
||
|
|
os.remove(x+'.vrt')
|
||
|
|
os.remove(x+'.xml')
|
||
|
|
|
||
|
|
for x in inputInterferograms:
|
||
|
|
os.remove(x)
|
||
|
|
os.remove(x+'.vrt')
|
||
|
|
os.remove(x+'.xml')
|
||
|
|
|
||
|
|
#log output info
|
||
|
|
log = '{} at {}\n'.format(os.path.basename(__file__), datetime.datetime.now())
|
||
|
|
log += '================================================================================================\n'
|
||
|
|
log += 'frame phase diff estimated: {}\n'.format(phaseDiffEst[1:])
|
||
|
|
log += 'frame phase diff used: {}\n'.format(phaseDiffUsed[1:])
|
||
|
|
log += 'frame phase diff used source: {}\n'.format(phaseDiffSource[1:])
|
||
|
|
log += 'frame phase diff samples used: {}\n'.format(numberOfValidSamples[1:])
|
||
|
|
log += '\n'
|
||
|
|
with open(os.path.join('../', logFile), 'a') as f:
|
||
|
|
f.write(log)
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|