Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
#!/usr/bin/env python3
#
# Author: Cunren Liang
# Copyright 2021
#
import os
import glob
import shutil
import datetime
import numpy as np
import isce, isceobj
from isceobj.Alos2Proc.Alos2ProcPublic import create_xml
def datesFromPairs(pairs):
'''get all dates from pairs
'''
dates = []
for p in pairs:
for x in p.split('_'):
if x not in dates:
dates.append(x)
dates.sort()
return dates
def least_sqares(H, S, W=None):
'''
#This can make use multiple threads (set environment variable: OMP_NUM_THREADS)
linear equations: H theta = s
W: weight matrix
'''
S.reshape(H.shape[0], 1)
if W is None:
#use np.dot instead since some old python versions don't have matmul
m1 = np.linalg.inv(np.dot(H.transpose(), H))
Z = np.dot( np.dot(m1, H.transpose()) , S)
else:
#use np.dot instead since some old python versions don't have matmul
m1 = np.linalg.inv(np.dot(np.dot(H.transpose(), W), H))
Z = np.dot(np.dot(np.dot(m1, H.transpose()), W), S)
return Z.reshape(Z.size)
def interp_2d(data1, numberRangeLooks1, numberRangeLooks2, numberAzimuthLooks1, numberAzimuthLooks2, width2=None, length2=None):
'''
interpolate data1 of numberRangeLooks1/numberAzimuthLooks1 to data2 of numberRangeLooks2/numberAzimuthLooks2
'''
length1, width1 = data1.shape
if width2 is None:
width2 = int(np.around(width1*numberRangeLooks1/numberRangeLooks2))
if length2 is None:
length2 = int(np.around(length1*numberAzimuthLooks1/numberAzimuthLooks2))
#number of range looks input
nrli = numberRangeLooks1
#number of range looks output
nrlo = numberRangeLooks2
#number of azimuth looks input
nali = numberAzimuthLooks1
#number of azimuth looks output
nalo = numberAzimuthLooks2
index1 = np.linspace(0, width1-1, num=width1, endpoint=True)
index2 = np.linspace(0, width2-1, num=width2, endpoint=True) * nrlo/nrli + (nrlo-nrli)/(2.0*nrli)
data2 = np.zeros((length2, width2), dtype=data1.dtype)
for i in range(length1):
f = interp1d(index1, data1[i,:], kind='cubic', fill_value="extrapolate")
data2[i, :] = f(index2)
index1 = np.linspace(0, length1-1, num=length1, endpoint=True)
index2 = np.linspace(0, length2-1, num=length2, endpoint=True) * nalo/nali + (nalo-nali)/(2.0*nali)
for j in range(width2):
f = interp1d(index1, data2[0:length1, j], kind='cubic', fill_value="extrapolate")
data2[:, j] = f(index2)
return data2
def cmdLineParse():
'''
command line parser.
'''
import sys
import argparse
parser = argparse.ArgumentParser(description='least squares estimation')
parser.add_argument('--idir', dest='idir', type=str, required=True,
help = 'input directory where each pair (YYYYMMDD_YYYYMMDD) is located. only folders are recognized')
parser.add_argument('--odir', dest='odir', type=str, required=True,
help = 'output directory for estimated result of each date')
parser.add_argument('--zro_date', dest='zro_date', type=str, default=None,
help = 'date in least squares estimation whose ionospheric phase is assumed to be zero. format: YYYYMMDD. default: first date')
parser.add_argument('--exc_date', dest='exc_date', type=str, nargs='+', default=[],
help = 'pairs involving these dates are excluded in least squares estimation. a number of dates seperated by blanks. format: YYYYMMDD YYYYMMDD YYYYMMDD...')
parser.add_argument('--exc_pair', dest='exc_pair', type=str, nargs='+', default=[],
help = 'pairs excluded in least squares estimation. a number of pairs seperated by blanks. format: YYYYMMDD-YYYYMMDD YYYYMMDD-YYYYMMDD...')
parser.add_argument('--tsmax', dest='tsmax', type=float, default=None,
help = 'maximum time span in years of pairs used in least squares estimation. default: None')
parser.add_argument('--nrlks1', dest='nrlks1', type=int, default=1,
help = 'number of range looks of input. default: 1')
parser.add_argument('--nalks1', dest='nalks1', type=int, default=1,
help = 'number of azimuth looks of input. default: 1')
parser.add_argument('--nrlks2', dest='nrlks2', type=int, default=1,
help = 'number of range looks of output. default: 1')
parser.add_argument('--nalks2', dest='nalks2', type=int, default=1,
help = 'number of azimuth looks of output. default: 1')
parser.add_argument('--width2', dest='width2', type=int, default=None,
help = 'width of output result. default: None, determined by program')
parser.add_argument('--length2', dest='length2', type=int, default=None,
help = 'length of output result. default: None, determined by program')
parser.add_argument('--merged_geom', dest='merged_geom', type=str, default=None,
help = 'a merged geometry file for getting width2/length2, e.g. merged/geom_reference/hgt.rdr. if provided, --width2/--length2 will be overwritten')
parser.add_argument('--interp', dest='interp', action='store_true', default=False,
help='interpolate estimated result to nrlks2/nalks2 sample size')
parser.add_argument('--msk_overlap', dest='msk_overlap', action='store_true', default=False,
help='mask output with overlap of all acquisitions')
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
idir = inps.idir
odir = inps.odir
dateZero = inps.zro_date
dateExcluded = inps.exc_date
pairExcluded = inps.exc_pair
tsmax = inps.tsmax
numberRangeLooks1 = inps.nrlks1
numberAzimuthLooks1 = inps.nalks1
numberRangeLooks2 = inps.nrlks2
numberAzimuthLooks2 = inps.nalks2
width2 = inps.width2
length2 = inps.length2
mergedGeom = inps.merged_geom
interp = inps.interp
maskOverlap = inps.msk_overlap
#######################################################
#all pair folders in order
pairDirs = sorted(glob.glob(os.path.join(os.path.abspath(idir), '*_*')))
pairDirs = [x for x in pairDirs if os.path.isdir(x)]
#all pairs in order
pairsAll = [os.path.basename(x) for x in pairDirs]
#all dates in order
datesAll = datesFromPairs(pairsAll)
#select pairs
pairs = []
for x in pairsAll:
dateReference, dateSecondary = x.split('_')
timeReference = datetime.datetime.strptime(dateReference, "%Y%m%d")
timeSecondary = datetime.datetime.strptime(dateSecondary, "%Y%m%d")
ts = np.absolute((timeSecondary - timeReference).total_seconds()) / (365.0 * 24.0 * 3600)
if (dateReference in dateExcluded) and (dateSecondary in dateExcluded):
continue
if (x in pairExcluded):
continue
if tsmax is not None:
if ts > tsmax:
continue
pairs.append(x)
dates = datesFromPairs(pairs)
if dateZero is not None:
if dateZero not in dates:
raise Exception('zro_date provided by user not in the dates involved in least squares estimation.')
else:
dateZero = dates[0]
print(f'all pairs ({len(pairsAll)}):\n{pairsAll}')
print(f'all dates ({len(datesAll)}):\n{datesAll}')
print(f'used pairs ({len(pairs)}):\n{pairs}')
print(f'used dates ({len(dates)}):\n{dates}')
####################################################################################
print('\nSTEP 1. read files')
####################################################################################
ndate = len(dates)
npair = len(pairs)
ionfile = os.path.join(idir, pairs[0], 'ion_cal', 'filt.ion')
img = isceobj.createImage()
img.load(ionfile+'.xml')
width = img.width
length = img.length
ionPairs = np.zeros((npair, length, width), dtype=np.float32)
flag = np.ones((length, width), dtype=np.float32)
#this is reserved for use
wls = False
stdPairs = np.ones((npair, length, width), dtype=np.float32)
for i in range(npair):
ionfile = os.path.join(idir, pairs[i], 'ion_cal', 'filt.ion')
ionPairs[i, :, :] = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
#flag of valid/invalid is defined by amplitde image
amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
flag *= (amp!=0)
####################################################################################
print('\nSTEP 2. do least squares')
####################################################################################
import copy
from numpy.linalg import matrix_rank
dates2 = copy.deepcopy(dates)
dates2.remove(dateZero)
#observation matrix
H0 = np.zeros((npair, ndate-1))
for k in range(npair):
dateReference = pairs[k].split('_')[0]
dateSecondary = pairs[k].split('_')[1]
if dateReference != dateZero:
dateReference_i = dates2.index(dateReference)
H0[k, dateReference_i] = 1
if dateSecondary != dateZero:
dateSecondary_i = dates2.index(dateSecondary)
H0[k, dateSecondary_i] = -1
rank = matrix_rank(H0)
if rank < ndate-1:
raise Exception('dates to be estimated are not fully connected by the pairs used in least squares')
else:
print('number of pairs to be used in least squares: {}'.format(npair))
print('number of dates to be estimated: {}'.format(ndate-1))
print('observation matrix rank: {}'.format(rank))
ts = np.zeros((ndate-1, length, width), dtype=np.float32)
for i in range(length):
if (i+1) % 50 == 0 or (i+1) == length:
print('processing line: %6d of %6d' % (i+1, length), end='\r')
if (i+1) == length:
print()
for j in range(width):
#observed signal
S0 = ionPairs[:, i, j]
if wls == False:
#observed signal
S = S0
H = H0
else:
#add weight
#https://stackoverflow.com/questions/19624997/understanding-scipys-least-square-function-with-irls
#https://stackoverflow.com/questions/27128688/how-to-use-least-squares-with-weight-matrix-in-python
wgt = (stdPairs[:, i, j])**2
W = np.sqrt(1.0/wgt)
H = H0 * W[:, None]
S = S0 * W
#do least-squares estimation
#[theta, residuals, rank, singular] = np.linalg.lstsq(H, S)
#make W full matrix if use W here (which is a slower method)
#'using W before this' is faster
theta = least_sqares(H, S, W=None)
ts[:, i, j] = theta
# #dump raw estimate
# cdir = os.getcwd()
# os.makedirs(odir, exist_ok=True)
# os.chdir(odir)
# for i in range(ndate-1):
# file_name = 'filt_ion_'+dates2[i]+ml2+'.ion'
# ts[i, :, :].astype(np.float32).tofile(file_name)
# create_xml(file_name, width, length, 'float')
# file_name = 'filt_ion_'+dateZero+ml2+'.ion'
# (np.zeros((length, width), dtype=np.float32)).astype(np.float32).tofile(file_name)
# create_xml(file_name, width, length, 'float')
# os.chdir(cdir)
####################################################################################
print('\nSTEP 3. interpolate ionospheric phase')
####################################################################################
from scipy.interpolate import interp1d
width1 = width
length1 = length
if width2 is None:
width2 = int(width1 * numberRangeLooks1 / numberRangeLooks2)
if length2 is None:
length2 = int(length1 * numberAzimuthLooks1 / numberAzimuthLooks2)
if mergedGeom is not None:
from osgeo import gdal
ds = gdal.Open(mergedGeom + ".vrt", gdal.GA_ReadOnly)
width2 = ds.RasterXSize
length2 = ds.RasterYSize
os.makedirs(odir, exist_ok=True)
for idate in range(ndate-1):
print('interplate {}'.format(dates2[idate]))
ionrectfile = os.path.join(odir, dates2[idate]+'.ion')
if interp and ((numberRangeLooks1 != numberRangeLooks2) or (numberAzimuthLooks1 != numberAzimuthLooks2)):
ionrect = interp_2d(ts[idate, :, :], numberRangeLooks1, numberRangeLooks2, numberAzimuthLooks1, numberAzimuthLooks2,
width2=width2, length2=length2)
#mask with overlap of all acquistions
if maskOverlap:
if idate == 0:
flagrect = interp_2d(flag, numberRangeLooks1, numberRangeLooks2, numberAzimuthLooks1, numberAzimuthLooks2,
width2=width2, length2=length2)
ionrect *= (flagrect>0.5)
ionrect.astype(np.float32).tofile(ionrectfile)
create_xml(ionrectfile, width2, length2, 'float')
else:
ionrect = ts[idate, :, :]
if maskOverlap:
ionrect *= flag
ionrect.astype(np.float32).tofile(ionrectfile)
create_xml(ionrectfile, width1, length1, 'float')
ionrectfile = os.path.join(odir, dateZero+'.ion')
if interp and ((numberRangeLooks1 != numberRangeLooks2) or (numberAzimuthLooks1 != numberAzimuthLooks2)):
(np.zeros((length2, width2), dtype=np.float32)).astype(np.float32).tofile(ionrectfile)
create_xml(ionrectfile, width2, length2, 'float')
else:
(np.zeros((length1, width1), dtype=np.float32)).astype(np.float32).tofile(ionrectfile)
create_xml(ionrectfile, width1, length1, 'float')