Coverage for sarvey/filtering.py: 55%

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"""Filtering module for SARvey.""" 

31import time 

32import multiprocessing 

33import matplotlib.pyplot as plt 

34import numpy as np 

35from scipy.interpolate import griddata 

36import gstools as gs 

37from logging import Logger 

38 

39from mintpy.utils import ptime 

40 

41import sarvey.utils as ut 

42 

43 

44def launchSpatialFiltering(parameters: tuple): 

45 """Launch_spatial_filtering. 

46 

47 Launches the spatial filtering to estimate the atmospheric phase screen with low-pass filtering. 

48 

49 Parameters 

50 ---------- 

51 parameters: tuple 

52 Tuple containing the following parameters: 

53 

54 idx_range: np.ndarray 

55 range of indices for the time series 

56 num_time: int 

57 number of time steps 

58 residuals: np.ndarray 

59 residual phase (size: num_points x num_ifgs) 

60 coord_utm1: np.ndarray 

61 coordinates in UTM of the first-order points for which the residuals are given (size: num_points_p1 x 2) 

62 coord_utm2: np.ndarray 

63 coordinates in UTM of the new points which shall be interpolated (size: num_points_p2 x 2) 

64 bins: np.ndarray 

65 bin edges for the variogram 

66 bool_plot: bool 

67 boolean flag to plot intermediate results 

68 logger: Logger 

69 Logging handler 

70 

71 Returns 

72 ------- 

73 idx_range: np.ndarray 

74 range of indices for the time series 

75 aps1: np.ndarray 

76 atmospheric phase screen for the known points (size: num_points_p1 x num_ifgs) 

77 aps2: np.ndarray 

78 atmospheric phase screen for the new points (size: num_points_p2 x num_ifgs) 

79 """ 

80 # Unpack the parameters 

81 (idx_range, num_time, residuals, coord_utm1, coord_utm2, bins, bool_plot, logger) = parameters 

82 

83 x = coord_utm1[:, 1] 

84 y = coord_utm1[:, 0] 

85 x_new = coord_utm2[:, 1] 

86 y_new = coord_utm2[:, 0] 

87 

88 aps1 = np.zeros((coord_utm1.shape[0], num_time), dtype=np.float32) 

89 aps2 = np.zeros((coord_utm2.shape[0], num_time), dtype=np.float32) 

90 

91 prog_bar = ptime.progressBar(maxValue=num_time) 

92 

93 for i in range(num_time): 

94 field = residuals[:, i].astype(np.float32) 

95 

96 # 1) estimate the variogram of the field 

97 bin_center, vario = gs.vario_estimate(pos=[x, y], field=field, bin_edges=bins) 

98 

99 # 2) fit model to empirical variogram 

100 model = gs.Stable(dim=2) 

101 try: 

102 model.fit_variogram(x_data=bin_center, y_data=vario, nugget=True, max_eval=1500) 

103 except RuntimeError as err: 

104 logger.error(msg="\nIMAGE {}: Not able to fit variogram! {}".format(idx_range[i], err)) 

105 if bool_plot: 

106 fig, ax = plt.subplots(2, figsize=[10, 5]) 

107 sca1 = ax[0].scatter(x, y, c=field) 

108 plt.colorbar(sca1, ax=ax[0], pad=0.03, shrink=0.5) 

109 ax[0].set_title("Not able to fit variogram! - PS1 residuals") 

110 ax[1].scatter(bin_center, vario) 

111 ax[1].set_xlabel("distance in [m]") 

112 ax[1].set_ylabel("semi-variogram") 

113 plt.close(fig) 

114 prog_bar.update(value=i + 1, every=1, suffix='{}/{} images'.format(i + 1, num_time)) 

115 continue 

116 

117 # 3) estimate parameters of kriging 

118 sk = gs.krige.Simple( 

119 model=model, 

120 cond_pos=[x, y], 

121 cond_val=field, 

122 ) 

123 

124 # 4) evaluate the kriging model at ORIGINAL locations 

