Skip to content

stripmapStack: fix reference date ignored in offset inversion #28129

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
huchangyang wants to merge 1 commit into
base: main
Choose a base branch
from
Open

stripmapStack: fix reference date ignored in offset inversion #28129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions contrib/stack/stripmapStack/Stack.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -474,19 +474,21 @@ def denseOffsets_Network(self, pairs, stackReference, secondaryDates, config_pre
self.runf.write(self.text_cmd + 'stripmapWrapper.py -c '+ configName+'\n')


def invertMisregPoly(self):
def invertMisregPoly(self, stackReference):

pairDirs = os.path.join(self.workDir, 'refineSecondaryTiming/pairs/')
dateDirs = os.path.join(self.workDir, 'refineSecondaryTiming/dates/')
cmd = self.text_cmd + 'invertMisreg.py -i ' + pairDirs + ' -o ' + dateDirs
cmd = self.text_cmd + 'invertMisreg.py -i ' + pairDirs + ' -o ' + dateDirs + ' -r ' + stackReference
self.runf.write(cmd + '\n')


def invertDenseOffsets(self):
def invertDenseOffsets(self, stackReference=None):

pairDirs = os.path.join(self.workDir, self.dense_offsets_folder, 'pairs')
dateDirs = os.path.join(self.workDir, self.dense_offsets_folder, 'dates')
cmd = self.text_cmd + 'invertOffsets.py -i ' + pairDirs + ' -o ' + dateDirs
if stackReference is not None:
cmd += ' -r ' + stackReference
self.runf.write(cmd + '\n')

def rubbersheet(self, secondaryDates, config_prefix):
Expand Down
30 changes: 23 additions & 7 deletions contrib/stack/stripmapStack/invertMisreg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ def createParser():
help='Directory with the overlap directories that has calculated misregistration for each pair')
parser.add_argument('-o', '--output', type=str, dest='output', required=True,
help='output directory to save misregistration for each date with respect to the stack Reference date')
parser.add_argument('-r', '--reference', type=str, dest='reference', default=None,
help='Reference date (YYYYMMDD format). If not provided, the first date in sorted dateList will be used.')
# parser.add_argument('-f', '--misregFileName', type=str, dest='misregFileName', default='misreg.txt',
# help='misreg file name that contains the calculated misregistration for a pair')

Expand Down Expand Up @@ -93,7 +95,7 @@ def getPolyInfo(filename):
#return np.size(azCoefs), np.size(rgCoefs), np.shape(azCoefs), np.shape(rgCoefs)

######################################
def design_matrix(pairDirs):
def design_matrix(pairDirs, referenceDate=None):
'''Make the design matrix for the inversion. '''
tbase,dateList,dateDict = date_list(pairDirs)
numDates = len(dateDict)
Expand Down Expand Up @@ -125,12 +127,20 @@ def design_matrix(pairDirs):
Laz[ni,:] = azOff[:]
Lrg[ni,:] = rgOff[:]

A = A[:,1:]
B = B[:,:-1]
if referenceDate is None:
refIdx = 0
elif referenceDate not in dateList:
print('Warning: Reference date {} not found in dateList. Using first date as reference.'.format(referenceDate))
refIdx = 0
else:
refIdx = dateList.index(referenceDate)

A = np.delete(A, refIdx, axis=1)
B = np.delete(B, refIdx, axis=1)

# ind=~np.isnan(Laz)
# return A[ind[:,0],:],B[ind[:,0],:],Laz[ind,:], Lrg[ind]
return A, B, Laz, Lrg
return A, B, Laz, Lrg, refIdx

######################################
def main(iargs=None):
Expand All @@ -143,7 +153,11 @@ def main(iargs=None):

tbase, dateList, dateDict = date_list(pairDirs)

A, B, Laz, Lrg = design_matrix(pairDirs)
referenceDate = inps.reference
if referenceDate is None:
referenceDate = dateList[0]

A, B, Laz, Lrg, refIdx = design_matrix(pairDirs, referenceDate=referenceDate)
A1 = np.linalg.pinv(A)
A1 = np.array(A1,np.float32)

Expand All @@ -158,8 +172,10 @@ def main(iargs=None):
RMSE_az = np.sqrt(np.sum(residual_az**2)/len(residual_az))
RMSE_rg = np.sqrt(np.sum(residual_rg**2)/len(residual_rg))

Saz = np.vstack((np.zeros((1,Saz.shape[1]), dtype=np.float32), Saz))
Srg = np.vstack((np.zeros((1,Srg.shape[1]), dtype=np.float32), Srg))
zero_row_az = np.zeros((1, Saz.shape[1]), dtype=np.float32)
zero_row_rg = np.zeros((1, Srg.shape[1]), dtype=np.float32)
Saz = np.vstack([Saz[:refIdx], zero_row_az, Saz[refIdx:]])
Srg = np.vstack([Srg[:refIdx], zero_row_rg, Srg[refIdx:]])

