Coverage for sarvey/unwrapping.py: 76%

1#!/usr/bin/env python

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

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/>.

30"""Unwrapping module for SARvey."""

31import multiprocessing

32from os.path import join, dirname

33import time

34from typing import Union

36import matplotlib.pyplot as plt

37import numpy as np

38from kamui import unwrap_arbitrary

39from scipy.sparse import csr_matrix

40from scipy.sparse.csgraph import structural_rank

41from scipy.sparse.linalg import lsqr

42from scipy.optimize import minimize

43from logging import Logger

45from mintpy.utils import ptime

47import sarvey.utils as ut

48from sarvey.ifg_network import IfgNetwork

49from sarvey.objects import Network, NetworkParameter, AmplitudeImage

52def objFuncTemporalCoherence(x, *args):

53 """Compute temporal coherence from parameters and phase. To be used as objective function for optimization.

55 Parameters

56 ----------

57 x: np.ndarray

58 Search space for the DEM error in a 1D grid.

59 args: tuple

60 Additional arguments: (design_mat, obs_phase, scale_vel, scale_demerr).

62 Returns

63 -------

64 1 - gamma: float

65 """

66 (design_mat, obs_phase, scale_vel, scale_demerr) = args

68 # equalize the gradients in both directions

69 x[0] *= scale_demerr

70 x[1] *= scale_vel

72 pred_phase = np.matmul(design_mat, x)

73 res = (obs_phase - pred_phase.T).ravel()

74 gamma = np.abs(np.mean(np.exp(1j * res)))

75 return 1 - gamma

78def gridSearchTemporalCoherence(*, demerr_grid: np.ndarray, vel_grid: np.ndarray, design_mat: np.ndarray,

79 obs_phase: np.ndarray):

80 """Grid search which maximizes the temporal coherence as the objective function.

82 Parameters

83 ----------

84 demerr_grid: np.ndarray

85 Search space for the DEM error in a 2D grid.

86 vel_grid: np.ndarray

87 Search space for the velocity in a 2D grid.

88 design_mat: np.ndarray

89 Design matrix for estimating parameters from arc phase.

90 obs_phase: np.ndarray

91 Observed phase of the arc.

93 Returns

94 -------

95 demerr: float

96 estimated DEM error.

97 vel: float

98 estimated velocity.

99 gamma: float

100 estimated temporal coherence.

101 """

102 demerr_grid_flat = demerr_grid.flatten()

103 vel_grid_flat = vel_grid.flatten()

104 gamma_flat = np.array(

105 [1 - objFuncTemporalCoherence(np.array([demerr_grid_flat[i], vel_grid_flat[i]]),

106 design_mat, obs_phase, 1, 1)

107 for i in range(demerr_grid_flat.shape[0])])

108 gamma = gamma_flat.reshape(demerr_grid.shape)

109 idx_max_gamma = np.argmax(gamma_flat)

110

111 # return demerr_grid_flat[idx_max_gamma], vel_grid_flat[idx_max_gamma], gamma_flat[idx_max_gamma]

112 return demerr_grid_flat[idx_max_gamma], vel_grid_flat[idx_max_gamma], gamma

113

114

115def findOptimum(*, obs_phase: np.ndarray, design_mat: np.ndarray, val_range: np.ndarray):

116 """Find optimal value within a one dimensional search space that fits to the observed phase.

117

118 Parameters

119 ----------

120 obs_phase: np.ndarray

121 Observed phase of the arc.

122 design_mat: np.ndarray

123 Design matrix for estimating parameters from arc phase.

124 val_range: np.ndarray

125 Range of possible values for the solution. Can be either for DEM error or velocity.

126

127 Returns

128 -------

129 opt_val: scipy.optimize.minimize return value

130 gamma: float

131 pred_phase: np.ndarray

132 """

133 pred_phase = design_mat[:, np.newaxis] * val_range[np.newaxis, :] # broadcasting

134 if len(obs_phase.shape) == 2:

135 # step densification

136 res = obs_phase[:, np.newaxis, :] - pred_phase.T

137 res = np.moveaxis(res, 0, 1)

138 res = res.reshape((pred_phase.shape[1], -1)) # combine residuals from all arcs

139 else:

140 # step consistency check

141 res = obs_phase - pred_phase.T

142

143 gamma = np.abs(np.mean(np.exp(1j * res), axis=1))

144 max_idx = np.argmax(gamma)

145 opt_val = val_range[max_idx]

146 return opt_val, gamma[max_idx], pred_phase[:, max_idx]

147

148

149def oneDimSearchTemporalCoherence(*, demerr_range: np.ndarray, vel_range: np.ndarray, obs_phase: np.ndarray,

150 design_mat: np.ndarray):

151 """One dimensional search for maximum temporal coherence that fits the observed arc phase.

152

153 Parameters

154 ----------

155 demerr_range: np.ndarray

156 Search space for the DEM error in a 1D grid.

157 vel_range: np.ndarray

158 Search space for the velocity in a 1D grid.

159 design_mat: np.ndarray

160 Design matrix for estimating parameters from arc phase.

161 obs_phase: np.ndarray

162 Observed phase of the arc.

163

164 Returns

165 -------

166 demerr: float

167 vel: float

168 gamma: float

169 """

