-
Notifications
You must be signed in to change notification settings - Fork 7
/
03_example_random_ifgs.py
executable file
·113 lines (85 loc) · 7.62 KB
/
03_example_random_ifgs.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Sep 18 11:55:59 2020
"""
import numpy as np
import matplotlib.pyplot as plt
import pickle
import sys
from pathlib import Path
import syinterferopy
from syinterferopy.random_generation import create_random_synthetic_ifgs
#%% 0: Things to set
srtm_tools_dir = Path('/home/matthew/university_work/15_my_software_releases/SRTM-DEM-tools-3.0.0') # SRTM-DEM-tools will be needed for this example. It can be downloaded from https://github.com/matthew-gaddes/SRTM-DEM-tools
gshhs_dir = Path("/home/matthew/university_work/data/coastlines/gshhg/shapefile/GSHHS_shp") # coastline information, available from: http://www.soest.hawaii.edu/pwessel/gshhg/
SRTM_dem_settings = {'SRTM1_or3' : 'SRTM3', # 1 arc second (SRTM1) is not yet supported.
'water_mask_resolution' : 'f', # 'c', 'i', 'h' or 'f' (lowest to highest resolution, highest can be very slow)
'SRTM3_tiles_folder' : Path('./SRTM3/'), # folder for DEM tiles.
'download' : True, # If tile is not available locally, try to download it
'void_fill' : True, # some tiles contain voids which can be filled (slow)
'gshhs_dir' : gshhs_dir} # srmt-dem-tools needs access to data about coastlines
synthetic_ifgs_settings = {'defo_sources' : ['no_def', 'dyke', 'sill', 'mogi'], # deformation patterns that will be included in the dataset.
'n_ifgs' : 7, # the number of synthetic interferograms to generate
'n_pix' : 224, # number of 3 arc second pixels (~90m) in x and y direction
'outputs' : ['uuu', 'uud', 'www', 'wwd', 'rid'], # channel outputs. uuu = unwrapped across all 3, uu
'intermediate_figure' : False, # if True, a figure showing the steps taken during creation of each ifg is displayed.
'coh_threshold' : 0.7, # if 1, there are no areas of incoherence, if 0 all of ifg is incoherent.
'min_deformation' : 0.05, # deformation pattern must have a signals of at least this many metres.
'max_deformation' : 0.25, # deformation pattern must have a signal no bigger than this many metres.
'snr_threshold' : 2.0, # signal to noise ratio (deformation vs turbulent and topo APS) to ensure that deformation is visible. A lower value creates more subtle deformation signals.
'turb_aps_mean' : 0.02, # turbulent APS will have, on average, a maximum strenghto this in metres (e.g 0.02 = 2cm)
'turb_aps_length' : 5000} # turbulent APS will be correlated on this length scale, in metres.
volcano_dems = [{'name': 'Vulsini', 'centre': (11.93, 42.6), 'side_length' : (40e3, 40e3)}, # centre is lon then lat, side length is x then y, in km.
{'name': 'Campi Flegrei', 'centre': (14.139, 40.827), 'side_length' : (40e3, 40e3)},
{'name': 'Witori', 'centre': (150.516, -5.576), 'side_length' : (40e3, 40e3)},
{'name': 'Lolobau', 'centre': (151.158, -4.92), 'side_length' : (40e3, 40e3)},
{'name': 'Krakatau', 'centre': (105.423, -6.102), 'side_length' : (40e3, 40e3)},
{'name': 'Batur', 'centre': (115.375, -8.242), 'side_length' : (40e3, 40e3)},
{'name': 'Taal', 'centre': (120.993, 14.002), 'side_length' : (40e3, 40e3)},
{'name': 'Aira', 'centre': (130.657, 31.593), 'side_length' : (40e3, 40e3)},
{'name': 'Asosan', 'centre': (131.104, 32.884), 'side_length' : (40e3, 40e3)},
{'name': 'Naruko', 'centre': (140.734, 38.729), 'side_length' : (40e3, 40e3)},
{'name': 'Towada', 'centre': (140.88, 40.51), 'side_length' : (40e3, 40e3)}]
#%% Import srtm_dem_tools
sys.path.append(str(srtm_tools_dir)) #
import srtm_dem_tools
from srtm_dem_tools.constructing import SRTM_dem_make_batch
#%% 1: Create a list of locations (in this case subaerial volcanoes) to make interferograms for, and make them.
np.random.seed(0) # 0 used in the example
try:
print('Trying to open a .pkl of the DEMs... ', end = '')
with open('example_03_dems.pkl', 'rb') as f:
volcano_dems2 = pickle.load(f)
f.close()
print('Done. ')
except:
print('Failed. Generating them from scratch, which can be slow. ')
ed_username = input(f'Please enter your USGS Earthdata username: ') # needed to download SRTM3 tiles
ed_password = input(f'Please enter your USGS Earthdata password (NB characters will be visible! ): ')
SRTM_dem_settings['ed_username'] = ed_username # append to the dict of dem_settings so it can be passed to SRTM_dem_make quickly.
SRTM_dem_settings['ed_password'] = ed_password
volcano_dems2 = SRTM_dem_make_batch(volcano_dems, **SRTM_dem_settings) # make the DEMS
with open(f'example_03_dems.pkl', 'wb') as f:
pickle.dump(volcano_dems2, f)
print('Saved the dems as a .pkl for future use. ')
#%% 2: Create the synthetic interferograms
X_all, Y_class, Y_loc, Y_source_kwargs = create_random_synthetic_ifgs(volcano_dems2, **synthetic_ifgs_settings)
#%% Plot one interferogram in the different channel formats.
for ifg_n in range(X_all['uuu'].shape[0]):
units = np.array([['rad', 'rad', 'rad', 'rad', 'intensity'],
['rad', 'rad', 'rad', 'rad', 'intensity'],
['rad', 'm', 'rad', 'm', 'm']])
fig1, axes = plt.subplots(3,5, figsize = (14,7))
fig1.subplots_adjust(hspace = 0.4, left = 0.02, right = 0.98, bottom = 0.05, top = 0.9)
fig1.canvas.manager.set_window_title(f'Channel format of data {ifg_n}')
fig1.suptitle(f"Label: {synthetic_ifgs_settings['defo_sources'][int(Y_class[ifg_n,0])]} \n"
f"u:unwrapped d:DEM w:wrapped r:real i:imaginary")
for format_n, key in enumerate(X_all):
axes[0, format_n].set_title(key)
for channel_n in range(3):
image = axes[channel_n, format_n].imshow(X_all[key][ifg_n,:,:,channel_n])
cbar = plt.colorbar(image, ax = axes[channel_n, format_n], orientation = 'horizontal', pad = 0.2 )
cbar.ax.set_xlabel(units[channel_n, format_n])
if format_n == 0:
axes[channel_n, format_n].set_ylabel(f'Channel {channel_n}')