125 fld_sk, _ = sk((x, y), return_var=True) 

126 aps1[:, i] = fld_sk 

127 

128 # 5) evaluate the kriging model at NEW locations 

129 fld_sk_new, var_sk_new = sk((x_new, y_new), return_var=True) 

130 aps2[:, i] = fld_sk_new 

131 

132 prog_bar.update(value=i + 1, every=1, suffix='{}/{} images'.format(i + 1, num_time)) 

133 

134 # 5) show results 

135 if bool_plot: 

136 min_val = np.min(field) 

137 max_val = np.max(field) 

138 

139 fig, ax = plt.subplots(2, 2, figsize=[10, 5]) 

140 

141 cur_ax = ax[0, 0] 

142 sca1 = cur_ax.scatter(x, y, c=field, vmin=min_val, vmax=max_val) 

143 plt.colorbar(sca1, ax=cur_ax, pad=0.03, shrink=0.5) 

144 cur_ax.set_title("PS1 residuals") 

145 

146 cur_ax = ax[0, 1] 

147 cur_ax = model.plot(x_max=bin_center[-1], ax=cur_ax) 

148 cur_ax.scatter(bin_center, vario) 

149 cur_ax.set_xlabel("distance in [m]") 

150 cur_ax.set_ylabel("semi-variogram") 

151 

152 if coord_utm2 is not None: 

153 cur_ax = ax[1, 0] 

154 sca2 = cur_ax.scatter(x_new, y_new, c=fld_sk_new, vmin=min_val, vmax=max_val) 

155 plt.colorbar(sca2, ax=cur_ax, pad=0.03, shrink=0.5) 

156 cur_ax.set_title("PS2 prediction of atmospheric effect") 

157 

158 cur_ax = ax[0, 1] 

159 sca4 = cur_ax.scatter(x_new, y_new, c=var_sk_new) 

160 plt.colorbar(sca4, ax=cur_ax, pad=0.03, shrink=0.5) 

161 cur_ax.set_title("Variance of predicted atmospheric effect") 

162 

163 plt.close(fig) 

164 

165 return idx_range, aps1, aps2 

166 

167 

168def estimateAtmosphericPhaseScreen(*, residuals: np.ndarray, coord_utm1: np.ndarray, coord_utm2: np.ndarray, 

169 num_cores: int = 1, bool_plot: bool = False, 

170 logger: Logger) -> tuple[np.ndarray, np.ndarray]: 

171 """Estimate_atmospheric_phase_screen. 

172 

173 Estimates the atmospheric phase screen from a stack of phase time series for a sparse set of points. 

174 Kriging is used to estimate the spatial dependence and to interpolate the phase screen over a set of new points. 

175 

176 Parameters 

177 ---------- 

178 residuals: np.ndarray 

179 residual phase (size: num_points1 x num_images) 

180 coord_utm1: np.ndarray 

181 coordinates in UTM of the points for which the residuals are given (size: num_points1 x 2) 

182 coord_utm2: np.ndarray 

183 coordinates in UTM of the new points which shall be interpolated (size: num_points2 x 2) 

184 num_cores: int 

185 Number of cores 

186 bool_plot: bool 

187 boolean flag to plot intermediate results (default: False) 

188 logger: Logger 

189 Logging handler 

190 

191 Returns 

192 ------- 

193 aps1: np.ndarray 

194 atmospheric phase screen for the known points (size: num_points1 x num_images) 

195 aps2: np.ndarray 

196 atmospheric phase screen for the new points (size: num_points2 x num_images) 

197 """ 

198 msg = "#" * 10 

199 msg += " ESTIMATE ATMOSPHERIC PHASE SCREEN (KRIGING) " 

200 msg += "#" * 10 

201 logger.info(msg=msg) 

202 

203 start_time = time.time() 

204 

205 num_points1 = residuals.shape[0] 

206 num_points2 = coord_utm2.shape[0] 

207 num_time = residuals.shape[1] # can be either num_ifgs or num_images 

208 