170 demerr, gamma_demerr, pred_phase_demerr = findOptimum(

171 obs_phase=obs_phase,

172 design_mat=design_mat[:, 0],

173 val_range=demerr_range

174 )

175

176 vel, gamma_vel, pred_phase_vel = findOptimum(

177 obs_phase=obs_phase,

178 design_mat=design_mat[:, 1],

179 val_range=vel_range

180 )

181

182 if gamma_vel > gamma_demerr:

183 demerr, gamma_demerr, pred_phase_demerr = findOptimum(

184 obs_phase=obs_phase - pred_phase_vel,

185 design_mat=design_mat[:, 0],

186 val_range=demerr_range

187 )

188 vel, gamma_vel, pred_phase_vel = findOptimum(

189 obs_phase=obs_phase - pred_phase_demerr,

190 design_mat=design_mat[:, 1],

191 val_range=vel_range

192 )

193 else:

194 vel, gamma_vel, pred_phase_vel = findOptimum(

195 obs_phase=obs_phase - pred_phase_demerr,

196 design_mat=design_mat[:, 1],

197 val_range=vel_range

198 )

199 demerr, gamma_demerr, pred_phase_demerr = findOptimum(

200 obs_phase=obs_phase - pred_phase_vel,

201 design_mat=design_mat[:, 0],

202 val_range=demerr_range

203 )

204

205 # improve initial estimate with gradient descent approach

206 scale_demerr = demerr_range.max()

207 scale_vel = vel_range.max()

208

209 demerr, vel, gamma = gradientSearchTemporalCoherence(

210 scale_vel=scale_vel,

211 scale_demerr=scale_demerr,

212 obs_phase=obs_phase,

213 design_mat=design_mat,

214 x0=np.array([demerr / scale_demerr,

215 vel / scale_vel]).T

216 )

217

218 pred_phase = np.matmul(design_mat, np.array([demerr, vel]))

219 res = (obs_phase - pred_phase.T).ravel()

220 gamma = np.abs(np.mean(np.exp(1j * res)))

221 return demerr, vel, gamma

222

223

224def gradientSearchTemporalCoherence(*, scale_vel: float, scale_demerr: float, obs_phase: np.ndarray,

225 design_mat: np.ndarray, x0: np.ndarray):

226 """GradientSearchTemporalCoherence.

227

228 Parameters

229 ----------

230 scale_demerr: float

231 Scaling factor for DEM error to equalize the axis of the search space.

232 scale_vel: float

233 Scaling factor for velocity to equalize the axis of the search space.

234 design_mat: np.ndarray

235 Design matrix for estimating parameters from arc phase.

236 obs_phase: np.ndarray

237 Observed phase of the arc.

238 x0: np.ndarray

239 Initial values for optimization.

240

241 Returns

242 -------

243 demerr: float

244 vel: float

245 gamma: float

246 """

247 opt_res = minimize(

248 objFuncTemporalCoherence,

249 x0,

250 args=(design_mat, obs_phase, scale_vel, scale_demerr),

251 bounds=((-1, 1), (-1, 1)),

252 method='L-BFGS-B'

253 )

254 gamma = 1 - opt_res.fun

255 demerr = opt_res.x[0] * scale_demerr

256 vel = opt_res.x[1] * scale_vel

257 return demerr, vel, gamma

258

259

260def launchAmbiguityFunctionSearch(parameters: tuple):

261 """Wrap for launching ambiguity function for temporal unwrapping in parallel.

262

263 Parameters

264 ----------

265 parameters: tuple

266 Arguments for temporal unwrapping in parallel.

267

268 Returns

269 -------

270 arc_idx_range: np.ndarray

271 demerr: np.ndarray

272 vel: np.ndarray

273 gamma: np.ndarray

274 """

275 (arc_idx_range, num_arcs, phase, slant_range, loc_inc, ifg_net_obj, wavelength, velocity_bound, demerr_bound,

276 num_samples) = parameters

277

278 demerr = np.zeros((num_arcs, 1), dtype=np.float32)

279 vel = np.zeros((num_arcs, 1), dtype=np.float32)

280 gamma = np.zeros((num_arcs, 1), dtype=np.float32)

281

282 design_mat = np.zeros((ifg_net_obj.num_ifgs, 2), dtype=np.float32)

283

284 demerr_range = np.linspace(-demerr_bound, demerr_bound, num_samples)

285 vel_range = np.linspace(-velocity_bound, velocity_bound, num_samples)

286

287 # prog_bar = ptime.progressBar(maxValue=num_arcs)

288

289 factor = 4 * np.pi / wavelength

290

291 for k in range(num_arcs):

292 design_mat[:, 0] = factor * ifg_net_obj.pbase_ifg / (slant_range[k] * np.sin(loc_inc[k]))

