Coverage for sarvey/ifg_network.py: 70%

1#!/usr/bin/env python 

2 

3# SARvey - A multitemporal InSAR time series tool for the derivation of displacements. 

4# 

5# Copyright (C) 2021-2024 Andreas Piter (IPI Hannover, piter@ipi.uni-hannover.de) 

6# 

7# This software was developed together with FERN.Lab (fernlab@gfz-potsdam.de) in the context 

8# of the SAR4Infra project with funds of the German Federal Ministry for Digital and 

9# Transport and contributions from Landesamt fuer Vermessung und Geoinformation 

10# Schleswig-Holstein and Landesbetrieb Strassenbau und Verkehr Schleswig-Holstein. 

11# 

12# This program is free software: you can redistribute it and/or modify it under 

13# the terms of the GNU General Public License as published by the Free Software 

14# Foundation, either version 3 of the License, or (at your option) any later 

15# version. 

16# 

17# Important: This package uses PyMaxFlow. The core of PyMaxflows library is the C++ 

18# implementation by Vladimir Kolmogorov. It is also licensed under the GPL, but it REQUIRES that you 

19# cite [BOYKOV04] (see LICENSE) in any resulting publication if you use this code for research purposes. 

20# This requirement extends to SARvey. 

21# 

22# This program is distributed in the hope that it will be useful, but WITHOUT 

23# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 

24# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 

25# details. 

26# 

27# You should have received a copy of the GNU Lesser General Public License along 

28# with this program. If not, see <https://www.gnu.org/licenses/>. 

29 

30"""IfgNetwork module for SARvey.""" 

31import datetime 

32import h5py 

33import os 

34import matplotlib.pyplot as plt 

35import numpy as np 

36from typing import Union 

37import warnings 

38from logging import Logger 

39from scipy.spatial import Delaunay 

40 

41 

42class IfgNetwork: 

43 """Abstract class/interface for different types of interferogram networks.""" 

44 

45 ifg_list: Union[list, np.ndarray] = None 

46 

47 def __init__(self): 

48 """Init.""" 

49 self.pbase = None 

50 self.tbase = None 

51 self.num_images = None 

52 self.ifg_list = list() # is later converted to np.array 

53 self.pbase_ifg = None 

54 self.tbase_ifg = None 

55 self.num_ifgs = None 

56 self.dates = list() 

57 

58 def plot(self): 

59 """Plot the network of interferograms.""" 

60 fig = plt.figure(figsize=(15, 5)) 

61 axs = fig.subplots(1, 3) 

62 dt = [datetime.date.fromisoformat(d) for d in self.dates] 

63 axs[0].plot(dt, self.pbase, 'ko') 

64 for idx in self.ifg_list: 

65 xx = np.array([dt[idx[0]], dt[idx[1]]]) 

66 yy = np.array([self.pbase[idx[0]], self.pbase[idx[1]]]) 

67 axs[0].plot(xx, yy, 'k-') 

68 axs[0].set_ylabel('perpendicular baseline [m]') 

69 axs[0].set_xlabel('temporal baseline [years]') 

70 axs[0].set_title('Network of interferograms') 

71 fig.autofmt_xdate() 

72 

73 axs[1].hist(self.tbase_ifg * 365.25, bins=100) 

74 axs[1].set_ylabel('Absolute frequency') 

75 axs[1].set_xlabel('temporal baseline [days]') 

76 

77 axs[2].hist(self.pbase_ifg, bins=100) 

78 axs[2].set_ylabel('Absolute frequency') 

79 axs[2].set_xlabel('perpendicular baseline [m]') 

80 return fig 

81 

82 def getDesignMatrix(self): 

83 """Compute the design matrix for the smallbaseline network.""" 

84 a = np.zeros((self.num_ifgs, self.num_images)) 

85 for i in range(len(self.ifg_list)): 

86 a[i, self.ifg_list[i][0]] = 1 

87 a[i, self.ifg_list[i][1]] = -1 

88 return a 

89 

90 def open(self, *, path: str): 

91 """Read stored information from already existing.h5 file. 

92 

93 Parameter 

94 ----------- 

95 path: str 

96 path to existing file to read from. 

97 """ 

98 with h5py.File(path, 'r') as f: 

99 self.num_images = f.attrs["num_images"] 

100 self.num_ifgs = f.attrs["num_ifgs"] 

101 

102 self.tbase_ifg = f['tbase_ifg'][:] 

103 self.pbase_ifg = f['pbase_ifg'][:] 

104 self.tbase = f['tbase'][:] 

105 self.pbase = f['pbase'][:] 

106 self.ifg_list = f['ifg_list'][:] 

107 try: 

108 self.dates = f['dates'][:] 

109 self.dates = [date.decode("utf-8") for date in self.dates] 

110 except KeyError as ke: 

111 self.dates = None 

112 print(f"IfgNetwork is in old dataformat. Cannot read 'dates'! {ke}") 

113 

114 f.close() 

115 

116 def writeToFile(self, *, path: str, logger: Logger): 

117 """Write all existing data to .h5 file. 

118 

119 Parameters 

120 ---------- 

121 path: str 

122 path to filename 

123 logger: Logger 

124 Logging handler. 

125 """ 

