Coverage for sarvey/objects.py: 88%
363 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-17 12:36 +0000
1#!/usr/bin/env python
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/>.
30"""Objects module for SARvey."""
31import os
32from os.path import join, dirname, exists, basename
33from typing import Optional
34import h5py
35import matplotlib.pyplot as plt
36import numpy as np
37from pyproj import Proj, CRS
38from pyproj.aoi import AreaOfInterest
39from pyproj.database import query_utm_crs_info
40from logging import Logger
42from miaplpy.objects.slcStack import slcStack
43from mintpy.utils import readfile
44from mintpy.utils.plot import auto_flip_direction
46from sarvey.ifg_network import IfgNetwork
49class AmplitudeImage:
50 """AmplitudeImage."""
52 def __init__(self, *, file_path: str):
53 """Init.
55 Parameters
56 ----------
57 file_path: str
58 path to filename
59 """
60 self.width = None
61 self.length = None
62 self.file_path = file_path
63 self.background_map = None
64 self.orbit_direction = None
66 def prepare(self, *, slc_stack_obj: slcStack, img: np.ndarray, logger: Logger):
67 """Read the SLC stack, compute the mean amplitude image and store it into a file.
69 Parameters
70 ----------
71 slc_stack_obj: slcStack
72 object of class slcStack from MiaplPy
73 img: np.ndarray
74 amplitude image, e.g. the mean over time
75 logger: Logger
76 Logging handler
77 """
78 self.orbit_direction = slc_stack_obj.metadata["ORBIT_DIRECTION"]
79 self.length = slc_stack_obj.length
80 self.width = slc_stack_obj.width
82 self.background_map = img
84 logger.info(msg="write data to {}...".format(self.file_path))
86 if exists(self.file_path):
87 os.remove(self.file_path)
89 with h5py.File(self.file_path, 'w') as f:
90 f.create_dataset('background_map', data=self.background_map)
91 f.attrs["ORBIT_DIRECTION"] = self.orbit_direction
92 f.attrs["LENGTH"] = self.length
93 f.attrs["WIDTH"] = self.width
95 def open(self):
96 """Open."""
97 # print("read from {}".format(self.file_path))
99 with h5py.File(self.file_path, 'r') as f:
100 self.background_map = f["background_map"][:]
101 self.orbit_direction = f.attrs["ORBIT_DIRECTION"]
102 self.length = f.attrs["LENGTH"]
103 self.width = f.attrs["WIDTH"]
105 def plot(self, *, ax: plt.Axes = None, logger: Logger):
106 """Plot the mean amplitude image as a background map.
108 Parameters
109 ----------
110 ax: plt.Axes
111 axes for plotting (default: None, a new figure will be created).
112 logger: Logger
113 Logging handler.
115 Return
116 ------
117 ax: plt.Axes
118 axes object.
119 """
120 if self.background_map is None:
121 try:
122 self.open()
123 except OSError as e:
124 logger.error(msg="Could not open file: {}".format(e))
125 fig = plt.figure(figsize=(15, 5))
126 ax = fig.add_subplot()
127 logger.error(msg="Orbit direction not available.")
128 return ax
130 if ax is None:
131 fig = plt.figure(figsize=(15, 5))
132 ax = fig.add_subplot()
133 ax.imshow(self.background_map, cmap=plt.cm.get_cmap("gray"))
134 meta = {"ORBIT_DIRECTION": self.orbit_direction}
135 auto_flip_direction(meta, ax=ax, print_msg=False)
137 ax.set_xlabel("Range")
138 ax.set_ylabel("Azimuth")
140 return ax
143class CoordinatesUTM:
144 """Coordinates in UTM for all pixels in the radar image."""
146 def __init__(self, *, file_path: str, logger: Logger):
147 """Init.
149 Parameters
150 ----------
151 file_path: str
152 path to filename
153 logger: Logger
154 Logging handler.
155 """
156 self.file_path = file_path
157 self.coord_utm = None
158 self.logger = logger
160 def prepare(self, *, input_path: str):
161 """Read the slc stack, computes the mean amplitude image and stores it into a file.
163 Parameters
164 ----------
165 input_path: str
166 path to slcStack.h5 file.
167 """
168 log = self.logger
169 lat = readfile.read(input_path, datasetName='latitude')[0]
170 lon = readfile.read(input_path, datasetName='longitude')[0]
172 log.info(msg="Transform coordinates from latitude and longitude (WGS84) to North and East (UTM).")
173 # noinspection PyTypeChecker
174 utm_crs_list = query_utm_crs_info(
175 datum_name="WGS 84",
176 area_of_interest=AreaOfInterest(
177 west_lon_degree=np.nanmin(lon.ravel()),
178 south_lat_degree=np.nanmin(lat.ravel()),
179 east_lon_degree=np.nanmax(lon.ravel()),
180 north_lat_degree=np.nanmax(lat.ravel())),
181 contains=True)
182 utm_crs = CRS.from_epsg(utm_crs_list[0].code)
183 lola2utm = Proj(utm_crs)
184 self.coord_utm = np.array(lola2utm(lon, lat))
186 log.info(msg="write data to {}...".format(self.file_path))
188 if exists(self.file_path):
189 os.remove(self.file_path)
191 with h5py.File(self.file_path, 'w') as f:
192 f.create_dataset('coord_utm', data=self.coord_utm)
194 def open(self):
195 """Open."""
196 with h5py.File(self.file_path, 'r') as f:
197 self.coord_utm = f["coord_utm"][:]
200class BaseStack:
201 """Class for 3D image-like data stacks."""
203 def __init__(self, *, file: str = None, logger: Logger):
204 """Init.
206 Parameters
207 ----------
208 file: str
209 path to filename
210 logger: Logger
211 Logging handler.
212 """
213 self.file = file
214 self.logger = logger
215 self.metadata = None
216 self.num_time = None
217 self.length = None
218 self.width = None
219 self.f = None
221 def close(self, *, print_msg: bool = True):
222 """Close."""
223 try:
224 self.f.close()
225 if print_msg:
226 self.logger.info(msg='close file: {}'.format(basename(self.file)))
227 except Exception as e:
228 self.logger.exception(msg=e)
229 pass
230 return None
232 def getShape(self, *, dataset_name: str):
233 """Open file and read shape of dataset."""
234 with h5py.File(self.file, 'r') as f:
235 dshape = f[dataset_name].shape
236 return dshape
238 def read(self, *, dataset_name: str, box: Optional[tuple] = None, print_msg: bool = True):
239 """Read dataset from slc file.
241 Parameters
242 ----------
243 dataset_name: str
244 name of dataset
245 box: tuple
246 tuple of 4 int, indicating x0,y0,x1,y1 of range, or
247 tuple of 6 int, indicating x0,y0,z0,x1,y1,z1 of range
248 print_msg: bool
249 print message.
251 Returns
252 -------
253 data: np.ndarray
254 2D or 3D dataset
255 """
256 if print_msg:
257 self.logger.info(msg='reading box {} from file: {} ...'.format(box, self.file))
259 with h5py.File(self.file, 'r') as f:
260 self.metadata = dict(f.attrs)
262 ds = f[dataset_name]
263 if len(ds.shape) == 3:
264 self.length, self.width, self.num_time = ds.shape
265 else:
266 self.length, self.width = ds.shape
268 # Get Index in space/2_3 dimension
269 if box is None:
270 box = [0, 0, self.width, self.length]
272 if len(ds.shape) == 3:
273 if len(box) == 4:
274 data = ds[box[1]:box[3], box[0]:box[2], :]
275 if len(box) == 6:
276 data = ds[box[1]:box[4], box[0]:box[3], box[2]:box[5]]
277 else:
278 if len(box) == 6:
279 raise IndexError("Cannot read 3D box from 2D data.")
280 data = ds[box[1]:box[3], box[0]:box[2]]
282 for key, value in self.metadata.items():
283 try:
284 self.metadata[key] = value.decode('utf8')
285 except Exception:
286 self.metadata[key] = value
287 return data
289 def prepareDataset(self, dataset_name: str, dshape: tuple, dtype: object,
290 metadata: Optional[dict], mode: str = "w", chunks: [tuple, bool] = True):
291 """PrepareDataset. Creates a dataset in file with specified size without writing any data.
293 Parameters
294 ----------
295 dataset_name: str
296 name of dataset.
297 dshape: tuple
298 shape of dataset.
299 dtype: object
300 data type of dataset.
301 metadata: dict
302 metadata of dataset (e.g. WAVELENGTH, ORBIT_DIRECTION, etc.). Usually the same as in slcStack.h5.
303 mode: str
304 open mode ('w' for writing new file or 'a' for appending to existing file).
305 chunks: tuple
306 chunk size ('True'/'False' or tuple specifying the dimension of the chunks)
307 """
308 with h5py.File(self.file, mode) as f:
309 self.logger.info(msg="Prepare dataset: {d:<25} of {t:<25} in size of {s}".format(
310 d=dataset_name,
311 t=str(dtype),
312 s=dshape))
314 f.create_dataset(dataset_name,
315 shape=dshape,
316 dtype=dtype,
317 chunks=chunks)
319 # write attributes
320 metadata = dict(metadata)
321 for key in metadata.keys():
322 f.attrs[key] = metadata[key]
324 return
326 def writeToFileBlock(self, *, data: np.ndarray, dataset_name: str, block: Optional[tuple] = None, mode: str = 'a',
327 print_msg: bool = True):
328 """Write data to existing HDF5 dataset in disk block by block.
330 Parameters
331 ----------
332 data: np.ndarray
333 1/2/3D matrix.
334 dataset_name: str
335 dataset name.
336 block: list
337 the list can contain 2, 4 or 6 integers indicating: [zStart, zEnd, yStart, yEnd, xStart, xEnd].
338 mode: str
339 open mode ('w' for writing new file or 'a' for appending to existing file).
340 print_msg: bool
341 print message.
343 Returns
344 --------
345 file: str
346 path to file
347 """
348 if block is None:
349 # data shape
350 if isinstance(data, list):
351 shape = (len(data),)
352 else:
353 shape = data.shape
355 if len(shape) == 1:
356 block = [0, shape[0]]
357 elif len(shape) == 2:
358 block = [0, shape[0],
359 0, shape[1]]
360 elif len(shape) == 3:
361 block = [0, shape[0],
362 0, shape[1],
363 0, shape[2]]
365 with h5py.File(self.file, mode) as f:
367 if print_msg:
368 self.logger.info(msg="writing dataset /{:<25} block: {}".format(dataset_name, block))
369 if len(block) == 6:
370 f[dataset_name][block[0]:block[1],
371 block[2]:block[3],
372 block[4]:block[5]] = data
374 elif len(block) == 4:
375 f[dataset_name][block[0]:block[1],
376 block[2]:block[3]] = data
378 elif len(block) == 2:
379 f[dataset_name][block[0]:block[1]] = data
381 return self.file
383 def writeToFile(self, *, data: np.ndarray, dataset_name: str, metadata: Optional[dict] = None, mode: str = 'a',
384 chunks: [tuple, bool] = True):
385 """Write the whole dataset to the file (not block-by-block).
387 Parameters
388 ----------
389 data: np.ndarray
390 3D data array.
391 dataset_name: str
392 name of dataset.
393 metadata: dict
394 metadata of dataset (e.g. WAVELENGTH, ORBIT_DIRECTION, etc.). Usually the same as in slcStack.h5.
395 mode: str
396 mode for opening the h5 file (e.g. write: 'w' or append: 'a')
397 chunks: tuple
398 chunk size ('True'/'False' or tuple specifying the dimension of the chunks)
399 """
400 # 3D dataset
401 self.logger.info(msg='create HDF5 file: {} with w mode'.format(self.file))
402 self.f = h5py.File(self.file, mode)
403 if dataset_name not in self.f:
404 self.logger.info(msg='create dataset /{n} of {t:<10} in size of {s}.'.format(n=dataset_name,
405 t=str(data.dtype),
406 s=data.shape))
407 self.f.create_dataset(dataset_name, data=data, chunks=chunks)
408 else:
409 self.logger.info(msg='overwrite dataset /{n} of {t:<10} in size of {s}.'.format(n=dataset_name,
410 t=str(data.dtype),
411 s=data.shape))
412 self.f[dataset_name] = data
414 # Attributes
415 if metadata is not None:
416 metadata = dict(metadata)
417 for key, value in metadata.items():
418 self.f.attrs[key] = str(value)
420 self.f.close()
421 self.logger.info(msg='finished writing to {}'.format(self.file))
422 return
425class Points:
426 """Points class for storing information about the selected scatterers."""
428 file_path: str
429 point_id: np.array
430 coord_xy: np.array
431 num_points: int
432 phase: np.array
433 wavelength: float
434 length: int
435 width: int
436 times: None
438 # etc.
440 def __init__(self, *, file_path: str, logger: Logger):
441 """Init.
443 Parameters
444 ----------
445 file_path: str
446 ath to filename
447 logger: Logger
448 Logging handler.
449 """
450 self.ifg_net_obj = IfgNetwork() # use parent class here which doesn't know and care about 'star' or 'sb'
451 self.coord_utm = None
452 self.coord_lalo = None
453 self.height = None
454 self.slant_range = None
455 self.loc_inc = None
456 self.file_path = file_path
457 self.logger = logger
459 def prepare(self, *, point_id: np.ndarray, coord_xy: np.ndarray, input_path: str):
460 """Assign point_id and radar coordinates to the object.
462 Store the point_id and radar coordinates of the scatterers in the object (not file) and read further
463 attributes from external files (ifg_network.h5, slcStack.h5, geometryRadar.h5, coordinates_utm.h5).
465 Parameters
466 ----------
467 point_id: np.ndarray
468 point_id of the scatterers.
469 coord_xy: np.ndarray
470 radar coordinates of the scatterers.
471 input_path: str
472 path to input files (slcStack.h5, geometryRadar.h5).
473 """
474 self.point_id = point_id
475 self.coord_xy = coord_xy
476 self.num_points = self.coord_xy.shape[0]
477 self.phase = None
478 self.openExternalData(input_path=input_path)
480 def writeToFile(self):
481 """Write data to .h5 file (num_points, coord_xy, point_id, phase)."""
482 self.logger.info(msg="write data to {}...".format(self.file_path))
484 if exists(self.file_path):
485 os.remove(self.file_path)
487 with h5py.File(self.file_path, 'w') as f:
488 f.attrs["num_points"] = self.num_points
489 f.create_dataset('coord_xy', data=self.coord_xy)
490 f.create_dataset('point_id', data=self.point_id)
491 f.create_dataset('phase', data=self.phase)
493 def open(self, input_path: str, other_file_path: str = None):
494 """Read data from file.
496 Read stored information from already existing .h5 file. This can be the file of the object itself. If the
497 data should be read from another file, the path to this file can be given as 'other_file_path'. Thereby, a new
498 Points object can be created with the data of another Points object.
500 Parameters
501 ----------
502 input_path: str
503 path to input files (slcStack.h5, geometryRadar.h5).
504 other_file_path: str
505 path to other .h5 file (default: None).
506 """
507 # 1) read own data: coord_xy, phase, point_id, num_points, reference_point_idx
508 if other_file_path is not None:
509 path = other_file_path
510 else:
511 path = self.file_path
512 self.logger.info(msg="read from {}".format(path))
514 with h5py.File(path, 'r') as f:
515 self.num_points = f.attrs["num_points"]
516 self.coord_xy = f["coord_xy"][:]
517 self.point_id = f["point_id"][:]
518 self.phase = f["phase"][:]
520 self.openExternalData(input_path=input_path)
522 def openExternalData(self, *, input_path: str):
523 """Load data which is stored in slcStack.h5, geometryRadar.h5, ifg_network.h5 and coordinates_utm.h5."""
524 # 1) read IfgNetwork
525 self.ifg_net_obj.open(path=join(dirname(self.file_path), "ifg_network.h5"))
527 # 2) read metadata from slcStack
528 slc_stack_obj = slcStack(join(input_path, "slcStack.h5"))
529 slc_stack_obj.open(print_msg=False)
530 self.wavelength = np.float64(slc_stack_obj.metadata["WAVELENGTH"])
531 self.length = slc_stack_obj.length # y-coordinate axis (azimut)
532 self.width = slc_stack_obj.width # x-coordinate axis (range)
534 # 3) read from geometry file
535 mask = self.createMask()
537 geom_path = join(input_path, "geometryRadar.h5")
539 # load geometry data
540 loc_inc, meta = readfile.read(geom_path, datasetName='incidenceAngle')
541 loc_inc *= np.pi / 180 # in [rad]
542 slant_range = readfile.read(geom_path, datasetName='slantRangeDistance')[0]
543 height = readfile.read(geom_path, datasetName='height')[0]
544 lat = readfile.read(geom_path, datasetName='latitude')[0]
545 lon = readfile.read(geom_path, datasetName='longitude')[0]
547 self.loc_inc = loc_inc[mask].ravel()
548 self.slant_range = slant_range[mask].ravel()
549 self.height = height[mask].ravel()
550 self.coord_lalo = np.array([lat[mask].ravel(), lon[mask].ravel()]).transpose()
552 # 4) read UTM coordinates
553 coord_utm_obj = CoordinatesUTM(file_path=join(dirname(self.file_path), "coordinates_utm.h5"),
554 logger=self.logger)
555 coord_utm_obj.open()
556 self.coord_utm = coord_utm_obj.coord_utm[:, mask].transpose()
558 def createMask(self):
559 """Create a mask.
561 Create a mask in the size of the radar image which is used to read the geometry and SLC data for the selected
562 scatterers.
563 """
564 mask = np.zeros((self.length, self.width), dtype=np.bool_)
565 tmp = [tuple([c[0], c[1]]) for c in self.coord_xy]
566 for i in tmp:
567 mask[i] = True
568 return mask
570 def addPointsFromObj(self, *, new_point_id: np.ndarray, new_coord_xy: np.ndarray, new_phase: np.ndarray,
571 new_num_points: int, input_path: str):
572 """Add new points and their attributes to the existing data.
574 Parameters
575 ----------
576 new_point_id: np.ndarray
577 point_id of the new scatterers.
578 new_coord_xy: np.ndarray
579 radar coordinates of the new scatterers.
580 new_phase: np.ndarray
581 phase of the new scatterers.
582 new_num_points: int
583 number of new points.
584 input_path: str
585 path to input files (slcStack.h5, geometryRadar.h5).
586 """
587 self.point_id = np.append(self.point_id, new_point_id)
588 self.coord_xy = np.append(self.coord_xy, new_coord_xy, axis=0)
589 self.phase = np.append(self.phase, new_phase, axis=0)
590 self.num_points += new_num_points
592 # all data must be ordered, so that all external data can be loaded correctly
593 sort_idx = np.argsort(self.point_id)
594 self.point_id = self.point_id[sort_idx]
595 self.coord_xy = self.coord_xy[sort_idx, :]
596 self.phase = self.phase[sort_idx, :]
597 # refresh by reopening all external data
598 self.openExternalData(input_path=input_path)
600 def removePoints(self, mask: np.ndarray = None, *, keep_id: [np.ndarray, list], input_path: str):
601 """Remove all entries from specified points.
603 The possible options exist for removing the points:
604 a) Keep all points which are set to True in a 'mask' with size (num_points x 1). Or
605 b) Keep all points whose ID is listed in keep_id. The rest of the points will be removed.
607 Parameters
608 ----------
609 mask: np.ndarray
610 mask to select points to be kept, rest will be removed (default: None).
611 keep_id: np.ndarray
612 list of point_id to keep.
613 input_path: str
614 path to input files (slcStack.h5, geometryRadar.h5).
615 """
616 if mask is None:
617 mask = np.ones((self.num_points,), dtype=np.bool_)
618 for p in self.point_id:
619 if p not in keep_id:
620 mask[self.point_id == p] = False
621 self.point_id = self.point_id[mask]
622 self.coord_xy = self.coord_xy[mask, :]
623 self.phase = self.phase[mask, :]
624 self.num_points = mask[mask].shape[0]
625 # refresh by reopening all external data
626 self.openExternalData(input_path=input_path)
629class Network:
630 """Spatial network of PS candidates."""
632 def __init__(self, *, file_path: str, logger: Logger):
633 """Init.
635 Parameters
636 ----------
637 file_path: str
638 absolute path to working directory for creating/loading 'psNetwork.h5'
639 logger: Logger
640 Logging handler.
641 """
642 self.num_arcs = None
643 self.gamma = None
644 self.arcs = None
645 self.slant_range = None
646 self.loc_inc = None
647 self.phase = None
648 self.vel = None
649 self.demerr = None
650 self.ifg_net_obj = None
651 self.width = None
652 self.length = None
653 self.wavelength = None
654 self.file_path = file_path
655 self.logger = logger
657 def writeToFile(self):
658 """Write all existing data to psNetwork.h5 file."""
659 self.logger.info(msg="write data to {}...".format(self.file_path))
661 if exists(self.file_path):
662 os.remove(self.file_path)
664 with h5py.File(self.file_path, 'w') as f:
665 f.attrs["num_arcs"] = self.num_arcs
666 f.create_dataset('arcs', data=self.arcs)
667 f.create_dataset('phase', data=self.phase)
668 f.create_dataset('loc_inc', data=self.loc_inc)
669 f.create_dataset('slant_range', data=self.slant_range)
671 def open(self, *, input_path: str):
672 """Read stored information from existing .h5 file."""
673 with h5py.File(self.file_path, 'r') as f:
674 self.num_arcs = f.attrs["num_arcs"]
675 self.arcs = f["arcs"][:]
676 self.phase = f["phase"][:]
677 self.loc_inc = f["loc_inc"][:]
678 self.slant_range = f["slant_range"][:]
679 self.openExternalData(input_path=input_path)
681 def openExternalData(self, *, input_path: str):
682 """Read data from slcStack.h5 and IfgNetwork.h5 files."""
683 slc_stack_obj = slcStack(join(input_path, "slcStack.h5"))
684 slc_stack_obj.open(print_msg=False)
685 self.wavelength = np.float64(slc_stack_obj.metadata["WAVELENGTH"])
686 self.length = slc_stack_obj.length # y-coordinate axis (azimut)
687 self.width = slc_stack_obj.width # x-coordinate axis (range)
689 # 3) read IfgNetwork
690 self.ifg_net_obj = IfgNetwork()
691 self.ifg_net_obj.open(path=join(dirname(self.file_path), "ifg_network.h5"))
693 def computeArcObservations(self, *, point_obj: Points, arcs: np.ndarray):
694 """Compute the phase observations for each arc.
696 Compute double difference phase observations, i.e. the phase differences for each arc in the network from the
697 phase of the two scatterers connected by the arc.
699 Parameters
700 ----------
701 point_obj: Points
702 object of class Points.
703 arcs: np.ndarray
704 Array with the indices of the points connected by an arc.
705 """
706 self.arcs = arcs
707 self.num_arcs = self.arcs.shape[0]
708 self.logger.info(msg="no. arcs:\t{}".format(self.num_arcs))
710 self.phase = np.zeros((self.num_arcs, point_obj.ifg_net_obj.num_ifgs))
711 self.loc_inc = np.zeros((self.num_arcs,))
712 self.slant_range = np.zeros((self.num_arcs,))
713 for idx, arc in enumerate(self.arcs):
714 self.phase[idx, :] = np.angle(
715 np.exp(1j * point_obj.phase[arc[0], :]) * np.conjugate(np.exp(1j * point_obj.phase[arc[1], :])))
716 self.loc_inc[idx] = np.mean([point_obj.loc_inc[arc[0]], point_obj.loc_inc[arc[1]]])
717 self.slant_range[idx] = np.mean([point_obj.slant_range[arc[0]], point_obj.slant_range[arc[1]]])
719 self.logger.info(msg="ifg arc observations created.")
721 def removeArcs(self, *, mask: np.ndarray):
722 """Remove arcs from the list of arcs in the network.
724 Parameter
725 ---------
726 mask: np.ndarray
727 mask to select arcs to be kept, rest will be removed.
728 """
729 self.demerr = self.demerr[mask]
730 self.vel = self.vel[mask]
731 self.phase = self.phase[mask, :]
732 self.loc_inc = self.loc_inc[mask]
733 self.slant_range = self.slant_range[mask]
734 self.arcs = np.array(self.arcs)
735 self.arcs = self.arcs[mask, :]
736 self.gamma = self.gamma[mask]
737 self.num_arcs = mask[mask].shape[0]
740class NetworkParameter(Network):
741 """Spatial Network with the estimated parameters of each arc in the network."""
743 def __init__(self, *, file_path: str, logger: Logger):
744 """Init."""
745 super().__init__(file_path=file_path, logger=logger)
746 self.gamma = None
747 self.vel = None
748 self.demerr = None
749 self.slant_range = None
750 self.loc_inc = None
751 self.phase = None
752 self.arcs = None
753 self.num_arcs = None
754 self.logger = logger
756 def prepare(self, *, net_obj: Network, demerr: np.ndarray, vel: np.ndarray, gamma: np.ndarray):
757 """Prepare.
759 Parameter
760 -----------
761 net_obj: Network
762 object of class Network.
763 demerr: np.ndarray
764 estimated DEM error for each arc in the network.
765 vel: np.ndarray
766 estimated velocity for each arc in the network.
767 gamma: np.ndarray
768 estimated temporal coherence for each arc in the network.
769 """
770 self.num_arcs = net_obj.num_arcs
771 self.arcs = net_obj.arcs
772 self.phase = net_obj.phase
773 self.loc_inc = net_obj.loc_inc
774 self.slant_range = net_obj.slant_range
775 self.demerr = demerr
776 self.vel = vel
777 self.gamma = gamma
779 def writeToFile(self):
780 """Write DEM error, velocity and temporal coherence to file."""
781 super().writeToFile()
783 with h5py.File(self.file_path, 'r+') as f: # append existing file
784 f.create_dataset('demerr', data=self.demerr)
785 f.create_dataset('vel', data=self.vel)
786 f.create_dataset('gamma', data=self.gamma)
788 def open(self, *, input_path: str):
789 """Read data from file."""
790 super().open(input_path=input_path)
792 with h5py.File(self.file_path, 'r') as f:
793 self.demerr = f["demerr"][:]
794 self.vel = f["vel"][:]
795 self.gamma = f["gamma"][:]