Coverage for sarvey/sarvey_plot.py: 0%
236 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"""Plot module for SARvey."""
31import argparse
32import time
33import os
34from os.path import join, basename, dirname
35import matplotlib
36import matplotlib.pyplot as plt
37from matplotlib import colormaps
38import numpy as np
39import logging
40from logging import Logger
41import sys
43from mintpy.utils import ptime
44from mintpy.objects.colors import ColormapExt
45from mintpy.utils.plot import auto_flip_direction
47from sarvey.ifg_network import IfgNetwork
48from sarvey.objects import Points, AmplitudeImage, BaseStack
49from sarvey import console
50from sarvey import viewer
51from sarvey.config import loadConfiguration
52import sarvey.utils as ut
54try:
55 matplotlib.use('QtAgg')
56except ImportError as e:
57 print(e)
59EXAMPLE = """Example:
60 sarvey_plot outputs/p2_coh60_ts.h5 -t # plot average velocity and time series
61 sarvey_plot outputs/p2_coh80_ts.h5 -m -a # plot velocity map and DEM correction interactively
62 sarvey_plot outputs/p2_coh80_ts.h5 -r -n 0 5 # plot residuals for image 0 to 5
63 sarvey_plot outputs/p2_coh80_ifg_wr.h5 -p -n 0 1 -a # plot wrapped phase of final point selection for interferogram 0
64 sarvey_plot outputs/p1_ifg_wr.h5 -p -n 0 1 -a # plot wrapped phase of the first order network
65 sarvey_plot -i -a outputs/ifg_stack.h5 # interactively plot interferograms
66 sarvey_plot -i outputs/ifg_stack.h5 # store interferograms as png files
67 [...]
68"""
71def plotMap(*, obj_name: str, save_path: str, interactive: bool = False, input_path: str, logger: Logger):
72 """Plot results from sarvey as map in radar coordinates.
74 Plot the velocity map, DEM error, squared sum of residuals, temporal coherence and spatiotemporal consistency.
76 Parameters
77 ----------
78 obj_name : str
79 Path to the Points object file.
80 save_path : str
81 Path to the directory where the figures are saved.
82 interactive : bool
83 If True, the plots will be shown interactively.
84 input_path : str
85 Path to the inputs directory containing slcStack.h5 and geometryRadar.h5.
86 logger : Logger
87 Logger object.
88 """
89 console.showLogoSARvey(logger=logger, step="Plot map")
91 scatter_size = 1
93 point_obj = Points(file_path=obj_name, logger=logger)
94 point_obj.open(input_path=input_path)
96 bmap_obj = AmplitudeImage(file_path=join(dirname(obj_name), "background_map.h5"))
97 vel, demerr, _, coherence, omega, v_hat = ut.estimateParameters(obj=point_obj, ifg_space=False)
99 ax = bmap_obj.plot(logger=logger)
100 sc = ax.scatter(point_obj.coord_xy[:, 1], point_obj.coord_xy[:, 0], c=demerr, s=scatter_size,
101 cmap=colormaps["jet_r"])
102 plt.colorbar(sc, label="[m]", pad=0.03, shrink=0.5)
103 plt.title("DEM correction")
104 plt.ylabel('Azimuth')
105 plt.xlabel('Range')
106 plt.tight_layout()
107 plt.gcf().savefig(join(save_path, "map_dem_correction.png"), dpi=300)
108 if interactive:
109 plt.show()
110 else:
111 plt.close(plt.gcf())
113 ax = bmap_obj.plot(logger=logger)
114 sc = ax.scatter(point_obj.coord_xy[:, 1], point_obj.coord_xy[:, 0], c=omega, s=scatter_size,
115 cmap=colormaps["autumn_r"])
116 plt.colorbar(sc, label="", pad=0.03, shrink=0.5)
117 plt.title("Squared sum of residuals")
118 plt.ylabel('Azimuth')
119 plt.xlabel('Range')
120 plt.tight_layout()
121 plt.gcf().savefig(join(save_path, "map_squared_sum_of_residuals.png"), dpi=300)
122 if interactive:
123 plt.show()
124 else:
125 plt.close(plt.gcf())
127 v_range = np.max(np.abs(vel * 100))
129 ax = bmap_obj.plot(logger=logger)
130 sc = ax.scatter(point_obj.coord_xy[:, 1], point_obj.coord_xy[:, 0], c=vel * 100, s=scatter_size,
131 cmap=colormaps["jet_r"],
132 vmin=-v_range, vmax=v_range)
133 plt.colorbar(sc, label="[cm / year]", pad=0.03, shrink=0.5)
134 plt.title("Mean Velocity")
135 plt.ylabel('Azimuth')
136 plt.xlabel('Range')
137 plt.tight_layout()
138 plt.gcf().savefig(join(save_path, "map_velocity.png"), dpi=300)
139 if interactive:
140 plt.show()
141 else:
142 plt.close(plt.gcf())
144 ax = bmap_obj.plot(logger=logger)
145 sc = ax.scatter(point_obj.coord_xy[:, 1], point_obj.coord_xy[:, 0], c=coherence, vmin=0, vmax=1, s=scatter_size,
146 cmap=colormaps["autumn"])
147 plt.colorbar(sc, label="[-]", pad=0.03, shrink=0.5)
148 plt.title("Temporal coherence")
149 plt.ylabel('Azimuth')
150 plt.xlabel('Range')
151 plt.tight_layout()
152 plt.gcf().savefig(join(save_path, "map_coherence.png"), dpi=300)
153 if interactive:
154 plt.show()
155 else:
156 plt.close(plt.gcf())
158 stc = ut.spatiotemporalConsistency(coord_utm=point_obj.coord_utm, phase=point_obj.phase,
159 wavelength=point_obj.wavelength,
160 min_dist=50, max_dist=np.inf, knn=40)
162 ax = bmap_obj.plot(logger=logger)
163 sc = ax.scatter(point_obj.coord_xy[:, 1], point_obj.coord_xy[:, 0], c=stc * 100, s=scatter_size,
164 cmap=colormaps["autumn_r"])
165 plt.colorbar(sc, label="[cm]", pad=0.03, shrink=0.5)
166 plt.title("Spatiotemporal consistency")
167 plt.ylabel('Azimuth')
168 plt.xlabel('Range')
169 plt.tight_layout()
170 plt.gcf().savefig(join(save_path, "map_spatiotemporal_consistency.png"), dpi=300)
171 if interactive:
172 plt.show()
173 else:
174 plt.close(plt.gcf())
177def plotTS(*, obj_name: str, input_path: str, logger: Logger):
178 """Plot the derived displacement time series.
180 Parameters
181 ----------
182 obj_name : str
183 Path to the Points object file.
184 input_path : str
185 Path to the inputs directory containing slcStack.h5 and geometryRadar.h5.
186 logger : Logger
187 Logger object.
188 """
189 console.showLogoSARvey(logger=logger, step="Plot time series")
191 point_obj = Points(file_path=obj_name, logger=logger)
192 point_obj.open(input_path=input_path)
193 if point_obj.phase.shape[1] == point_obj.ifg_net_obj.num_ifgs:
194 logger.warning(msg="File contains ifg phase and not phase time series. Cannot display.")
195 else:
196 viewer.TimeSeriesViewer(point_obj=point_obj, logger=logger, input_path=input_path)
197 plt.show()
200def plotPhase(*, obj_name: str, save_path: str, image_range: tuple, interactive: bool = False, input_path: str,
201 logger: Logger):
202 """Plot the phase of a Points object file in geographic coordinates (WGS84).
204 Plot the phase of each interferogram or each image depending on the domain of the input file.
206 Parameters
207 ----------
208 obj_name : str
209 Path to the Points object file.
210 save_path : str
211 Path to the directory where the figures are saved.
212 image_range : tuple
213 Range of images to be plotted.
214 interactive : bool
215 If True, the plots will be shown interactively.
216 input_path : str
217 Path to the inputs directory containing slcStack.h5 and geometryRadar.h5.
218 logger : Logger
219 Logger object.
220 """
221 console.showLogoSARvey(logger=logger, step="Plot phase")
223 point_obj = Points(file_path=obj_name, logger=logger)
224 point_obj.open(input_path=input_path)
226 if image_range is None:
227 viewer.plotIfgs(phase=point_obj.phase, coord=point_obj.coord_lalo, ttl="Phase")
228 else:
229 viewer.plotIfgs(phase=point_obj.phase[:, image_range[0]:image_range[1]], coord=point_obj.coord_lalo,
230 ttl="Phase")
232 plt.gcf().savefig(join(save_path, "{}_phase.png".format(basename(obj_name)[:-3])), dpi=300)
234 if interactive:
235 plt.show()
238def plotResidualPhase(*, obj_name: str, save_path: str, image_range: tuple, interactive: bool = False,
239 input_path: str, logger: Logger):
240 """Plot the residual phase of a Points object file in geographic coordinates (WGS84).
242 The residuals are derived by substracting the phase contributions based on the estimated parameters.
244 Parameters
245 ----------
246 obj_name : str
247 Path to the Points object file.
248 save_path : str
249 Path to the directory where the figures are saved.
250 image_range : tuple
251 Range of images to be plotted.
252 interactive : bool
253 If True, the plots will be shown interactively.
254 input_path : str
255 Path to the inputs directory containing slcStack.h5 and geometryRadar.h5.
256 logger : Logger
257 Logger object.
258 """
259 console.showLogoSARvey(logger=logger, step="Plot residual phase")
261 point_obj = Points(file_path=obj_name, logger=logger)
262 point_obj.open(input_path=input_path)
264 if point_obj.phase.shape[1] == point_obj.ifg_net_obj.num_ifgs:
265 v_hat = ut.estimateParameters(obj=point_obj, ifg_space=True)[-1]
266 else:
267 v_hat = ut.estimateParameters(obj=point_obj, ifg_space=False)[-1]
269 if image_range is None:
270 viewer.plotIfgs(phase=v_hat, coord=point_obj.coord_lalo, ttl="Residual phase")
271 else:
272 viewer.plotIfgs(phase=v_hat[:, image_range[0]:image_range[1]], coord=point_obj.coord_lalo, ttl="Residual phase")
274 plt.gcf().savefig(join(save_path, "{}_residual_phase.png".format(basename(obj_name)[:-3])), dpi=300)
276 if interactive:
277 plt.show()
280def plotAllIfgs(*, obj_name: str, save_path: str, interactive: bool = False, logger: Logger):
281 """Plot all interferograms inside the ifg_stack.h5 file.
283 If interactivity is enabled, the plots are shown and figures are not saved. Otherwise, the figures are
284 not shown but saved as png files.
285 If the ifg_network.h5 file is available, the baselines are displayed in the title of each interferogram.
287 Parameters
288 ----------
289 obj_name : str
290 Path to the ifg_stack.h5 file.
291 save_path : str
292 Path to the directory where the figures are saved.
293 interactive : bool
294 If True, the plots will be shown interactively.
295 logger : Logger
296 Logger object.
297 """
298 console.showLogoSARvey(logger=logger, step="Plot interferograms")
300 if obj_name.split("/")[-1] != "ifg_stack.h5":
301 logger.warning(msg="Cannot plot ifgs from {}".format(obj_name))
302 return
304 ifg_stack_obj = BaseStack(file=obj_name, logger=logger)
305 ifgs = ifg_stack_obj.read(dataset_name="ifgs")
307 path_ifg_net = join(dirname(obj_name), "ifg_network.h5")
308 if os.path.exists(path_ifg_net):
309 ifg_net_obj = IfgNetwork()
310 ifg_net_obj.open(path=path_ifg_net)
311 else:
312 logger.warning(msg="'ifg_network.h5' is not available in the same directory as 'ifg_stack.h5'. "
313 "No baseline information available.")
314 ifg_net_obj = None
316 num_ifgs = ifgs.shape[2]
318 prog_bar = ptime.progressBar(maxValue=num_ifgs)
319 start_time = time.time()
320 logger.info(msg="plot and save figures of ifgs.")
321 for i in range(num_ifgs):
322 fig = plt.figure(figsize=[15, 5])
323 ax = fig.add_subplot()
324 ifg = np.angle(ifgs[:, :, i])
325 ifg[ifg == 0] = np.nan
326 im = plt.imshow(ifg, cmap=ColormapExt('cmy').colormap, interpolation='nearest', vmin=-np.pi, vmax=np.pi)
327 auto_flip_direction(ifg_stack_obj.metadata, ax=ax, print_msg=False)
328 ax.set_xlabel("Range")
329 ax.set_ylabel("Azimuth")
330 plt.colorbar(im, ax=ax, label="[rad]", pad=0.03, shrink=0.5)
331 if ifg_net_obj is not None:
332 if ifg_net_obj.dates is not None:
333 date1 = ifg_net_obj.dates[ifg_net_obj.ifg_list[i][0]]
334 date2 = ifg_net_obj.dates[ifg_net_obj.ifg_list[i][1]]
335 ttl = "{date1} - {date2}\nbaselines: {tbase} days, {pbase} m".format(
336 date1=date1,
337 date2=date2,
338 tbase=int(np.round(ifg_net_obj.tbase_ifg[i] * 365.25)),
339 pbase=int(np.round(ifg_net_obj.pbase_ifg[i]))
340 )
341 else:
342 ttl = "baselines: {tbase} days, {pbase} m".format(
343 tbase=int(np.round(ifg_net_obj.tbase_ifg[i] * 365.25)),
344 pbase=int(np.round(ifg_net_obj.pbase_ifg[i]))
345 )
346 plt.title(ttl)
347 plt.tight_layout()
348 if interactive:
349 plt.show()
350 else:
351 fig.savefig(join(save_path, "{}_ifg".format(i)), dpi=300)
352 plt.close(fig)
353 prog_bar.update(value=i + 1, every=1, suffix='{}/{} ifgs'.format(i + 1, num_ifgs))
354 prog_bar.close()
355 m, s = divmod(time.time() - start_time, 60)
356 logger.debug(msg='time used: {:02.0f} mins {:02.1f} secs.'.format(m, s))
359def createParser():
360 """Create_parser."""
361 parser = argparse.ArgumentParser(
362 description='Plot results from MTI\n\n',
363 formatter_class=argparse.RawTextHelpFormatter,
364 epilog=EXAMPLE)
366 parser.add_argument('input_file', help='Path to the input file')
368 parser.add_argument('-t', '--plot-ts', default=False, dest="plotTS", action="store_true",
369 help='Creates an interactive time series viewer.')
371 parser.add_argument('-p', '--plot-phase', default=False, dest="plotPhase", action="store_true",
372 help='Plots the phase.')
374 parser.add_argument('-r', '--plot-residual-phase', default=False, dest="plot_res_phase", action="store_true",
375 help='Plots the residual phase after substracting known components.')
377 parser.add_argument('-m', '--plot-map', default=False, dest="plotMap", action="store_true",
378 help='Plots the velocity map and DEM correction.')
380 parser.add_argument('-i', '--plot-all-ifgs', default=False, dest="plotAllIfgs", action="store_true",
381 help='Plots all ifgs.')
383 parser.add_argument('-n', '--image_range', default=None, dest="image_range", nargs=2, type=int,
384 help='Reduces the number of phase images to the given range. Has no effect on -m and -t. Tuple.'
385 '(default: all images).')
387 parser.add_argument('-a', '--interactive', default=False, dest="interactive", action="store_true",
388 help='Enables interactive visualisation of figures. Is always ON for plot-ts.')
390 parser.add_argument('-w', '--workdir', default=None, dest="workdir",
391 help='Working directory (default: current directory).')
393 return parser
396def main(iargs=None):
397 """Run Main."""
398 parser = createParser()
399 args = parser.parse_args(iargs)
401 if args.workdir is None:
402 args.workdir = os.getcwd()
404 # initiate logger
405 logging_level = logging.getLevelName('DEBUG')
407 log_format = logging.Formatter('%(asctime)s - %(levelname)s - %(message)s')
408 logger = logging.getLogger(__name__)
410 current_datetime = time.strftime("%Y-%m-%d-%H-%M-%S", time.localtime())
411 log_filename = f"sarvey_plot_{current_datetime}.log"
412 if not os.path.exists(os.path.join(os.getcwd(), "logfiles")):
413 os.mkdir(os.path.join(os.getcwd(), "logfiles"))
414 file_handler = logging.FileHandler(filename=os.path.join(os.getcwd(), "logfiles", log_filename))
415 file_handler.setFormatter(log_format)
416 logger.addHandler(file_handler)
418 console_handler = logging.StreamHandler(sys.stdout)
419 console_handler.setFormatter(log_format)
420 logger.addHandler(console_handler)
421 logger.setLevel(logging_level)
423 logger.info("Working directory: {}".format(args.workdir))
424 args.input_file = join(args.workdir, args.input_file)
426 config_file_path = os.path.abspath(join(args.workdir, dirname(args.input_file), "config.json"))
428 if not os.path.exists(config_file_path):
429 # check if any config file is available in upper directory (backward compatibility)
430 files = np.array([os.path.abspath(f) for f in os.listdir(join(dirname(config_file_path), ".."))
431 if os.path.isfile(f)])
432 potential_configs = np.array([(basename(f).split(".")[-1] == "json") and ("config" in basename(f))
433 for f in files])
434 if potential_configs[potential_configs].shape[0] == 0:
435 raise FileNotFoundError(f"Backup configuration file not found: {config_file_path}!")
436 else:
437 logger.warning(msg=f"Backup configuration file not found: {config_file_path}!")
438 logger.warning(msg=f"Other configuration files automatically detected: {files[potential_configs]}!")
439 logger.warning(msg=f"Automatically selected configuration file: {files[potential_configs][0]}!")
440 config_file_path = files[potential_configs][0]
442 config = loadConfiguration(path=config_file_path)
444 folder_name = "p1" if "p1" in basename(args.input_file) else basename(args.input_file)[:5]
445 folder_name = "ifgs" if "ifg_stack" in basename(args.input_file) else folder_name
447 save_path = join(dirname(args.input_file), "pic", folder_name)
448 if not os.path.exists(save_path):
449 if not args.plotTS: # not needed for interactive time series
450 os.mkdir(save_path)
452 selected = False
453 if args.plotTS:
454 # todo: read input_path from config file in same directory as file to be able to load height from geometryRadar
455 plotTS(obj_name=args.input_file, input_path=config.general.input_path, logger=logger)
456 selected = True
458 if args.plotPhase:
459 plotPhase(obj_name=args.input_file, save_path=save_path, image_range=args.image_range,
460 interactive=args.interactive, input_path=config.general.input_path, logger=logger)
461 selected = True
463 if args.plot_res_phase:
464 plotResidualPhase(obj_name=args.input_file, save_path=save_path, image_range=args.image_range,
465 interactive=args.interactive, input_path=config.general.input_path, logger=logger)
466 selected = True
468 if args.plotMap:
469 plotMap(obj_name=args.input_file, save_path=save_path, interactive=args.interactive,
470 input_path=config.general.input_path, logger=logger)
471 selected = True
473 if args.plotAllIfgs:
474 plotAllIfgs(obj_name=args.input_file, save_path=save_path, interactive=args.interactive, logger=logger)
475 selected = True
477 if not selected:
478 logger.info(msg="No action chosen.")
480 # close log-file to avoid problems with deleting the files
481 if logger.hasHandlers():
482 for handler in logger.handlers[:]:
483 logger.removeHandler(handler)
484 handler.flush()
485 handler.close()
488if __name__ == '__main__':
489 main()