print('')
print('Rank of design matrix: ' + str(np.linalg.matrix_rank(A)))
Expand Down
47 changes: 35 additions & 12 deletions contrib/stack/stripmapStack/invertOffsets.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@ def createParser():
help='Directory with the pair directories that includes dense offsets for each pair')
parser.add_argument('-o', '--output', type=str, dest='output', required=True,
help='output directory to save dense-offsets for each date with respect to the stack Reference date')
parser.add_argument('-r', '--reference', type=str, dest='reference', default=None,
help='Reference date (YYYY-MM-DD HH:MM:SS format from h5 file, or YYYYMMDD format). If not provided, the first date in sorted dateList will be used.')

return parser

Expand Down Expand Up @@ -104,7 +106,7 @@ def date_list(h5file):

#####################################

def design_matrix(h5File):
def design_matrix(h5File, referenceDate=None):
tbase,dateList,dateDict, references, secondarys = date_list(h5File)
numDates = len(dateDict)
numPairs = len(references)
Expand All @@ -118,17 +120,31 @@ def design_matrix(h5File):
A[ni,ndxt2] = 1
B[ni,ndxt1:ndxt2] = tbase[ndxt1+1:ndxt2+1]-tbase[ndxt1:ndxt2]

#print('A',A)
#print('%%%%%%%%%%%%%%% %%%%%')
A = A[:,1:]
B = B[:,:-1]
if referenceDate is not None and referenceDate not in dateList:
referenceDate = datetime.datetime(
*time.strptime(referenceDate, "%Y%m%d")[0:5]
).strftime("%Y-%m-%d %H:%M:%S")

return A, B
if referenceDate is None:
refIdx = 0
elif referenceDate not in dateList:
print('Warning: Reference date {} not found in dateList. Using first date as reference.'.format(referenceDate))
refIdx = 0
else:
refIdx = dateList.index(referenceDate)

A = np.delete(A, refIdx, axis=1)
B = np.delete(B, refIdx, axis=1)

return A, B, refIdx

def invert_wlq(inps,h5File):
tbase,dateList,dateDict, references, secondarys = date_list(h5File)
referenceDate = inps.reference
if referenceDate is None:
referenceDate = dateList[0]
numPairs = len(references)
A,B = design_matrix(h5File)
A,B,refIdx = design_matrix(h5File, referenceDate=referenceDate)

h5 = h5py.File(h5File,'r')
data = h5['/platform-track/observations'].get('offset-azimuth')
Expand Down Expand Up @@ -168,9 +184,12 @@ def invert_wlq(inps,h5File):
#Cm = np.vstack((np.zeros((1,ts.shape[1]), dtype=np.float32), ts))

ts = ts.reshape([NumValidPixels,Npar]).T
#ts = np.vstack((np.zeros((1,ts.shape[1]), dtype=np.float32), ts))
ds[1:,j,ind] = ts
dsq[1:,j,ind] = Cm
zero_row = np.zeros((1, ts.shape[1]), dtype=np.float32)
ts = np.vstack([ts[:refIdx], zero_row, ts[refIdx:]])
ds[:,j,ind] = ts
Cm_zero = np.zeros((1, Cm.shape[1]), dtype=np.float32)
Cm = np.vstack([Cm[:refIdx], Cm_zero, Cm[refIdx:]])
dsq[:,j,ind] = Cm

dateListE = [d.encode("ascii", "ignore") for d in dateList]
dateListE = np.array(dateListE)
Expand Down Expand Up @@ -219,8 +238,11 @@ def invert_wlq(inps,h5File):
def invert(inps,h5File):

tbase,dateList,dateDict, references, secondarys = date_list(h5File)
referenceDate = inps.reference
if referenceDate is None:
referenceDate = dateList[0]
numPairs = len(references)
A,B = design_matrix(h5File)
A,B,refIdx = design_matrix(h5File, referenceDate=referenceDate)

h5 = h5py.File(h5File,'r')
data = h5['/platform-track/observations'].get('offset-azimuth')
Expand All @@ -245,7 +267,8 @@ def invert(inps,h5File):
#dsr[:,i,:] = L_residual
dst[i,:] = np.absolute(np.sum(np.exp(1j*L_residual),0))/Nz

ts = np.vstack((np.zeros((1,ts.shape[1]), dtype=np.float32), ts))
zero_row = np.zeros((1, ts.shape[1]), dtype=np.float32)
ts = np.vstack([ts[:refIdx], zero_row, ts[refIdx:]])
ds[:,i,:] = ts

dateListE = [d.encode("ascii", "ignore") for d in dateList]
Expand Down
4 changes: 2 additions & 2 deletions contrib/stack/stripmapStack/stackStripMap.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,7 +182,7 @@ def slcStack(inps, acquisitionDates, stackReferenceDate, secondaryDates, pairs,
i+=1
runObj = run()
runObj.configure(inps, 'run_{:02d}_invertMisreg'.format(i))
runObj.invertMisregPoly()
runObj.invertMisregPoly(stackReferenceDate)
runObj.finalize()

i+=1
Expand All @@ -203,7 +203,7 @@ def slcStack(inps, acquisitionDates, stackReferenceDate, secondaryDates, pairs,
i+=1
runObj = run()
runObj.configure(inps, 'run_{:02d}_invertDenseOffsets'.format(i))
runObj.invertDenseOffsets()
runObj.invertDenseOffsets(stackReferenceDate)
runObj.finalize()

i+=1
Expand Down