126 logger.info(msg="write IfgNetwork to {}".format(path)) 

127 

128 if os.path.exists(path): 

129 os.remove(path) 

130 

131 dates = np.array(self.dates, dtype=np.string_) 

132 

133 with h5py.File(path, 'w') as f: 

134 f.attrs["num_images"] = self.num_images 

135 f.attrs["num_ifgs"] = self.num_ifgs 

136 

137 f.create_dataset('tbase_ifg', data=self.tbase_ifg) 

138 f.create_dataset('pbase_ifg', data=self.pbase_ifg) 

139 f.create_dataset('tbase', data=self.tbase) 

140 f.create_dataset('pbase', data=self.pbase) 

141 f.create_dataset('ifg_list', data=self.ifg_list) 

142 f.create_dataset('dates', data=dates) 

143 

144 

145class StarNetwork(IfgNetwork): 

146 """Star network of interferograms (single-reference).""" 

147 

148 def configure(self, *, pbase: np.ndarray, tbase: np.ndarray, ref_idx: int, dates: list): 

149 """Create list of interferograms containing the indices of the images and computes baselines. 

150 

151 Parameter 

152 --------- 

153 pbase: np.ndarray 

154 Perpendicular baselines of the SAR acquisitions. 

155 tbase: np.ndarray 

156 Temporal baselines of the SAR acquisitions. 

157 ref_idx: int 

158 Index of the reference image. 

159 dates: list 

160 Dates of the acquisitions. 

161 """ 

162 self.pbase = pbase 

163 self.tbase = tbase / 365.25 

164 self.num_images = pbase.shape[0] 

165 self.dates = dates 

166 

167 for i in range(self.num_images): 

168 if i == ref_idx: 

169 continue 

170 self.ifg_list.append((ref_idx, i)) 

171 

172 self.pbase_ifg = np.delete(self.pbase - self.pbase[ref_idx], ref_idx) 

173 self.tbase_ifg = np.delete(self.tbase - self.tbase[ref_idx], ref_idx) 

174 self.num_ifgs = self.num_images - 1 

175 

176 

177class SmallTemporalBaselinesNetwork(IfgNetwork): 

178 """Small temporal baselines network of interferograms without restrictions on the perpendicular baselines.""" 

179 

180 def configure(self, *, pbase: np.ndarray, tbase: np.ndarray, num_link: int = None, dates: list): 

181 """Create list of interferograms containing the indices of the images and computes baselines. 

182 

183 Parameter 

184 ----------- 

185 pbase: np.ndarray 

186 Perpendicular baselines of the SAR acquisitions. 

187 tbase: np.ndarray 

188 Temporal baselines of the SAR acquisitions. 

189 num_link: int 

190 Number of consecutive links in time connecting acquisitions. 

191 dates: list 

192 Dates of the acquisitions. 

193 """ 

194 self.pbase = pbase 

195 self.tbase = tbase / 365.25 

196 self.num_images = pbase.shape[0] 

197 self.dates = dates 

198 

199 for i in range(self.num_images): 

200 for j in range(num_link): 

201 if i + j + 1 >= self.num_images: 

202 continue 

203 self.ifg_list.append((i, i + j + 1)) 

204 

205 self.ifg_list = [(i, j) for i, j in self.ifg_list if i != j] # remove connections to itself, e.g. (0, 0) 

206 

207 self.pbase_ifg = np.array([self.pbase[idx[1]] - self.pbase[idx[0]] for idx in self.ifg_list]) 

208 self.tbase_ifg = np.array([self.tbase[idx[1]] - self.tbase[idx[0]] for idx in self.ifg_list]) 

209 self.num_ifgs = self.pbase_ifg.shape[0] 

210 

211 

212class SmallBaselineNetwork(IfgNetwork): 

213 """Small baseline network of interferograms restricting both temporal and spatial baselines.""" 

214 

215 def configure(self, *, pbase: np.ndarray, tbase: np.ndarray, num_link: int, max_tbase: int, dates: list): 

216 """Create list of interferograms containing the indices of the images and computes baselines. 

217 

218 Parameter 

219 ----------- 

220 pbase: np.ndarray 

221 perpendicular baselines of the SAR acquisitions. 

222 tbase: np.ndarray 

223 temporal baselines of the SAR acquisitions. 

224 max_tbase: int 

225 maximum temporal baseline in [days] (default: None). 

226 num_link: int 

227 number of links within the range of maximum temporal baseline. 

228 dates: list 

229 Dates of the acquisitions. 

230 """ 

231 self.pbase = pbase 

232 self.tbase = tbase / 365.25 

233 self.num_images = pbase.shape[0] 

234 self.dates = dates 

235 flag_restrict_to_max_tbase = False 

236 

237 # in this section use tbase in [days] (function argument, not self.) 

238 for i in range(self.num_images - 1): 

239 

240 if i + 1 < self.num_images - 1: 

241 # always use one connection to nearest neighbour in time 

242 self.ifg_list.append((i, i + 1)) 

243 else: 

244 self.ifg_list.append((i, i + 1)) 