293 design_mat[:, 1] = factor * ifg_net_obj.tbase_ifg

294

295 demerr[k], vel[k], gamma[k] = oneDimSearchTemporalCoherence(

296 demerr_range=demerr_range,

297 vel_range=vel_range,

298 obs_phase=phase[k, :],

299 design_mat=design_mat

300 )

301

302 return arc_idx_range, demerr, vel, gamma

303

304

305def temporalUnwrapping(*, ifg_net_obj: IfgNetwork, net_obj: Network, wavelength: float, velocity_bound: float,

306 demerr_bound: float, num_samples: int, num_cores: int = 1, logger: Logger) -> \

307 tuple[np.ndarray, np.ndarray, np.ndarray]:

308 """Solve ambiguities for every arc in spatial Network object.

309

310 Parameters

311 ----------

312 ifg_net_obj: IfgNetwork

313 The IfgNetwork object.

314 net_obj: Network

315 The Network object.

316 wavelength: float

317 The wavelength.

318 velocity_bound: float

319 The velocity bound.

320 demerr_bound: float

321 The DEM error bound.

322 num_samples: int

323 The number of samples for the search space.

324 num_cores: int

325 Number of cores to be used. Default is 1.

326 logger: Logger

327 Logging handler.

328

329 Returns

330 -------

331 demerr: np.ndarray

332 vel: np.ndarray

333 gamma: np.ndarray

334 """

335 msg = "#" * 10

336 msg += " TEMPORAL UNWRAPPING: AMBIGUITY FUNCTION "

337 msg += "#" * 10

338 logger.info(msg=msg)

339

340 start_time = time.time()

341

342 if num_cores == 1:

343 args = (

344 np.arange(net_obj.num_arcs), net_obj.num_arcs, net_obj.phase,

345 net_obj.slant_range, net_obj.loc_inc, ifg_net_obj, wavelength, velocity_bound, demerr_bound, num_samples)

346 arc_idx_range, demerr, vel, gamma = launchAmbiguityFunctionSearch(parameters=args)

347 else:

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

349 pool = multiprocessing.Pool(processes=num_cores)

350

351 demerr = np.zeros((net_obj.num_arcs, 1), dtype=np.float32)

352 vel = np.zeros((net_obj.num_arcs, 1), dtype=np.float32)

353 gamma = np.zeros((net_obj.num_arcs, 1), dtype=np.float32)

354

355 num_cores = net_obj.num_arcs if num_cores > net_obj.num_arcs else num_cores # avoids having more samples then

356 # cores

357 idx = ut.splitDatasetForParallelProcessing(num_samples=net_obj.num_arcs, num_cores=num_cores)

358

359 args = [(

360 idx_range,

361 idx_range.shape[0],

362 net_obj.phase[idx_range, :],

363 net_obj.slant_range[idx_range],

364 net_obj.loc_inc[idx_range],

365 ifg_net_obj,

366 wavelength,

367 velocity_bound,

368 demerr_bound,

369 num_samples) for idx_range in idx]

370

371 results = pool.map(func=launchAmbiguityFunctionSearch, iterable=args)

372

373 # retrieve results

374 for i, demerr_i, vel_i, gamma_i in results:

375 demerr[i] = demerr_i

376 vel[i] = vel_i

377 gamma[i] = gamma_i

378

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

380 logger.info(msg="Finished temporal unwrapping.")

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

382 return demerr, vel, gamma

383

384

385def launchSpatialUnwrapping(parameters: tuple) -> tuple[np.ndarray, np.ndarray]:

386 """LaunchSpatialUnwrapping.

387

388 Parameters

389 ----------

390 parameters: tuple

391 idx_range, num_ifgs, num_points, edges, phase

392

393 Returns

394 -------

395 idx_range: np.ndarray

396 unw_phase: np.ndarray

397 """

398 # Unpack the parameters

399 (idx_range, num_ifgs, num_points, method, edges, phase) = parameters

400

401 prog_bar = ptime.progressBar(maxValue=num_ifgs)

402

403 unw_phase = np.zeros((num_points, num_ifgs), dtype=np.float32)

404

405 # Perform the PUMA phase unwrapping

406 for i in range(num_ifgs):

407 if method == "puma":

408 unw_phase[:, i] = unwrap_arbitrary(

409 psi=phase[:, i],

410 edges=edges,

411 simplices=None,

412 method="gc",

413 period=2*np.pi,

414 start_i=0,

415 p=0.2

416 )

417 else:

418 unw_phase[:, i] = unwrap_arbitrary(

419 psi=phase[:, i],

420 edges=edges,

421 simplices=None, # todo: compute simplices for ILP

422 method="ilp",

423 period=2*np.pi,

424 start_i=0,

425 )

426 prog_bar.update(value=i + 1, every=1,

427 suffix='{}/{} ifgs unwrapped. '.format(i + 1, num_ifgs))

428

429 unw_phase = unw_phase - np.mean(unw_phase, axis=0)

430 return idx_range, unw_phase

431

432