209 bins = gs.variogram.standard_bins(pos=(coord_utm1[:, 1], coord_utm1[:, 0]), 

210 dim=2, latlon=False, mesh_type='unstructured', bin_no=30, max_dist=None) 

211 

212 if num_cores == 1: 

213 args = (np.arange(0, num_time), num_time, residuals, coord_utm1, coord_utm2, bins, bool_plot, logger) 

214 _, aps1, aps2 = launchSpatialFiltering(parameters=args) 

215 else: 

216 logger.info(msg="start parallel processing with {} cores.".format(num_cores)) 

217 pool = multiprocessing.Pool(processes=num_cores) 

218 

219 aps1 = np.zeros((num_points1, num_time), dtype=np.float32) 

220 aps2 = np.zeros((num_points2, num_time), dtype=np.float32) 

221 

222 num_cores = num_time if num_cores > num_time else num_cores # avoids having more samples than cores 

223 idx = ut.splitDatasetForParallelProcessing(num_samples=num_time, num_cores=num_cores) 

224 

225 args = [( 

226 idx_range, 

227 idx_range.shape[0], 

228 residuals[:, idx_range], 

229 coord_utm1, 

230 coord_utm2, 

231 bins, 

232 False, 

233 logger) for idx_range in idx] 

234 

235 results = pool.map(func=launchSpatialFiltering, iterable=args) 

236 

237 # retrieve results 

238 for i, aps1_i, aps2_i in results: 

239 aps1[:, i] = aps1_i 

240 aps2[:, i] = aps2_i 

241 

242 m, s = divmod(time.time() - start_time, 60) 

243 logger.debug(msg='time used: {:02.0f} mins {:02.1f} secs.\n'.format(m, s)) 

244 

245 return aps1, aps2 

246 

247 

248def simpleInterpolation(*, residuals: np.ndarray, coord_utm1: np.ndarray, coord_utm2: np.ndarray, 

249 interp_method: str = "linear"): 

250 """SimpleInterpolation. 

251 

252 Simple interpolation of atmospheric phase screen using scipy's griddata function with options "linear" or "cubic". 

253 For pixels outside the convex hull of the input points, the nearest neighbor is used. 

254 

255 Parameters 

256 ---------- 

257 residuals: np.ndarray 

258 residual phase (size: num_points x num_ifgs) 

259 coord_utm1: np.ndarray 

260 coordinates in UTM of the points for which the residuals are given (size: num_points_p1 x 2) 

261 coord_utm2: np.ndarray 

262 coordinates in UTM of the new points which shall be interpolated (size: num_points_p2 x 2) 

263 interp_method: str 

264 interpolation method (default: "linear"; options: "linear", "cubic") 

265 

266 Returns 

267 ------- 

268 aps1: np.ndarray 

269 atmospheric phase screen for the known points (size: num_points_p1 x num_images) 

270 aps2: np.ndarray 

271 atmospheric phase screen for the new points (size: num_points_p2 x num_images) 

272 """ 

273 num_points2 = coord_utm2.shape[0] 

274 num_images = residuals.shape[1] 

275 

276 aps1 = np.zeros_like(residuals, dtype=np.float32) 

277 aps2 = np.zeros((num_points2, num_images), dtype=np.float32) 

278 for i in range(num_images): 

279 aps1[:, i] = griddata(coord_utm1, residuals[:, i], coord_utm1, method=interp_method) 

280 aps2[:, i] = griddata(coord_utm1, residuals[:, i], coord_utm2, method=interp_method) 

281 # interpolation with 'linear' or 'cubic' yields nan values for pixel that need to be extrapolated. 

282 # interpolation with 'knn' solves this problem. 

283 mask_extrapolate = np.isnan(aps2[:, i]) 

284 aps2[mask_extrapolate, i] = griddata( 

285 coord_utm1, 

286 residuals[:, i], 

287 coord_utm2[mask_extrapolate, :], 

288 method='nearest' 

289 ) 

290 

291 return aps1, aps2