245 break 

246 # compute index corresponding to max_tbase for current time 

247 diff = np.abs(tbase - (tbase[i] + max_tbase)) 

248 max_idx = np.where(diff == diff.min())[0][0] 

249 self.ifg_list.append((i, max_idx)) 

250 

251 if max_idx == i: # no further images between i and max_idx 

252 flag_restrict_to_max_tbase = True 

253 continue 

254 

255 # spread the rest of the links over the remaining time steps in between 

256 links = np.floor(np.arange(i, max_idx, (max_idx - i) / (num_link - 1)))[1:].astype(int) 

257 for link in links: 

258 self.ifg_list.append((i, link)) 

259 self.ifg_list = np.unique(self.ifg_list, axis=0) 

260 

261 if flag_restrict_to_max_tbase: 

262 warnings.warn(f"Cannot restrict ifgs to maximum temporal baseline of {max_tbase} days.") 

263 

264 self.ifg_list = [(i, j) for i, j in self.ifg_list if i != j] # remove connections to itself, e.g. (0, 0) 

265 

266 self.pbase_ifg = np.array([self.pbase[idx[1]] - self.pbase[idx[0]] for idx in self.ifg_list]) 

267 self.tbase_ifg = np.array([self.tbase[idx[1]] - self.tbase[idx[0]] for idx in self.ifg_list]) 

268 self.num_ifgs = self.pbase_ifg.shape[0] 

269 

270 

271class DelaunayNetwork(IfgNetwork): 

272 """Delaunay network of interferograms which restricts both the temporal and perpendicular baselines.""" 

273 

274 def configure(self, *, pbase: np.ndarray, tbase: np.ndarray, dates: list): 

275 """Create list of interferograms containing the indices of the images and computes baselines. 

276 

277 Parameter 

278 ----------- 

279 pbase: np.ndarray 

280 perpendicular baselines of the SAR acquisitions, array 

281 tbase: np.ndarray 

282 temporal baselines of the SAR acquisitions, array 

283 dates: list 

284 Dates of the acquisitions, list. 

285 """ 

286 self.pbase = pbase 

287 self.tbase = tbase / 365.25 

288 self.num_images = pbase.shape[0] 

289 self.dates = dates 

290 scale = 0.25 

291 

292 network = Delaunay(points=np.stack([self.pbase, self.tbase * 365.25 * scale]).T) 

293 for p1, p2, p3 in network.simplices: 

294 self.ifg_list.append((p1, p2)) 

295 self.ifg_list.append((p1, p3)) 

296 self.ifg_list.append((p2, p3)) 

297 self.ifg_list = np.unique(self.ifg_list, axis=0) 

298 

299 self.pbase_ifg = np.array([self.pbase[idx[1]] - self.pbase[idx[0]] for idx in self.ifg_list]) 

300 self.tbase_ifg = np.array([self.tbase[idx[1]] - self.tbase[idx[0]] for idx in self.ifg_list]) 

301 self.num_ifgs = self.pbase_ifg.shape[0] 

302 

303 

304class SmallBaselineYearlyNetwork(IfgNetwork): 

305 """Small baseline network of interferograms with yearly connections.""" 

306 

307 def configure(self, *, pbase: np.ndarray, tbase: np.ndarray, num_link: int = None, dates: list): 

308 """Create list of interferograms containing the indices of the images and computes baselines. 

309 

310 Parameter 

311 ----------- 

312 pbase: np.ndarray 

313 perpendicular baselines of the SAR acquisitions, array 

314 tbase: np.ndarray 

315 temporal baselines of the SAR acquisitions, array 

316 num_link: int 

317 Number of consecutive links in time connecting acquisitions. 

318 dates: list 

319 Dates of the acquisitions, list. 

320 """ 

321 self.pbase = pbase 

322 self.tbase = tbase / 365.25 

323 self.num_images = pbase.shape[0] 

324 self.dates = dates 

325 

326 # add small temporal baselines 

327 for i in range(self.num_images): 

328 for j in range(num_link): 

329 if i + j + 1 >= self.num_images: 

330 continue 

331 self.ifg_list.append((i, i + j + 1)) 

332 

333 # add yearly ifgs 

334 for i in range(self.num_images): 

335 # find index of image at roughly one year distance 

336 diff = np.abs(tbase - (tbase[i] + 365.25)) 

337 year_idx = np.where(diff == diff.min())[0][0] 

338 print(year_idx) 

339 if year_idx != self.num_images - 1: # avoid connections to the last image 

340 self.ifg_list.append((i, year_idx)) 

341 print("found!") 

342 

343 self.ifg_list = np.unique(self.ifg_list, axis=0) 

344 self.ifg_list = [(i, j) for i, j in self.ifg_list if i != j] # remove connections to itself, e.g. (0, 0) 

345 

346 self.pbase_ifg = np.array([self.pbase[idx[1]] - self.pbase[idx[0]] for idx in self.ifg_list]) 

347 self.tbase_ifg = np.array([self.tbase[idx[1]] - self.tbase[idx[0]] for idx in self.ifg_list]) 

348 self.num_ifgs = self.pbase_ifg.shape[0]