diff --git a/05_example_random_ts.py b/05_example_random_ts.py index 018a4c2..668c78b 100644 --- a/05_example_random_ts.py +++ b/05_example_random_ts.py @@ -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 @@ -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. @@ -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) diff --git a/syinterferopy/atmosphere.py b/syinterferopy/atmosphere.py index 14f82e7..789206d 100644 --- a/syinterferopy/atmosphere.py +++ b/syinterferopy/atmosphere.py @@ -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 #%% @@ -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 @@ -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) diff --git a/syinterferopy/aux.py b/syinterferopy/aux.py index 175e0e1..c18e314 100755 --- a/syinterferopy/aux.py +++ b/syinterferopy/aux.py @@ -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): diff --git a/syinterferopy/random_generation.py b/syinterferopy/random_generation.py index ac1ed41..29ff7b0 100755 --- a/syinterferopy/random_generation.py +++ b/syinterferopy/random_generation.py @@ -6,6 +6,97 @@ @author: matthew """ +import pdb + +#%% + + +def create_random_ts_1_volc(outdir, dem_dict, n_pix = 224, d_start = "20141231", d_stop = "20230801", + n_def_location = 10, n_tcs = 10, n_atms = 10, + topo_delay_var = 0.00005, 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_delay_var | float | variance of the delays. Higher value gives more variety in atmospheres, so bigger interfergoram atmospheres. + 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. + + """ + import numpy.ma as ma + import pickle + + 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 + from syinterferopy.aux import rescale_defo + + out_file_n = 0 + + defos_m, 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. + defos = rescale_defo(defos_m, magnitude = 1) # rescale deformation so maximum is awlays 1 + tcs, def_dates = generate_random_tcs(n_tcs, d_start, d_stop, min_def_rate = 1., max_def_rate = 2.) # generate random time courses for the deformation. + volc_outdir_name = dem_dict['name'].replace(' ', '_') + (outdir / volc_outdir_name).mkdir(parents = True, exist_ok = True) + + for defo_n, defo in enumerate(defos): + for tc_n, tc in enumerate(tcs): + #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) # time series of just the deformation. ((n-acq) - 1 x ny x nx) + + for atm_n in range(n_atms): + 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 atmospheres, turbulent. + atm_topos = atmosphere_topos(len(acq_dates)-1, dems[defo_n,], delay_grad = -0.0003, delay_var = topo_delay_var) # generate some random atmospheres, topographically correlated + 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}") + + # possible debug figure. + # n_cols = 10 + # f, axes = plt.subplots(4,n_cols) + # all_signals = np.concatenate((defo_ts[np.newaxis], atm_turbs[np.newaxis], atm_topos[np.newaxis], ts[np.newaxis]), axis = 0) + # signal_max = ma.max(all_signals) + # signal_min = ma.min(all_signals) + + # for col_n in range(n_cols): + # axes[0, col_n].matshow(defo_ts[col_n,:,:], vmin = signal_min, vmax = signal_max) + # axes[1, col_n].matshow(atm_turbs[col_n,:,:], vmin = signal_min, vmax = signal_max) + # axes[2, col_n].matshow(atm_topos[col_n,:,:], vmin = signal_min, vmax = signal_max) + # axes[3, col_n].matshow(ts[col_n,:,:], vmin = signal_min, vmax = signal_max) + # for ax in np.ravel(axes): + # ax.set_xticks([]) + # ax.set_yticks([]) + + with open(outdir / volc_outdir_name / f'defo_{defo_n:06d}_tc_{tc_n:06d}.pkl', 'wb') as f: + pickle.dump(tss, f) + pickle.dump(tc, f) + pickle.dump(dems[defo_n], f) + out_file_n += 1 + + + #%% diff --git a/syinterferopy/temporal.py b/syinterferopy/temporal.py index 9f45ada..30c61f7 100644 --- a/syinterferopy/temporal.py +++ b/syinterferopy/temporal.py @@ -11,13 +11,16 @@ #%% -def generate_random_tcs(n_tcs = 100, d_start = "20141231", d_stop = "20230801"): +def generate_random_tcs(n_tcs = 100, d_start = "20141231", d_stop = "20230801", + min_def_rate = 1., max_def_rate = 2. ): """ Generate n_tcs random time courses (ie temporal behavious of deformation. ) Inputs: n_tcs | int | Number of time series to generate. d_start | string | start date, form YYYYMMDD d_end | string | start date, form YYYYMMDD + min_def_rate | float | m/yr minimum deformation rate. + max_def_rate | float | m/yr maximum deformation rate. Returns: tcs | rank 2 array | time courses as row vectors. n_tcs rows. def_dates | list of datetimes | all the dates that we have deformation for. @@ -28,7 +31,7 @@ def generate_random_tcs(n_tcs = 100, d_start = "20141231", d_stop = "20230801"): from syinterferopy.temporal import tc_uniform_inflation for ts_n in range(n_tcs): - def_rate = np.random.uniform(0.05, 1.) # random deformation rate. + def_rate = np.random.uniform(min_def_rate, max_def_rate) # random deformation rate. tc_def, def_dates = tc_uniform_inflation(def_rate, d_start, d_stop) # generate time course if ts_n == 0: # if the first time. tcs = np.zeros((n_tcs, tc_def.shape[0])) # initialise array to store results