433def spatialUnwrapping(*, num_ifgs: int, num_points: int, phase: np.ndarray, edges: np.ndarray, method: str,

434 num_cores: int, logger: Logger):

435 """Spatial unwrapping of interferograms for a set of points.

436

437 Parameters

438 ----------

439 num_ifgs: int

440 Number of interferograms.

441 num_points: int

442 Number of points.

443 phase: np.ndarray

444 Phase of the interferograms at the points.

445 edges: np.ndarray

446 Edges/arcs of the graph.

447 method: str

448 Method for spatial unwrapping (puma or ilp).

449 num_cores: int

450 Number of cores to be used in multiprocessing.

451 logger: Logger

452 Logging handler.

453

454 Returns

455 -------

456 unw_phase: np.ndarray

457 Unwrapped phase of the interferograms at the points.

458 """

459 msg = "#" * 10

460 msg += f" SPATIAL UNWRAPPING: {method} "

461 msg += "#" * 10

462 logger.info(msg=msg)

463

464 start_time = time.time()

465

466 if num_cores == 1:

467 parameters = (

468 np.arange(num_ifgs),

469 num_ifgs,

470 num_points,

471 method,

472 edges,

473 phase

474 )

475 idx_range, unw_phase = launchSpatialUnwrapping(parameters=parameters)

476 else:

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

478 pool = multiprocessing.Pool(processes=num_cores)

479

480 unw_phase = np.zeros((num_points, num_ifgs), dtype=np.float32)

481 num_cores = num_ifgs if num_cores > num_ifgs else num_cores

482 # avoids having more samples than cores

483 idx = ut.splitDatasetForParallelProcessing(num_samples=num_ifgs, num_cores=num_cores)

484

485 args = [(

486 idx_range,

487 idx_range.shape[0],

488 num_points,

489 method,

490 edges,

491 phase[:, idx_range]) for idx_range in idx]

492 results = pool.map(func=launchSpatialUnwrapping, iterable=args)

493

494 # retrieve results

495 for i, phase in results:

496 unw_phase[:, i] = phase

497

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

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

500

501 return unw_phase

502

503

504def spatialParameterIntegrationIterative(*,

505 val_arcs: np.ndarray,

506 all_arcs: np.ndarray,

507 coord_xy: np.ndarray,

508 all_weights: np.ndarray,

509 spatial_ref_idx: int = 0,

510 res_tol: float = 1e-3,

511 max_rm_fraction: float = 0.001,

512 logger: Logger):

513 """Unwrapping double-difference arc parameters spatially.

514

515 The parameters at the arcs are integrated spatially to the points. The integration is done iteratively using

516 least-squares by removing the arcs with the highest residuals in each iteration.

517 The integration stops when the sum of the residuals is below a threshold.

518 Function is adopted from StaMPS software (Hooper et al., 2007).

519

520 Parameters

521 ----------

522 val_arcs: np.ndarray

523 Value at the arcs (e.g. DEM error, velocity).

524 all_arcs: np.ndarray

525 Arcs of the spatial network.

526 coord_xy: np.ndarray

527 Radar coordinates of the points in the spatial network.

528 all_weights: np.ndarray

529 Weights of the arcs (e.g. temporal coherence from temporal unwrapping)

530 spatial_ref_idx: int

531 Index of the spatial reference point (default = 0). Can be arbitrary.

532 res_tol: float

533 Threshold on the sum of the residual phase (default = 1e-3). Convergence criterion.

534 max_rm_fraction: float

535 Fraction of the arcs that are removed in each iteration (default = 0.001).

536 logger: Logger

537 Logging handler

538

539 Returns

540 -------

541 val_points: np.ndarray

542 Estimated parameters at the points resulting from the integration of the parameters at the arcs.

543 """

544 all_arcs = np.array(all_arcs)

545 num_points = coord_xy.shape[0]

546 num_arcs = all_arcs.shape[0]

547

548 # create design matrix

549 a = np.zeros((num_arcs, num_points))

550 for i in range(num_arcs):

551 a[i, all_arcs[i][0]] = 1

552 a[i, all_arcs[i][1]] = -1

553

554 # find the number of arcs per point

555 arcs_per_point = np.zeros(num_points, )

556

557 for i in range(num_points):

558 arcs_per_point[i] = np.where(a[:, i] != 0)[0].shape[0]

559

560 # remove reference point from design matrix

561 all_a = csr_matrix(all_weights * np.delete(a, spatial_ref_idx, 1))

562

563 # don't even start if the network is not connected

564 if structural_rank(all_a) < all_a.shape[1]:

565 logger.exception(msg="Spatial point network is not connected. Phase cannot be unwrapped!")

566 raise Exception

567

568 # set n_bad to maximum fraction of bad edges that can be removed

569 n_bad = np.ceil(num_arcs * max_rm_fraction).astype(np.int16)

570

571 # initialize output

572 val_points = np.zeros((num_points,))

573 points_idx = np.ones((num_points,), dtype=bool)

574 points_idx[spatial_ref_idx] = False

575 x_hat = np.zeros((num_points - 1,))

576

