Skip to content

Commit

Permalink
Lightly tested version with ability to make random time series.
Browse files Browse the repository at this point in the history
  • Loading branch information
@matthew-gaddes
matthew-gaddes committed Aug 25, 2023
1 parent b3065e4 commit 1f5fab9
Show file tree
Hide file tree
Showing 5 changed files with 179 additions and 154 deletions.
141 changes: 5 additions & 136 deletions 05_example_random_ts.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,53 +21,7 @@
import pdb

import syinterferopy
from syinterferopy.random_generation import create_random_synthetic_ifgs

#%% Debug functions


def matrix_show(matrix, title=None, ax=None, fig=None, save_path = None, vmin0 = False):
"""Visualise a matrix
Inputs:
matrix | r2 array or masked array
title | string
ax | matplotlib axes
save_path | string or None | if a string, save as a .png in this location.
vmin0 | boolean |
2017/10/18 | update so can be passed an axes and plotted in an existing figure
2017/11/13 | fix bug in how colorbars are plotted.
2017/12/01 | fix bug if fig is not None
"""
import matplotlib.pyplot as plt
import numpy as np

if ax is None:
fig, ax = plt.subplots()
matrix = np.atleast_2d(matrix) # make at least 2d so can plot column/row vectors

if isinstance(matrix[0,0], np.bool_): # boolean arrays will plot, but mess up the colourbar
matrix = matrix.astype(int) # so convert

if vmin0:
matrixPlt = ax.imshow(matrix,interpolation='none', aspect='auto', vmin = 0)
else:
matrixPlt = ax.imshow(matrix,interpolation='none', aspect='auto')
fig.colorbar(matrixPlt,ax=ax)
if title is not None:
ax.set_title(title)
fig.canvas.set_window_title(f"{title}")

if save_path is not None: # possibly save the figure
if title is None: # if not title is supplied, save with a default name
fig.savefig(f"{save_path}/matrix_show_output.png")
else:
fig.savefig(f"{save_path}/{title}.png") # or with the title, if it's supplied


plt.pause(1) # to force it to be shown when usig ipdb

from syinterferopy.random_generation import create_random_ts_1_volc