577 start_time = time.time()

578

579 arcs = all_arcs

580 obv_vec = val_arcs.reshape(-1, ) * all_weights.reshape(-1, )

581 a = all_a

582 weights = all_weights

583 num_arcs = obv_vec.size

584

585 r = None

586 num_arcs_save = None

587 arcs_save = None

588 a_save = None

589 weights_save = None

590 obv_vec_save = None

591 i = 0

592 while True:

593 if structural_rank(a) >= a.shape[1]:

594 x_hat[:] = lsqr(a, obv_vec)[0]

595

596 # store the current version of variables, being able to go back to previous iteration if too many arcs

597 # removed

598 a_save = a

599 obv_vec_save = obv_vec

600 weights_save = weights

601 arcs_save = arcs

602 num_arcs_save = num_arcs

603

604 # compute residuals

605 r = obv_vec - np.matmul(a.toarray(), x_hat)

606

607 else: # network is not connected anymore, remove less psPoints and try again

608 # x_hat = np.linalg.lstsq(a_save, obv_vec_save, rcond=None)[0] # unclear: I think it is not necessary to

609 # recompute the inversion.

610 n_bad = np.ceil(n_bad / 10).astype(np.int16) # remove less point

611

612 if np.all(np.abs(r) < res_tol):

613 break

614 else:

615 # drop arcs with the highest residuals, but only drop max one arc per point

616 ps_w_dropped_arc = np.zeros((num_points,))

617 good_arc_idx = np.ones((num_arcs_save,), dtype=bool)

618 r_sort_idx = np.abs(r).argsort()[::-1] # descending order, makes for loop easier

619

620 for j in range(n_bad): # remove arcs one by one

621 bad_arc_idx = r_sort_idx[j]

622 ps_idx0 = arcs_save[bad_arc_idx][0]

623 ps_idx1 = arcs_save[bad_arc_idx][1]

624 if (ps_w_dropped_arc[ps_idx0] == 0) and (ps_w_dropped_arc[

625 ps_idx1] == 0): # if arc not already dropped for either

626 # point of current arc drop current arc

627 good_arc_idx[bad_arc_idx] = False

628 # mark both psPoints from the arc as having an arc dropped

629 ps_w_dropped_arc[ps_idx0] = 1

630 ps_w_dropped_arc[ps_idx1] = 1

631

632 # update all variables for next iteration

633 arcs = arcs_save[good_arc_idx, :]

634 obv_vec = obv_vec_save[good_arc_idx]

635 a = a_save[good_arc_idx, :]

636 weights = weights_save[good_arc_idx]

637 num_arcs = obv_vec.size

638

639 i += 1

640

641 val_points[points_idx] = x_hat

642

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

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

645

646 return val_points

647

648

649def spatialParameterIntegration(*,

650 val_arcs: np.ndarray,

651 arcs: np.ndarray,

652 coord_xy: np.ndarray,

653 weights: np.ndarray,

654 spatial_ref_idx: int = 0,

655 logger: Logger):

656 """Unwrapping double-difference arc parameters spatially.

657

658 The parameters at the arcs are integrated spatially to the points. The integration is done using least-squares.

659

660 Parameters

661 ----------

662 val_arcs: np.ndarray

663 Value at the arcs (e.g. DEM error, velocity).

664 arcs: np.ndarray

665 Arcs of the spatial network.

666 coord_xy: np.ndarray

667 Radar coordinates of the points in the spatial network.

668 weights: np.ndarray

669 Weights of the arcs (e.g. temporal coherence from temporal unwrapping)

670 spatial_ref_idx: int

671 Index of the spatial reference point (default = 0). Can be arbitrary.

672 logger: Logger

673 Logging handler

674

675 Returns

676 -------

677 val_points: np.ndarray

678 Estimated parameters at the points resulting from the integration of the parameters at the arcs.

679 """

680 arcs = np.array(arcs)

681 num_points = coord_xy.shape[0]

682 num_arcs = arcs.shape[0]

683

684 # create design matrix

685 design_mat = np.zeros((num_arcs, num_points))

686 for i in range(num_arcs):

687 design_mat[i, arcs[i][0]] = 1

688 design_mat[i, arcs[i][1]] = -1

689

690 # remove reference point from design matrix

691 design_mat = csr_matrix(weights * np.delete(design_mat, spatial_ref_idx, 1))

692

693 # don't even start if the network is not connected

694 if structural_rank(design_mat) < design_mat.shape[1]:

695 raise Exception("Spatial point network is not connected. Cannot integrate parameters spatially!")

696

697 start_time = time.time()

698

699 obv_vec = val_arcs.reshape(-1, ) * weights.reshape(-1, )

700

701 x_hat = lsqr(design_mat, obv_vec)[0]

702

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

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

705

706 val_points = np.zeros((num_points,))

707 points_idx = np.ones((num_points,), dtype=bool)

708 points_idx[spatial_ref_idx] = False

709 val_points[points_idx] = x_hat

710

711 return val_points

712

713