#%% 0: Things to set
Expand All @@ -82,17 +36,6 @@ def matrix_show(matrix, title=None, ax=None, fig=None, save_path = None, vmin0 =
'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': 'Campi Flegrei', 'centre': (14.139, 40.827), 'side_length' : (40e3, 40e3)}] # centre is lon then lat, side length is x then y, in km.
Expand Down Expand Up @@ -138,85 +81,11 @@ def matrix_show(matrix, title=None, ax=None, fig=None, save_path = None, vmin0 =



#%%

# one function to call on each volcano
# dem should be large.
# trasnalte def and dem.
# random time series (x1000)
#

def create_random_ts_1_volc(dem_dict, n_pix = 224, d_start = "20141231", d_stop = "20230801",
n_def_location = 2, n_tcs = 5, n_atms = 5,
topo_aps_mean = 56.0, topo_aps_var = 2.0, turb_aps_mean = 0.02):
"""Create random time series of Sentinel-1 data for a given volcano.
n_def_location random crops of the DEM with the deformation placed somewhere randomly are made (the majority of the deformation must be on land. )
n_tcs patterns for the temporal evolution of the deformation are made.
n_atms instances of atmospheric noise are made (turbulent and topographically correlated)
Total number of time series = n_def_location * n_tcs * n_atms.
Inputs:
dem_dict | dict | DEM and vairous parameters.
n_pix | int | size out interferograms outp.
d_start | string | YYYYMMDD of time series start
d_stop | string | YYYYMMDD of time series end
n_def_locations | int | Number of random DEM crops with random deformation placement.
n_tcs | int | Number of random deformation time courses.
n_atms | int | Number of random time series atmospheres.
topo_aps_mean | float | rad/km of delay for the topographically correlated APS
topo_aps_var | float | rad/km. Sets the strength difference between topographically correlated APSs
turb_aps_mean | float | mean strength of turbulent atmospheres, in metres. Note the the atmosphere_turb funtion takes cmm, and the value provided in m is converted first
Returns:
Saves pickle file of the time series, the deformaiton time course, and hte DEM.
History:
2023_08_24 | MEG | Written.
"""


from syinterferopy.random_generation import generate_dems_and_defos
from syinterferopy.atmosphere import atmosphere_turb, atmosphere_topos
from syinterferopy.temporal import defo_to_ts, generate_random_temporal_baselines, sample_deformation_on_acq_dates
from syinterferopy.temporal import generate_uniform_temporal_baselines, generate_random_tcs

out_file_n = 0

defos, dems = generate_dems_and_defos(dem_dict, n_pix, min_deformation = 0.05, max_deformation = 0.25, n_def_location = n_def_location) # generate multiple crops of DEM with deformation places randomly on it.
tcs, def_dates = generate_random_tcs(n_tcs, d_start, d_stop) # generate random time courses for the deformation.


for defo_n, defo in enumerate(defos):
for tc_n, tc in enumerate(tcs):
for atm_n in range(n_atms):
#acq_dates, tbaselines = generate_random_temporal_baselines(d_start, d_stop)
acq_dates, tbaselines = generate_uniform_temporal_baselines(d_start, d_stop)
tc_resampled = sample_deformation_on_acq_dates(acq_dates, tc, def_dates) # get the deformation on only the days when there's a satellite acquisitions. Note that this is 0 on the first date.
defo_ts = defo_to_ts(defo, tc_resampled)

atm_turbs = atmosphere_turb(len(acq_dates)-1, dem_dict['lons_mg'][:n_pix,:n_pix], dem_dict['lats_mg'][:n_pix,:n_pix], mean_m = turb_aps_mean) # generate some random atmos
atm_topos = atmosphere_topos(len(acq_dates)-1, dems[defo_n,], topo_aps_mean, topo_aps_var)

ts = defo_ts + atm_turbs + atm_topos
if atm_n == 0:
tss = ma.zeros((n_atms, ts.shape[0], ts.shape[1], ts.shape[2])) # initiliase
tss[atm_n,] = ts
print(f"Deformation location: {defo_n} Time course: {tc_n} Atmosphere {atm_n}")


with open(f'file_{out_file_n:05d}.pkl', 'wb') as f:
pickle.dump(tss, f)
pickle.dump(tc, f)
pickle.dump(dems[defo_n], f)
out_file_n += 1



#%% 2: Make time series for that volcano.

create_random_ts_1_volc(volcano_dems2[0], n_def_location = 10, n_pix = 224)
create_random_ts_1_volc(outdir = Path('./05_example_outputs/'), dem_dict = volcano_dems2[0], n_pix = 224, d_start = "20141231", d_stop = "20230801",
n_def_location = 2, n_tcs = 2, n_atms = 2,
topo_delay_var = 0.00005, turb_aps_mean = 0.02)



Expand Down
74 changes: 58 additions & 16 deletions syinterferopy/atmosphere.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,36 +11,76 @@

#%%

def atmosphere_topos(n_ifgs, dem_m, strength_mean = 56.0, strength_var = 2.0):
""" Create multiple topographically correlated APS. See function below that is called to make each one.
def atmosphere_topos(n_ifgs, dem_m, delay_grad = -0.0003, delay_var = 0.00005, zero_delay_h = 10000.):
""" Create multiple topographically correlated APS for each acquisition.
Different to the function below that makes them for . See function below that is called to make each one.
Inputs:
n_ifgs | int | Number of interferogram atmospheres to generate.
strength_mean | float | rad/km of delay. default is 56.0, taken from Fig5 Pinel 2011 (Statovolcanoes...)
strength_var | float | variance of rad/km delay. Default is 2.0, which gives values similar to Fig 5 of above.
delay_grad | float | absolute delay gradient. Units are m delay / m height. With -0.0003, ~3000m elevation required for 1m of delay.
delay_var | float | variance of the delays. Higher value gives more variety in atmospheres, so bigger interfergoram atmospheres.
zero_delay_h | float | Height at which delays tend to zero. Not tested.
Returns:
atm_topos | r3 ma | interferograms with mask applied.
History:
2023_08_24 | MEG | Written
E.g. GACOS for vesuvius: Vesuvius 1281m high, 2.1m of delay.
sea level: (0m) 2.45
"""

import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt

n_atms = n_ifgs + 1 # n_ifgs is 1 more than n_acqs or n_atms

for n_atm in range(n_atms): # loop through making atms. Make n_acqs of them, which is 1 more than n_ifgs
atm_topo = atmosphere_topo(dem_m, strength_mean, strength_var, difference = False) #
# sanity check of how delay is defined.
# grad_av = (0 - 1281) / (2.45 - 2.1) # height on the y axis, delay on the x (as per Bekaert 2015). Makes sense as heiht is normally vertical.
# c = 10000 # height at which delays are 0, in m. Also called zero_delay_height.
# # heights = grad_av * delays + c
# # heights - c / grad_av = delays
# heights = np.arange(0, 10000, 1)
# delays = (heights - c) / grad_av
# f, ax = plt.subplots()
# ax.plot(delays, heights)


n_atms = n_ifgs + 1 # n_ifgs is 1 more than n_acqs or n_atms
delay_grads = delay_grad + (delay_var * np.random.randn(n_atms))

# # debug figure
# f, ax = plt.subplots(1,1)
# ax.set_xlabel('Delay (m)')
# ax.set_ylabel('Height (m)')
# heights = np.arange(0, 10000, 1)
# for delay_grad in delay_grads:
# delays = (heights - zero_delay_h) / (1 / delay_grad)
# ax.plot(delays, heights)

if n_atm == 0: # if the first time.
for n_atm in range(n_atms): # loop through making atms. Make n_acqs of them, which is 1 more than n_ifgs
atm_topo = (dem_m - zero_delay_h) / (1/delay_grads[n_atm])
#atm_topo = (dem_m) / (1/delay_grads[n_atm])
if n_atm == 0: # if the first time.
atm_topos = ma.zeros((n_atms, atm_topo.shape[0], atm_topo.shape[1])) # initialise array to store results
atm_topos[n_atm, :, :] = atm_topo #

atm_topos_diff = ma.diff(atm_topos, axis = 0) # atmoshperes in ifgs are the difference of two atmospheres.
# f, axes = plt.subplots(2,8)
# for i in range(8):
# axes[0,i].matshow(atm_topos[i,], vmin = np.min(atm_topos), vmax = np.max(atm_topos))
# axes[1,i].matshow(atm_topos_diff[i,], vmin = np.min(atm_topos_diff), vmax = np.max(atm_topos_diff))
dem_minimum_arg = np.unravel_index(np.argmin(dem_m), dem_m.shape) # find the lowest pixel in the DEM
atm_topo_offsets = atm_topos[:, dem_minimum_arg[0], dem_minimum_arg[1]] # use this as a reference pixel for the atmoshperes.

for atm_n, atm_topo in enumerate(atm_topos): # loop through the atmospheres
atm_topo -= atm_topo_offsets[atm_n] # and reference them.

atm_topos_diff = ma.diff(atm_topos, axis = 0) # atmoshperes in ifgs are the difference of two atmospheres.

# # debug figure
# f, axes = plt.subplots(3,n_atm)
# for i in range(n_atm):
# axes[0,i].plot(ma.ravel(atm_topos[i,]), ma.ravel(dem_m))
# axes[0,i].set_xlim(left = ma.min(atm_topos))
# axes[1,i].matshow(atm_topos[i,], vmin = np.min(atm_topos), vmax = np.max(atm_topos))
# axes[2,i].matshow(atm_topos_diff[i,], vmin = np.min(atm_topos_diff), vmax = np.max(atm_topos_diff))

return atm_topos_diff

#%%
Expand All @@ -62,11 +102,13 @@ def atmosphere_topo(dem_m, strength_mean = 56.0, strength_var = 2.0, difference
import numpy as np
import numpy.ma as ma

envisat_lambda = 0.056 #envisat/S1 wavelength in m
satellite_lambda = 0.055465763 # S1 wavelength in m
#satellite_lambda = 0.056 #envisat wavelength in m
dem = 0.001 * dem_m # convert from metres to km

if difference is False:
ph_topo = (strength_mean + strength_var * np.random.randn(1)) * dem
print(f"The use of this function with difference = False is discouraged. Use atmosphere_topos instead (and not atmosphere_topo)")
elif difference is True:
ph_topo_aq1 = (strength_mean + strength_var * np.random.randn(1)) * dem # this is the delay for one acquisition
ph_topo_aq2 = (strength_mean + strength_var * np.random.randn(1)) * dem # and for another
Expand All @@ -77,10 +119,10 @@ def atmosphere_topo(dem_m, strength_mean = 56.0, strength_var = 2.0, difference
import sys; sys.exit()

# convert from rad to m
ph_topo_m = (ph_topo / (4*np.pi)) * envisat_lambda # delay/elevation ratio is taken from a paper (pinel 2011) using Envisat data
ph_topo_m = (ph_topo / (4*np.pi)) * satellite_lambda # delay/elevation ratio is taken from a paper (pinel 2011) using Envisat data


if np.max(ph_topo_m) < 0: # ensure that it always start from 0, either increasing or decreasing
if np.max(ph_topo_m) < 0: # ensure that it always starts from 0, either increasing or decreasing
ph_topo_m -= np.max(ph_topo_m)
else:
ph_topo_m -= np.min(ph_topo_m)
Expand Down
20 changes: 20 additions & 0 deletions syinterferopy/aux.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,26 @@
@author: matthew
"""

import pdb

#%%

def rescale_defo(defos, magnitude = 1.):
""" Given deformation patterns in metres, rescale so that the maximum value is "magnitude".
Inputs:
defos | r3 ma | n_times x ny x nx, water masked.
magnitude | float | new maximum signals
Retrurns:
defos | r3 ma | As above, with new magnitude.
History:
2023_08_25 | MEG | Written.
"""
import numpy.ma as ma
defo_maxs = ma.max(defos, axis = (1,2))
for defo_n, defo in enumerate(defos):
defo /= defo_maxs[defo_n]
return defos

#%%

def lon_lat_to_ijk(lons_mg, lats_mg):
Expand Down
Loading

0 comments on commit 1f5fab9

Please sign in to comment.