714def computeNumArcsPerPoints(*, net_obj: Network, point_id: np.ndarray,

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

716 """Remove Points with less than specified number of arcs.

717

718 Parameters

719 ----------

720 net_obj: Network

721 The spatial Network object.

722 point_id: np.ndarray

723 ID of the points in the network.

724 logger: Logger

725 Logging handler.

726

727 Returns

728 -------

729 design_mat: np.ndarray

730 Design matrix of the spatial network

731 arcs_per_point: np.ndarray

732 Number of arcs that each point is connected with.

733 """

734 logger.info(msg="Removal of arcs and PSC that cannot be tested.")

735

736 num_points = point_id.shape[0]

737

738 # create design matrix

739 design_mat = np.zeros((net_obj.num_arcs, num_points))

740 for i in range(net_obj.num_arcs):

741 design_mat[i, net_obj.arcs[i][0]] = 1

742 design_mat[i, net_obj.arcs[i][1]] = -1

743

744 # find the number of arcs per point

745 arcs_per_point = np.zeros(num_points, )

746

747 for i in range(num_points):

748 arcs_per_point[i] = np.where(design_mat[:, i] != 0)[0].shape[0]

749

750 return design_mat, arcs_per_point

751

752

753def computeAvgCoherencePerPoint(*, net_obj: Network, point_id: np.ndarray, logger: Logger) -> np.ndarray:

754 """Compute the average coherence from all arcs that a point is connected with. Used to remove incoherent points.

755

756 Parameters

757 ----------

758 net_obj: Network

759 The Network object.

760 point_id: np.ndarray

761 ID of the points in the network.

762 logger: Logger

763 Logging handler.

764

765 Returns

766 -------

767 mean_gamma_point: np.ndarray

768 Average coherence per point

769 """

770 logger.info(msg="Removal of points whose arcs are incoherent in average.")

771

772 num_points = point_id.shape[0]

773

774 # create design matrix

775 a = np.zeros((net_obj.num_arcs, num_points))

776 for i in range(net_obj.num_arcs):

777 a[i, net_obj.arcs[i][0]] = net_obj.gamma[i]

778 a[i, net_obj.arcs[i][1]] = net_obj.gamma[i]

779

780 a[a == 0] = np.nan

781 mean_gamma_point = np.nanmean(a, axis=0)

782

783 return mean_gamma_point

784

785

786def removeArcsByPointMask(*, net_obj: Union[Network, NetworkParameter], point_id: np.ndarray, coord_xy: np.ndarray,

787 p_mask: np.ndarray, design_mat: np.ndarray,

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

789 """Remove all entries related to the arc observations connected to the points which have a False value in p_mask.

790

791 Parameters

792 ----------

793 net_obj: Network

794 The Network object.

795 point_id: np.ndarray

796 ID of the points in the network.

797 coord_xy: np.ndarray

798 Radar coordinates of the points in the spatial network.

799 p_mask: np.ndarray

800 Boolean mask with True for points to keep, and False for points to remove.

801 design_mat: np.ndarray

802 Design matrix describing the relation between arcs and points.

803 logger: Logger

804 Logging handler.

805

806 Returns

807 -------

808 net_obj: Network

809 Network object without the removed arcs and points.

810 point_id: np.ndarray

811 ID of the points in the network without the removed points.

812 coord_xy: np.ndarray

813 Radar coordinates of the points in the spatial network without the removed points.

814 design_mat: np.ndarray

815 Design matrix describing the relation between arcs and points without the removed points and arcs.

816 """

817 # find respective arcs

818 a_idx = list()

819 for p_idx in np.where(~p_mask)[0]:

820 a_idx.append(np.where(design_mat[:, p_idx] != 0)[0])

821

822 if len(a_idx) != 0:

823 a_idx = np.hstack(a_idx)

824 a_mask = np.ones((net_obj.num_arcs,), dtype=np.bool_)

825 a_mask[a_idx] = False

826 net_obj.removeArcs(mask=a_mask)

827 design_mat = design_mat[a_mask, :]

828 else:

829 a_idx = np.array(a_idx) # so I can check the size

830

831 # remove psPoints

832 point_id = point_id[p_mask]

833 design_mat = design_mat[:, p_mask]

834 coord_xy = coord_xy[p_mask, :]

835

836 # beside removing the arcs in "arcs", the tuple indices have to be changed to make them fit to new point indices

837 for p_idx in np.sort(np.where(~p_mask)[0])[::-1]:

838 net_obj.arcs[np.where((net_obj.arcs[:, 0] > p_idx)), 0] -= 1

839 net_obj.arcs[np.where((net_obj.arcs[:, 1] > p_idx)), 1] -= 1

840

841 logger.info(msg="Removed {} arc(s) connected to the removed point(s)".format(a_idx.size))

842 return net_obj, point_id, coord_xy, design_mat

843

844

845def removeGrossOutliers(*, net_obj: Network, point_id: np.ndarray, coord_xy: np.ndarray, min_num_arc: int = 3,

846 quality_thrsh: float = 0.0,

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

848 """Remove both gross outliers which have many low quality arcs and points which are not well connected.

849

850 Parameters

851 ----------

852 net_obj: Network

853 The spatial Network object.

854 point_id: np.ndarray

855 ID of the points in the network.

856 coord_xy: np.ndarray

857 Radar coordinates of the points in the spatial network.

858 min_num_arc: int

859 Threshold on the minimal number of arcs per point. Default = 3.

860 quality_thrsh: float

861 Threshold on the temporal coherence of the arcs. Default = 0.0.

862 logger: Logger

863 Logging handler.

864

865 Returns

866 -------

867 net_obj: Network

868 Network object without the removed arcs and points.

869 point_id: np.ndarray

870 ID of the points in the network without the removed points.

871 coord_xy: np.ndarray

872 Radar coordinates of the points in the spatial network without the removed points.

873 a: np.ndarray

874 Design matrix describing the relation between arcs and points without the removed points and arcs.

875 """

876 logger.info(msg="Detect points with low quality arcs (mean): < {}".format(quality_thrsh))

877 mean_gamma_point = computeAvgCoherencePerPoint(net_obj=net_obj,

878 point_id=point_id, logger=logger)

879 # not yet removed, because arcs are removed separately

880 p_mask_mean_coh = (mean_gamma_point >= quality_thrsh).ravel()

881 logger.info(msg="Detected {} point(s) with mean coherence of all connected arcs < {} ".format(

882 p_mask_mean_coh[~p_mask_mean_coh].shape[0], quality_thrsh))

883

884 logger.info(msg="Removal of low quality arcs: < {}".format(quality_thrsh))

885 a_mask = (net_obj.gamma >= quality_thrsh).ravel()

886 logger.info(msg="Removed {} arc(s)".format(a_mask[~a_mask].shape[0]))

887 net_obj.removeArcs(mask=a_mask)

888

889 design_mat, arcs_per_point = computeNumArcsPerPoints(net_obj=net_obj, point_id=point_id, logger=logger)

890

891 p_mask_num_arcs = (arcs_per_point >= min_num_arc).ravel()

892 logger.info(msg="Detected {} point(s) with less than {} arcs".format(p_mask_num_arcs[~p_mask_num_arcs].shape[0],

893 min_num_arc))

894

895 # remove them jointly

896 p_mask = p_mask_num_arcs & p_mask_mean_coh

897 logger.info(msg="Remove {} point(s)".format(p_mask[~p_mask].shape[0]))

898 net_obj, point_id, coord_xy, design_mat = removeArcsByPointMask(net_obj=net_obj, point_id=point_id,

899 coord_xy=coord_xy, p_mask=p_mask,

900 design_mat=design_mat, logger=logger)

901 return net_obj, point_id, coord_xy, design_mat

902

903

904def parameterBasedNoisyPointRemoval(*, net_par_obj: NetworkParameter, point_id: np.ndarray, coord_xy: np.ndarray,

905 design_mat: np.ndarray, rmse_thrsh: float = 0.02, num_points_remove: int = 1,

906 bmap_obj: AmplitudeImage = None, bool_plot: bool = False,

907 logger: Logger):

908 """Remove Points during spatial integration step if residuals at many connected arcs are high.

909

910 The idea is similar to outlier removal in DePSI, but without hypothesis testing.

911 It can be used as a preprocessing step to spatial integration.

912 The points are removed based on the RMSE computed from the residuals of the parameters (DEM error, velocity) per

913 arc. The point with the highest RMSE is removed in each iteration. The process stops when the maximum RMSE is below

914 a threshold.

915

916

917 Parameters

918 ----------

919 net_par_obj: NetworkParameter

920 The spatial NetworkParameter object containing the parameters estimates at each arc.

921 point_id: np.ndarray

922 ID of the points in the network.

923 coord_xy: np.ndarray

924 Radar coordinates of the points in the spatial network.

925 design_mat: np.ndarray

926 Design matrix describing the relation between arcs and points.

927 rmse_thrsh: float

928 Threshold for the RMSE of the residuals per point. Default = 0.02.

929 num_points_remove: int

930 Number of points to remove in each iteration. Default = 1.

931 bmap_obj: AmplitudeImage

932 Basemap object for plotting. Default = None.

933 bool_plot: bool

934 Plot the RMSE per point. Default = False.

935 logger: Logger

936 Logging handler.

937

938 Returns

939 -------

940 spatial_ref_id: int

941 ID of the spatial reference point.

942 point_id: np.ndarray

943 ID of the points in the network without the removed points.

944 net_par_obj: NetworkParameter

945 The NetworkParameter object without the removed points.

946 """

947 msg = "#" * 10

948 msg += " NOISY POINT REMOVAL BASED ON ARC PARAMETERS "

949 msg += "#" * 10

950 logger.info(msg=msg)

951

952 num_points = point_id.shape[0]

953

954 logger.info(msg="Selection of the reference PSC")

955 # select one of the two pixels which are connected via the arc with the highest quality

956 spatial_ref_idx = np.where(design_mat[np.argmax(net_par_obj.gamma), :] != 0)[0][0]

957 coord_xy = np.delete(arr=coord_xy, obj=spatial_ref_idx, axis=0)

958 spatial_ref_id = point_id[spatial_ref_idx]

959 point_id = np.delete(arr=point_id, obj=spatial_ref_idx, axis=0)

960 num_points -= 1

961

962 # remove reference point from design matrix

963 design_mat = net_par_obj.gamma * np.delete(arr=design_mat, obj=spatial_ref_idx, axis=1)

964

965 logger.info(msg="Spatial integration to detect noisy point")

966 start_time = time.time()

967

968 it_count = 0

969 while True:

970 logger.info(msg="ITERATION: {}".format(it_count))

971 design_mat = csr_matrix(design_mat)

972

973 if structural_rank(design_mat) < design_mat.shape[1]:

974 logger.error(msg="Singular normal matrix. Network is no longer connected!")

975 # point_id = np.sort(np.hstack([spatial_ref_id, point_id]))

976 # return spatial_ref_id, point_id, net_par_obj

977 raise ValueError

978 # demerr

979 obv_vec = net_par_obj.demerr.reshape(-1, )

980 demerr_points = lsqr(design_mat.toarray(), obv_vec * net_par_obj.gamma.reshape(-1, ))[0]

981 r_demerr = obv_vec - np.matmul(design_mat.toarray(), demerr_points)

982

983 # vel

984 obv_vec = net_par_obj.vel.reshape(-1, )

985 vel_points = lsqr(design_mat.toarray(), obv_vec * net_par_obj.gamma.reshape(-1, ))[0]

986 r_vel = obv_vec - np.matmul(design_mat.toarray(), vel_points)

987

988 rmse_demerr = np.zeros((num_points,))

989 rmse_vel = np.zeros((num_points,))

990 for p in range(num_points):

991 r_mask = design_mat[:, p].toarray() != 0

992 rmse_demerr[p] = np.sqrt(np.mean(r_demerr[r_mask.ravel()].ravel() ** 2))

993 rmse_vel[p] = np.sqrt(np.mean(r_vel[r_mask.ravel()].ravel() ** 2))

994

995 rmse = rmse_vel.copy()

996 max_rmse = np.max(rmse.ravel())

997 logger.info(msg="Maximum RMSE DEM correction: {:.2f} m".format(np.max(rmse_demerr.ravel())))

998 logger.info(msg="Maximum RMSE velocity: {:.4f} m / year".format(np.max(rmse_vel.ravel())))

999

1000 if bool_plot:

1001 # vel

1002 ax = bmap_obj.plot(logger=logger)

1003 sc = ax.scatter(coord_xy[:, 1], coord_xy[:, 0], c=rmse_vel * 1000, s=3.5,

1004 cmap=plt.cm.get_cmap("autumn_r"), vmin=0, vmax=rmse_thrsh * 1000)

1005 plt.colorbar(sc, pad=0.03, shrink=0.5)

1006 ax.set_title("{}. iteration\nmean velocity - RMSE per point in [mm / year]".format(it_count))

1007 fig = ax.get_figure()

1008 plt.tight_layout()

1009 fig.savefig(join(dirname(net_par_obj.file_path), "pic", f"step_1_rmse_vel_{it_count}th_iter.png"),

1010 dpi=300)

1011 plt.close(fig)

1012

1013 # demerr

1014 ax = bmap_obj.plot(logger=logger)

1015 sc = ax.scatter(coord_xy[:, 1], coord_xy[:, 0], c=rmse_demerr, s=3.5,

1016 cmap=plt.cm.get_cmap("autumn_r"))

1017 plt.colorbar(sc, pad=0.03, shrink=0.5)

1018 ax.set_title("{}. iteration\nDEM correction - RMSE per point in [m]".format(it_count))

1019 fig = ax.get_figure()

1020 plt.tight_layout()

1021 fig.savefig(join(dirname(net_par_obj.file_path), "pic",

1022 f"step_1_rmse_dem_correction_{it_count}th_iter.png"),

1023 dpi=300)

1024 plt.close(fig)

1025

1026 if max_rmse <= rmse_thrsh:

1027 logger.info(msg="No noisy pixels detected.")

1028 break

1029

1030 # remove point with highest rmse

1031 p_mask = np.ones((num_points,), dtype=np.bool_)

1032 p_mask[np.argsort(rmse)[::-1][:num_points_remove]] = False # see description of function removeArcsByPointMask

1033 net_par_obj, point_id, coord_xy, design_mat = removeArcsByPointMask(net_obj=net_par_obj, point_id=point_id,

1034 coord_xy=coord_xy, p_mask=p_mask,

1035 design_mat=design_mat.toarray(),

1036 logger=logger)

1037 num_points -= num_points_remove

1038 it_count += 1

1039

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

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

1042

1043 # add spatialRefIdx back to point_id

1044 point_id = np.sort(np.hstack([spatial_ref_id, point_id]))

1045 return spatial_ref_id, point_id, net_par_obj