From d194d510a8a1931b87dfb349a5c0bf75677c33a9 Mon Sep 17 00:00:00 2001 From: hugowatine Date: Thu, 22 May 2025 15:46:31 +0200 Subject: [PATCH 1/6] invert_ramp_topo_unw.py flatsim option --- AtmoCorr/invert_ramp_topo_unw.py | 52 ++++++++++++++++++++------------ 1 file changed, 32 insertions(+), 20 deletions(-) diff --git a/AtmoCorr/invert_ramp_topo_unw.py b/AtmoCorr/invert_ramp_topo_unw.py index 18fdbbb..4883dd6 100755 --- a/AtmoCorr/invert_ramp_topo_unw.py +++ b/AtmoCorr/invert_ramp_topo_unw.py @@ -10,7 +10,7 @@ invert_ramp_topo_unw.py ------------- Estimates atmospheric phase/elevation correlations or/and azimuthal and range ramps polynomial coefficients -on unwrapped interferograms (ROI_PAC/GAMMA/GTIFF). Temporal inversion of all coeffient with a strong weight for +on unwrapped interferograms (ROI_PAC/GAMMA/GTIFF/FLATSIM). Temporal inversion of all coeffient with a strong weight for small temporal baselines and lage cover interferograms. Reconstruction of the empirical phase correction. usage: invert_ramp_topo_unw.py --int_list= [--ref_zone= ] [--int_path=] [--out_path=] \ @@ -59,7 +59,7 @@ --mask PATH Mask in .r4 format. Keep only values > threshold_mask. [default:None] --threshold_mask Thresclean_r4.pyhold on mask: take only values > threshold_mask [default: -1] --cohpixel yes/no If Yes, use amplitude interferogram to weight and mask pixels (e.g Coherence, Colinearity, Amp Filter) [default: no] ---format VALUE Format input files: ROI_PAC, GAMMA, GTIFF [default: ROI_PAC] +--format VALUE Format input files: ROI_PAC, GAMMA, GTIFF, FLATSIM [default: ROI_PAC] --threshold_coh VALUE Threshold on cohpixel file [default:0] --ibeg_mask VALUE Start line number for the mask [default: None] --iend_mask VALUE Stop line number for the mask [default: None] @@ -1286,10 +1286,13 @@ def empirical_cor(kk): lines, cols = ds.RasterYSize, ds.RasterXSize - elif sformat == 'GTIFF': - folder = 'int_'+ str(date1) + '_' + str(date2) + '/' - infile = int_path + folder + prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.tiff' - + elif sformat == 'GTIFF' or sformat == 'FLATSIM': + if sformat == 'GTIFF': + infile = int_path + prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.geo.unw.tif' + if sformat == 'FLATSIM': + folder = 'int_'+ str(date1) + '_' + str(date2) + '/' + infile = int_path + folder + prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.tiff' + checkinfile(infile) ds = gdal.Open(infile, gdal.GA_ReadOnly) @@ -1321,9 +1324,12 @@ def empirical_cor(kk): rms_map[:ds.RasterYSize,:ds.RasterXSize] = ds_band1.ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize)[:mlines,:mcols] k = np.nonzero(np.logical_or(rms_map==0.0, rms_map==9999)) rms_map[k] = float('NaN') - elif sformat == 'GTIFF': - folder = 'int_'+ str(date1) + '_' + str(date2) + '/' - rmsfile= int_path + folder + 'CNES_Coh_geo_' + str(date1) + '_' + str(date2) + rlook + '.tiff' + elif sformat == 'GTIFF' or sformat == 'FLATSIM': + if sformat == 'GTIFF': + rmsfile= int_path + str(date1) + '_' + str(date2) + '.geo.cc.tif' + if sformat == 'FLATSIM': + folder = 'int_'+ str(date1) + '_' + str(date2) + '/' + rmsfile= int_path + folder + 'CNES_Coh_geo_' + str(date1) + '_' + str(date2) + rlook + '.tiff' checkinfile(rmsfile) ds = gdal.Open(rmsfile, gdal.GA_ReadOnly) ds_band1 = ds.GetRasterBand(1) @@ -1547,7 +1553,7 @@ def empirical_cor(kk): plt.colorbar(cax, cax=c) fig.tight_layout() - if sformat == 'ROI_PAC' or sformat == 'GTIFF': + if sformat == 'ROI_PAC' or sformat == 'GTIFF' or sformat == 'FLATSIM': fig.savefig(int_path + folder + idate +'corrections.png', format='PNG') else: fig.savefig(out_path + idate +'corrections.png', format='PNG') @@ -1594,10 +1600,14 @@ def apply_cor(kk, sp, sp_inv): rms_map = ds_band1.ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize) lines, cols = ds.RasterYSize, ds.RasterXSize - elif sformat == 'GTIFF': + elif sformat == 'GTIFF' or sformat == 'FLATSIM': + if sformat == 'GTIFF': + infile = int_path + prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.geo.unw.tif' + outfile = out_path + prefix + str(date1) + '_' + str(date2) + suffix + suffout + rlook + '.geo.unw.tif' + if sformat == 'FLATSIM': + infile = int_path + folder + prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.tiff' + outfile = out_path + folder + prefix + str(date1) + '_' + str(date2) + suffix + suffout + rlook + '.tiff' - infile = int_path + prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.geo.unw.tif' - outfile = out_path + prefix + str(date1) + '_' + str(date2) + suffix + suffout + rlook + '.geo.unw.tif' ds = gdal.Open(infile, gdal.GA_ReadOnly) ds_band2 = ds.GetRasterBand(1) los_map = ds_band2.ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize) @@ -1605,8 +1615,10 @@ def apply_cor(kk, sp, sp_inv): lines, cols = ds.RasterYSize, ds.RasterXSize if rmsf == 'yes': - rmsfile= int_path + str(date1) + '_' + str(date2) + 'geo.cc.tif' - rmsfile= int_path + str(date1) + '_' + str(date2) + '.geo.cc.tif' + if sformat == 'GTIFF': + rmsfile= int_path + str(date1) + '_' + str(date2) + 'geo.cc.tif' + if sformat == 'FLATSIM': + rmsfile= int_path + folder + 'CNES_Coh_geo_' + str(date1) + '_' + str(date2) + rlook + '.tiff' ds = gdal.Open(rmsfile, gdal.GA_ReadOnly) ds_band1 = ds.GetRasterBand(1) rms_map = ds_band1.ReadAsArray(0, 0, ds.RasterXSize, ds.RasterYSize) @@ -1687,7 +1699,7 @@ def apply_cor(kk, sp, sp_inv): dst_band2.FlushCache() del dst_ds, ds - elif sformat == 'GTIFF': + elif sformat == 'GTIFF' or sformat == 'FLATSIM': dst_ds = driver.Create(outfile, cols, lines, 1, gdal.GDT_Float32) dst_band2 = dst_ds.GetRasterBand(1) dst_band2.WriteArray(flatlos,0,0) @@ -1749,7 +1761,7 @@ def apply_cor(kk, sp, sp_inv): plt.colorbar(cax, cax=c) fig.tight_layout() - if sformat == 'ROI_PAC': + if sformat == 'ROI_PAC' or sformat == 'FLATSIM': fig.savefig( int_path + folder + idate + '_reconstruc_corrections.png', format='PNG') else: fig.savefig( out_path + idate + '_reconstruc_corrections.png', format='PNG') @@ -1784,7 +1796,7 @@ def apply_cor(kk, sp, sp_inv): int_path=arguments["--int_path"] if arguments["--out_path"] == None: - out_path=np.copy(int_path) + out_path=np.copy(int_path).item() else: out_path=arguments["--out_path"] + '/' makedirs(out_path) @@ -1984,7 +1996,7 @@ def apply_cor(kk, sp, sp_inv): # load ref to define mlines, mcols and format if ref is not None: ds_extension = os.path.splitext(ref)[1] - if sformat == 'GTIFF': + if sformat == 'GTIFF' or sformat == 'FLATSIM': geotiff = arguments["--ref"] georef = gdal.Open(ref) gt = georef.GetGeoTransform() @@ -2011,7 +2023,7 @@ def apply_cor(kk, sp, sp_inv): fid.close() else: - if sformat == 'GTIFF': + if sformat == 'GTIFF' or sformat == 'FLATSIM': geotiff = arguments["--geotiff"] georef = gdal.Open(radar) gt = georef.GetGeoTransform() From 04fc05877f92a8882a4465ee0fb6de2cd80e01e0 Mon Sep 17 00:00:00 2001 From: hugowatine Date: Fri, 23 May 2025 10:41:12 +0200 Subject: [PATCH 2/6] add check_coverage_tif.py script --- NSBAS-playground/utils/check_coverage_tif.py | 102 +++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100755 NSBAS-playground/utils/check_coverage_tif.py diff --git a/NSBAS-playground/utils/check_coverage_tif.py b/NSBAS-playground/utils/check_coverage_tif.py new file mode 100755 index 0000000..f4784d3 --- /dev/null +++ b/NSBAS-playground/utils/check_coverage_tif.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +############################################ +# +# PyGdalSAR: An InSAR post-processing package +# written in Python-Gdal +# +############################################ +# Author : Simon DAOUT (Oxford) / Hugo Watine (CRPG) +############################################ + +"""\ +check_coverage_tiff.py +------------- +Check coverage unwrapped ifgs, can be used for Flatsim ifgs + +Usage: check_coverage_tiff.py --int_list= --int_path= [--prefix=] [--suffix=] [--rlook=] + +Options: +-h --help Show this screen. +--int_list= Text file containing list of interferograms dates in two colums, $data1 $date2 +--int_path= Path to input interferograms directory +--prefix= Prefix name $prefix$date1_$date2$suffix_$rlookrlks.tiff [default: ''] +--suffix= Suffix name $prefix$date1_$date2$suffix_$rlookrlks.tiff [default: ''] +--rlook= look int. $prefix$date1_$date2$suffix_$rlookrlks.tiff [default: 0] +""" + +from osgeo import gdal +import numpy as np +import docopt +import warnings, sys +from os import getcwd, path +from osgeo import gdalconst + + +warnings.filterwarnings("ignore", category=FutureWarning) +warnings.filterwarnings("ignore", category=RuntimeWarning) + +# read arguments +arguments = docopt.docopt(__doc__) +int_list=arguments["--int_list"] +if arguments["--int_path"] == None: + int_path='./' +else: + int_path=arguments["--int_path"] + '/' +if arguments["--prefix"] == None: + prefix = '' +else: + prefix=arguments["--prefix"] +if arguments["--suffix"] == None: + suffix = '' +else: + suffix=arguments["--suffix"] +if arguments["--rlook"] == None: + rlook = '' +else: + rlook = '_' + arguments["--rlook"] + 'rlks' + +def check(kk): + date1, date2 = date_1[kk], date_2[kk] + idate = str(date1) + '_' + str(date2) + folder = 'int_'+ str(date1) + '_' + str(date2) + '/' + infile=int_path + folder + prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.tiff' + name = prefix + str(date1) + '_' + str(date2) + suffix + rlook + '.tiff' + + if path.exists(infile) is not False: + print("Open: {0} in {1}".format(name,folder)) + ds = gdal.Open(infile, gdalconst.GA_ReadOnly) + los = ds.GetRasterBand(1).ReadAsArray() + del ds + else: + print("File: {0} not found in {1}".format(name,folder)) + los = np.zeros((2,2)) + los[np.isnan(los)] = 0. + return los, name + +print('----------------------------------') +print('Check interferograms list:', int_list) +date_1,date_2=np.loadtxt(int_list,comments="#",unpack=True,usecols=(0,1),dtype='i,i') +Nifg=len(date_1) +print("number of interferogram: {}".format(Nifg)) +print() + +ListInterfero = "interf_pair_coverage.txt" +print('Save empty interferogram list:', ListInterfero) +wf = open(ListInterfero, 'w') + +for j in range(Nifg): + try: + los,name = check(j) + size = np.shape(los)[0]* np.shape(los)[1] + unw = np.count_nonzero(los==0) + cov = 1 - (unw / size) + except: + print("Cannot open {0}: !!!!".format(name)) + cov = 0.0 + print("Unwrapping coverage {0}: {1} ".format(name,cov)) + wf.write("%i %i %f\n" % (date_1[j], date_2[j], cov)) +wf.close() + +print('----------------------------------') +print() From f36e1c93757ea49deb971afefcbe8e243c74fa55 Mon Sep 17 00:00:00 2001 From: hugowatine Date: Tue, 3 Jun 2025 15:36:44 +0200 Subject: [PATCH 3/6] wieght option for nsb_plot_interferograms_network.py --- .../utils/nsb_plot_interferograms_network.py | 59 +++++++++++++++---- 1 file changed, 46 insertions(+), 13 deletions(-) diff --git a/NSBAS-playground/utils/nsb_plot_interferograms_network.py b/NSBAS-playground/utils/nsb_plot_interferograms_network.py index e024a03..ec79247 100755 --- a/NSBAS-playground/utils/nsb_plot_interferograms_network.py +++ b/NSBAS-playground/utils/nsb_plot_interferograms_network.py @@ -22,17 +22,16 @@ nsb_plot_interferograms_network Usage: - nsb_plot_interferograms_network [-o ] + nsb_plot_interferograms_network [-o ] [-w ] nsb_plot_interferograms_network -h | --help Options: -o OUT_PATH Write figure to file. -h --help Show this screen. + -w show the weight used for TS analysis in the third column of pair_file (default: no) """ - - import string, os from datetime import datetime @@ -48,9 +47,17 @@ if arguments["-o"]: mpl.use('pdf') +if arguments["-w"] == None: + w = 'n' +else : + w = 'y' +print(w) + import matplotlib.pyplot as plt from matplotlib import dates from matplotlib.dates import date2num, num2date +import matplotlib.cm as cm +import matplotlib.colors as colors import pydot #except ModuleNotFoundError as e: @@ -61,13 +68,25 @@ # Load graph from and graph = pydot.Dot("interferogram_network", graph_type="digraph") +weight = [] for line in open(arguments[""], "r"): lines = line.split() date, pbaseline, rest = lines[0], lines[1], lines[2] graph.add_node(pydot.Node(date.strip(), label=date.strip(), bperp=float(pbaseline))) for line in open(arguments[""], "r"): - date1, date2 = line.split()[0], line.split()[1] - graph.add_edge(pydot.Edge(date1.strip(), date2.strip())) + if w == 'n': + date1, date2 = line.split()[0], line.split()[1] + graph.add_edge(pydot.Edge(date1.strip(), date2.strip())) + weight.append(1) + else : + date1, date2, wei = line.split()[0], line.split()[1], float(line.split()[2]) + graph.add_edge(pydot.Edge(date1.strip(), date2.strip())) + weight.append(wei) + +if w == 'y': + norm = colors.Normalize(vmin=min(weight), vmax=max(weight)) + cmap = cm.rainbow + # Draw nodes x, y = [], [] for n in graph.get_nodes(): @@ -75,14 +94,28 @@ y.append(float(n.get_attributes()["bperp"])) ax.plot(x, y, "o", color='dodgerblue', mec='black', markersize=4, picker=5) # Draw arrows -for edge in graph.get_edges(): +for i, edge in enumerate(graph.get_edges()): master = graph.get_node(edge.get_source())[0] slave = graph.get_node(edge.get_destination())[0] x = date2num(datetime.strptime(master.get_label(), "%Y%m%d")) y = float(master.get_attributes()["bperp"]) dx = date2num(datetime.strptime(slave.get_label(), "%Y%m%d")) - x dy = float(slave.get_attributes()["bperp"]) - y - ax.arrow(x, y, dx, dy, linewidth=.5, color='black', alpha=.5) + + if w == 'y': + wei = weight[i] + color = cmap(norm(wei)) + else: + color = 'black' + + ax.arrow(x, y, dx, dy, linewidth=.5, color=color, alpha=.5) + +if w == 'y': + sm = cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + cbar = plt.colorbar(sm, ax=ax) + cbar.set_label('Weight') + # Register click usage to display date of nearest point def onpick(event): global ax, an @@ -100,12 +133,12 @@ def onpick(event): an.set_text(num2date(x).strftime("%Y/%m/%d")) an.set_visible(True) fig.canvas.draw() -an = ax.annotate("", - xy=(0, 0), xycoords="data", - xytext=(0, 0), textcoords="data", - arrowprops=dict(facecolor="black", width=1, frac=0.3), - bbox=dict(boxstyle="round", fc="w")) -an.set_visible(False) +#an = ax.annotate("", +# xy=(0, 0), xycoords="data", +# xytext=(0, 0), textcoords="data", +# arrowprops=dict(facecolor="black", width=1, frac=0.3), +# bbox=dict(boxstyle="round", fc="w")) +#an.set_visible(False) fig.canvas.mpl_connect("pick_event", onpick) # Show #fig.canvas.set_window_title("Interferogram network") From c8a0b5c4982609696786e3678d2e68a17be890ae Mon Sep 17 00:00:00 2001 From: hugowatine Date: Tue, 3 Jun 2025 16:29:41 +0200 Subject: [PATCH 4/6] change colorbar --- NSBAS-playground/utils/nsb_plot_interferograms_network.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/NSBAS-playground/utils/nsb_plot_interferograms_network.py b/NSBAS-playground/utils/nsb_plot_interferograms_network.py index ec79247..0fba3fc 100755 --- a/NSBAS-playground/utils/nsb_plot_interferograms_network.py +++ b/NSBAS-playground/utils/nsb_plot_interferograms_network.py @@ -51,7 +51,6 @@ w = 'n' else : w = 'y' -print(w) import matplotlib.pyplot as plt from matplotlib import dates @@ -85,7 +84,7 @@ if w == 'y': norm = colors.Normalize(vmin=min(weight), vmax=max(weight)) - cmap = cm.rainbow + cmap = cm.PuBu # Draw nodes x, y = [], [] From 279ee1c6563c1304b194652f1e4e17815e93d933 Mon Sep 17 00:00:00 2001 From: hugowatine Date: Tue, 19 Aug 2025 15:14:06 +0200 Subject: [PATCH 5/6] Add an acceleration fonction in the inversion [default: no] --- TimeSeries/invers_disp_pixel.py | 285 +++++++++++++++++++++----------- 1 file changed, 184 insertions(+), 101 deletions(-) diff --git a/TimeSeries/invers_disp_pixel.py b/TimeSeries/invers_disp_pixel.py index 4ca1484..4dd67f7 100755 --- a/TimeSeries/invers_disp_pixel.py +++ b/TimeSeries/invers_disp_pixel.py @@ -11,7 +11,7 @@ Temporal decomposition of the time series delays of selected pixels (used depl_cumule (BIP format) and images_retenues, output of invers_pixel). Usage: invers_disp_pixel.py --cols= --lines= [--cube=] [--list_images=] [--windowsize=] [--windowrefsize=] [--lectfile=] [--aps=] \ -[--linear=] [--steps=] [--postseismic=] [--seasonal=] [--seasonal_increase=] [--vector=] [--info=]\ +[--linear=] [--acceleration=] [--steps=] [--postseismic=] [--seasonal=] [--seasonal_increase=] [--vector=] [--info=]\ [--semianual=] [--bianual=] [--degreeday=] [--bperp=] [--imref=] [--cond=] [--slowslip=] [--ineq=] \ [--name=] [--scale=] [--plot=] [] [] [--bounds=] [--dateslim=] [--plot_dateslim=] [--color=] [--fillstyle=] [--ndatasets=] @@ -24,11 +24,12 @@ --cube PATH Path to displacement file [default: depl_cumul_flat] --list_images PATH Path to list images file made of 5 columns containing for each images 1) number 2) Doppler freq (not read) 3) date in YYYYMMDD format 4) numerical date 5) perpendicular baseline [default: images_retenues] --windowsize VALUE Number of pixels around the pixel defining the window [default: 0] ---windowrefsize VALUE Number of pixels around the referenced pixel defining the window [default: windowsize] +--windowrefsize VALUE Number of pixels around the referenced pixel defining the window [default: windowsize] --lectfile PATH Path to the lect.in file (output of invers_pixel) [default: lect.in] --aps PATH Path to the APS file giving the error associated to each dates [default: No weigthing] ---linear PATH Add a linear function to the inversion ---steps PATH Add heaviside functions to the inversion .Indicate steps time (e.g 2004.,2006.) +--linear PATH Add a linear function to the inversion +--acceleration=PATH Add an acceleration fonction in the inversion [default: no] +--steps PATH Add heaviside functions to the inversion .Indicate steps time (e.g 2004.,2006.) --postseismic PATH Add logarithmic transients to each steps step. Indicate characteristic time of the log function, must be a serie of values of the same lenght than steps (e.g 1.,1.). To not associate postseismic function to a given steps step, put None (e.g None,1.) --slowslip VALUE Add slow-slip function in the inversion (as defined by Larson et al., 2004). Indicate median and characteristic time of the events (e.g. 2004.,1,2006,0.5) [default: None] --seasonal PATH If yes, add seasonal terms in the inversion @@ -38,15 +39,15 @@ --bianual PATH If yes, add bianual terms in the inversion --vector PATH Path to the vector text files containing a value for each dates [default: None] --info PATH Path to extra file in r4 or tif format to plot is value on the selected pixel, e.g. aspect [default: None]. ---bperp PATH If yes, add term proportional to the perpendicular baseline in the inversion +--bperp PATH If yes, add term proportional to the perpendicular baseline in the inversion --imref VALUE Reference image number [default: 1] ---cond VALUE Condition value for optimization: Singular value smaller than cond are considered zero [default: 1.e-3] ---ineq VALUE If yes, add ineguality constrained in the inversion: use least square result to iterate the inversion. Force postseismic to be the same sign than steps [default: no]. +--cond VALUE Condition value for optimization: Singular value smaller than cond are considered zero [default: 1.e-6] +--ineq VALUE If yes, add inequality constrained in the inversion: use least square result to iterate the inversion. Force postseismic to be the same sign than steps [default: no]. --name Value Name output figures [default: None] ---scale Scaling value between input data and desired output [default: 1] +--scale Scaling value between input data and desired output [default: 1] --plot Display results [default: yes] -iref colum numbers of the reference pixel [default: None] -jref lign number of the reference pixel [default: None] +iref colum numbers of the reference pixel [default: None] +jref lign number of the reference pixel [default: None] --bounds yMin,yMax time series plots --dateslim Datemin,Datemax time series --plot_dateslim Datemin,Datemax time series for plot only @@ -98,6 +99,50 @@ # docopt (command line parser) import docopt +def write_envi_hdr(filename, shape, dtype='float32', interleave='bip'): + """ + Crée un fichier ENVI .hdr à partir d'un fichier. + + Parameters: + - filename: nom du fichier sans extension (.hdr sera ajouté) + - shape: tuple (lines, samples, bands) + - dtype: type numpy ('float32', 'int16', ...) + - interleave: 'bsq', 'bil' ou 'bip' + """ + + if len(shape) == 3: + lines, samples, bands = shape + else: + lines, samples = shape + bands = 1 + dtype_map = { + 'uint8': 1, + 'int16': 2, + 'int32': 3, + 'float32': 4, + 'float64': 5, + 'complex64': 6, + 'complex128': 9, + 'uint16': 12, + 'uint32': 13, + 'int64': 14, + 'uint64': 15, + } + + if dtype not in dtype_map: + raise ValueError(f"Unsupported data type for ENVI: {dtype}") + + hdr_content = f"""ENVI +samples = {samples} +lines = {lines} +bands = {bands} +data type = {dtype_map[dtype]} +interleave = {interleave} +byte order = 0 +""" + + with open(filename + '.hdr', 'w') as f: + f.write(hdr_content) ######################################################################## # Define basis functions @@ -152,6 +197,15 @@ def __init__(self,name,reduction,date): def g(self,t): func=(t-self.to) return func + +class acceleration(pattern): + def __init__(self,name,reduction,date): + pattern.__init__(self,name,reduction,date) + self.to=date + + def g(self,t): + func= (t-self.to)**2 + return func class sint(pattern): def __init__(self,name,reduction,date): @@ -321,13 +375,11 @@ def date2dec(dates): arguments = docopt.docopt(__doc__) if arguments["--list_images"] == None: - listim = "images_retenues" -else: - listim = arguments["--list_images"] + arguments["--list_images"] = "images_retenues" if arguments["--ineq"] == None: arguments["--ineq"] = 'no' if arguments["--cond"] == None: - rcond = 1e-3 + rcond = 1e-6 else: rcond = float(arguments["--cond"]) if arguments["--imref"] == None: @@ -358,6 +410,10 @@ def date2dec(dates): inter = 'no' else: inter = arguments["--linear"] +if arguments["--acceleration"] == None: + acc = 'no' +else: + acc = arguments["--acceleration"] if arguments["--seasonal"] == None: seasonal = 'no' else: @@ -483,13 +539,14 @@ def date2dec(dates): # Charger les colonnes pour idates et dates (colonnes 1 et 3) # Générer dynamiquement les colonnes pour chaque base (à partir de la colonne 5) columns = [1, 3] + list(range(5, 5 + ndata)) - + # Charger les données en utilisant `usecols=columns` data = np.loadtxt(arguments["--list_images"], comments='#', usecols=columns, unpack=True) # Extraire les données de dates idates = data[0].astype(int) dates = data[1] + # Extraire les bases et initialiser les références for i in range(ndata): @@ -655,6 +712,11 @@ def date2dec(dates): indexinter=index index = index + 1 +if acc =='yes': + basis.append(acceleration(name='acceleration', reduction='acc', date=datemin)) + indexacc=index + index += 1 + if degreeday=='yes': indexdd = index basis.append(stefan(name='degree-day',reduction='ddt',date=datemin,tcar1=ddt,tcar2=ddf)) @@ -739,9 +801,9 @@ def date2dec(dates): indexco = indexco.astype(int) indexsse = indexsse.astype(int) -eguality = False +equality = False if arguments["--seasonal_increase"] == 'yes' and arguments["--seasonal"] == 'yes': - eguality = True + equality = True arguments["--ineq"] = 'yes' # define size G matrix @@ -754,86 +816,88 @@ def date2dec(dates): for i in range((Mbasis)): basis[i].info() -## inversion procedure -def consInvert(A,b,sigmad,ineq='yes',cond=1.0e-3, iter=100,acc=1e-6, eguality=False): - '''Solves the constrained inversion problem. - - Minimize: - - ||Ax-b||^2 - - Subject to: - mmin < m < mmax - ''' +def consInvert(A, b, sigmad, ineq='yes', cond=1e-6, iter=60, acc=5e-4, equality=False): + """ + Résout Ax ≈ b sous contraintes d'inégalité (et éventuellement d’égalité). + Retourne : + fsoln : vecteur de solution + sigmam : incertitudes (diag de la covariance) + """ + global indexpo, indexco, indexpofull, pos, indexseas, indexseast if A.shape[0] != len(b): - raise ValueError('Incompatible dimensions for A and b') + raise ValueError('Dimensions incompatibles pour A et b') if ineq == 'no': - print('ineq=no: Least-squared inversion') - try: - Cd = np.diag(sigmad**2, k = 0) - fsoln = np.dot(np.linalg.inv(np.dot(np.dot(A.T,np.linalg.inv(Cd)),A)),np.dot(np.dot(A.T,np.linalg.inv(Cd)),b)) - print('LSQT solution:', fsoln) - except: - fsoln = lst.lstsq(A,b,rcond=None)[0] - + W = np.diag(1.0 / sigmad) + fsoln = np.linalg.lstsq(W @ A, W @ b, rcond=cond)[0] else: - print('ineq=yes: Iterative least-square inversion.') - if len(indexpo>0): - # invert first without post-seismic - Ain = np.delete(A,indexpo,1) - mtemp = lst.lstsq(Ain,b,rcond=cond)[0] - print('Prior obtained with SVD decomposition neglecting small eigenvectors inferior to {} (cond)'.format(cond)) - print('SVD solution:', mtemp) - - # rebuild full vector - for z in range(len(indexpo)): - mtemp = np.insert(mtemp,indexpo[z],0) - minit = np.copy(mtemp) - # # initialize bounds - mmin,mmax = -np.ones(len(minit))*np.inf, np.ones(len(minit))*np.inf - - # We here define bounds for postseismic to be the same sign than steps - # and steps inferior or egual to the coseimic initial - print('Impose postseismic to be the same sign than steps') - for i in range(len(indexco)): - if (pos[i] > 0.) and (minit[int(indexco[i])]>0.): - mmin[int(indexpofull[i])], mmax[int(indexpofull[i])] = 0, np.inf - mmin[int(indexco[i])], mmax[int(indexco[i])] = 0, minit[int(indexco[i])] - if (pos[i] > 0.) and (minit[int(indexco[i])]<0.): - mmin[int(indexpofull[i])], mmax[int(indexpofull[i])] = -np.inf , 0 - mmin[int(indexco[i])], mmax[int(indexco[i])] = minit[int(indexco[i])], 0 - + # Initialisation + if indexpo is not None and len(indexpo) > 0: + Ain = np.delete(A, indexpo, axis=1) + Win = np.diag(1.0 / np.delete(sigmad, indexpo)) + mtemp = np.linalg.lstsq(Win @ Ain, Win @ b, rcond=cond)[0] + + # Réinsérer les post-sismiques + for z in range(len(indexpo)): + mtemp = np.insert(mtemp, indexpo[z], 0.0) + minit = np.copy(mtemp) else: - minit = lst.lstsq(A,b,rcond=None)[0] - mmin,mmax = -np.ones(len(minit))*np.inf, np.ones(len(minit))*np.inf - - bounds=list(zip(mmin,mmax)) - def eq_cond(x, *args): - return (x[indexseast+1]/x[indexseast]) - (x[indexseas+1]/x[indexseas]) - - ####Objective function and derivative - _func = lambda x: np.sum(((np.dot(A,x)-b)/sigmad)**2) - _fprime = lambda x: 2*np.dot(A.T/sigmad, (np.dot(A,x)-b)/sigmad) - if eguality: - res = opt.fmin_slsqp(_func,minit,bounds=bounds,fprime=_fprime,eqcons=[eq_cond], \ - iter=iter,full_output=True,iprint=0,acc=acc) + W = np.diag(1.0 / sigmad) + minit = np.linalg.lstsq(W @ A, W @ b, rcond=cond)[0] + + # Définir les bornes + n = len(minit) + mmin = -np.ones(n) * np.inf + mmax = np.ones(n) * np.inf + + if indexpo is not None and indexco is not None: + for i in range(len(indexco)): + ico = int(indexco[i]) + ipo = int(indexpofull[i]) + if pos[i] > 0. and minit[ico] > 0.: + mmin[ipo], mmax[ipo] = 0, np.inf + mmin[ico], mmax[ico] = 0, minit[ico] + elif pos[i] > 0. and minit[ico] < 0.: + mmin[ipo], mmax[ipo] = -np.inf, 0 + mmin[ico], mmax[ico] = minit[ico], 0 + + bounds = list(zip(mmin, mmax)) + + # Fonction à minimiser + def _func(x): + return np.sum(((A @ x - b) / sigmad) ** 2) + + def _fprime(x): + return 2 * A.T @ ((A @ x - b) / sigmad**2) + + # Contraintes d'égalité + if equality: + def eq_cond(x): + return (x[indexseast + 1] / x[indexseast]) - (x[indexseas + 1] / x[indexseas]) + res = opt.fmin_slsqp(_func, minit, bounds=bounds, fprime=_fprime, + eqcons=[eq_cond], iter=iter, acc=acc, + full_output=True, iprint=0) else: - res = opt.fmin_slsqp(_func,minit,bounds=bounds,fprime=_fprime, \ - iter=iter,full_output=True,iprint=0,acc=acc) - fsoln = res[0] - print('Optimization:', fsoln) - + res = opt.fmin_slsqp(_func, minit, bounds=bounds, fprime=_fprime, + iter=iter, acc=acc, full_output=True, iprint=0) + + fsoln, fx, its, imode, msg = res + if imode != 0: + logger.warning("SLSQP did not converge: {msg}") + fsoln = minit # ou np.full_like(minit, np.nan) + + # Calcul de l'incertitude try: - varx = np.linalg.inv(np.dot(A.T,A)) - res2 = np.sum(pow((b-np.dot(A,fsoln)),2)) - scale = 1./(A.shape[0]-A.shape[1]) - sigmam = np.sqrt(scale*res2*np.diag(varx)) - except: - sigmam = np.ones((A.shape[1]))*float('NaN') - print('model errors:', sigmam) - return fsoln,sigmam + varx = np.linalg.pinv(A.T @ A) + res2 = np.sum((b - A @ fsoln) ** 2) + scale = 1. / (A.shape[0] - A.shape[1]) + sigmam = np.sqrt(scale * res2 * np.diag(varx)) + except np.linalg.LinAlgError: + sigmam = np.full(A.shape[1], np.nan) + + return fsoln, sigmam + # plot diplacements maps if Npix > 2: @@ -925,15 +989,16 @@ def eq_cond(x, *args): names.append(kernels[l].reduction) print('basis functions:', names) - print(len(k), N) - print(G[:6,:6]) - print(taby) - print(sigmad[k]) - m,sigmam = consInvert(G,taby,sigmad[k],cond=rcond, ineq=arguments["--ineq"], eguality=eguality) + #print(len(k), N) + #print(G[:6,:6]) + #print(taby) + #print(sigmad[k]) + m,sigmam = consInvert(G,taby,sigmad[k],cond=rcond, ineq=arguments["--ineq"], equality=equality) #sys.exit() # forward model in original order mdisp[k] = np.dot(G,m) + if bperp=='yes': bperperr[k] = np.dot(G[:,indexbperp[0]:indexbperp[0]+ndata],m[indexbperp[0]:indexbperp[0]+ndata]) else: @@ -959,8 +1024,10 @@ def eq_cond(x, *args): disp_seas = np.zeros((N)) disp_seast = np.zeros((N)) - if inter=='yes': + if inter=='yes' and acc=='no': lin[k] = np.dot(G[:,indexinter],m[indexinter]) + if inter=='yes' and acc=='yes': + lin[k] = np.dot(G[:,indexinter],m[indexinter]) + np.dot(G[:, indexacc], m[indexacc]) # plot data and model minus bperp error and seasonal terms if seasonal=='yes': @@ -989,7 +1056,7 @@ def eq_cond(x, *args): # convert between 0 and 2pi if phit<0: phit = phit + 2*np.pi - + print(phi,phit) print(m[indexseas+1]/m[indexseas],m[indexseast+1]/m[indexseast] ) @@ -1036,7 +1103,10 @@ def eq_cond(x, *args): ) # plot data and model minus bperp error and linear term - if inter=='yes': + if inter=='yes' and acc=='no': + ax3.plot(x,disp-bperperr-lin,markers[jj],color=color,fillstyle=fillstyle,markersize=4,label='detrended data') + ax3.errorbar(x,disp-bperperr-lin,yerr = sigmad, ecolor=color,fmt='none', alpha=0.5) + elif inter=='yes' and acc=='yes': ax3.plot(x,disp-bperperr-lin,markers[jj],color=color,fillstyle=fillstyle,markersize=4,label='detrended data') ax3.errorbar(x,disp-bperperr-lin,yerr = sigmad, ecolor=color,fmt='none', alpha=0.5) @@ -1046,6 +1116,7 @@ def eq_cond(x, *args): ax2.errorbar(x,disp-disp_seas-bperperr,yerr = sigmad, ecolor=color,fmt='none', alpha=0.3) # create synthetic time + #t = np.arange(xmin, xmax, 0.1) t = np.array([xmin + datetime.timedelta(days=d) for d in range(0, 2920)]) tdec = np.array([float(date.strftime('%Y')) + float(date.strftime('%j'))/365.1 for date in t]) mseas = np.zeros(len(tdec)) @@ -1058,8 +1129,10 @@ def eq_cond(x, *args): G[:,Mbasis+l]=np.interp(tdec,tabx,kernels[l].g(k)) model = np.dot(G,m) - if inter=='yes': + if inter=='yes' and acc=='no': model_lin = np.dot(G[:,indexinter],m[indexinter]) + elif inter=='yes' and acc=='yes': + model_lin = np.dot(G[:,indexinter],m[indexinter]) + np.dot(G[:, indexacc], m[indexacc]) # we need to substrat bperp model to the whole model if bperp=='yes': @@ -1093,15 +1166,24 @@ def eq_cond(x, *args): # plot model if inter=='yes': - if seasonal=='yes' and seasonalt=='no': + if seasonal=='yes' and seasonalt=='no' and acc=='no': ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexinter],amp,phi)) ax3.plot(t,model-model_lin-model_bperp,'-r') - elif seasonal=='yes' and seasonalt=='yes': + elif seasonal=='yes' and seasonalt=='no' and acc=='yes': + ax.plot(t,model-model_bperp,'-r',label='Acc: {:.2f}, Rate: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexacc], m[indexinter],amp,phi)) + ax3.plot(t,model-model_lin-model_bperp,'-r') + elif seasonal=='yes' and seasonalt=='yes' and acc=='no': ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexinter],ampt,phit,amp,phi)) ax3.plot(t,model-model_lin-model_bperp,'-r') - elif seasonal=='no' and seasonalt=='yes': + elif seasonal=='yes' and seasonalt=='yes' and acc=='yes': + ax.plot(t,model-model_bperp,'-r',label='Acc: {:.2f}, Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}, Amp: {:.2f}, Phi: {:.2f}'.format(m[indexacc], m[indexinter],ampt,phit,amp,phi)) + ax3.plot(t,model-model_lin-model_bperp,'-r') + elif seasonal=='no' and seasonalt=='yes' and acc=='no': ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}'.format(m[indexinter],ampt,phit)) ax3.plot(t,model-model_lin-model_bperp,'-r') + elif seasonal=='no' and seasonalt=='yes' and acc=='yes': + ax.plot(t,model-model_bperp,'-r',label='Acc: {:.2f}, Rate: {:.2f}, Ampt: {:.2f}, Phit: {:.2f}'.format(m[indexacc], m[indexinter],ampt,phit)) + ax3.plot(t,model-model_lin-model_bperp,'-r') else: ax.plot(t,model-model_bperp,'-r',label='Rate: {:.2f}'.format(m[indexinter])) ax3.plot(t,model-model_lin-model_bperp,'-r') @@ -1149,3 +1231,4 @@ def eq_cond(x, *args): if plot == 'yes': plt.show() sys.exit() + From 38b68773db5490fc9b215b7308c8d535f5e15d62 Mon Sep 17 00:00:00 2001 From: hugowatine Date: Tue, 19 Aug 2025 15:21:07 +0200 Subject: [PATCH 6/6] remove write_envi_hdr fct --- TimeSeries/invers_disp_pixel.py | 45 --------------------------------- 1 file changed, 45 deletions(-) diff --git a/TimeSeries/invers_disp_pixel.py b/TimeSeries/invers_disp_pixel.py index 4dd67f7..e9ffe01 100755 --- a/TimeSeries/invers_disp_pixel.py +++ b/TimeSeries/invers_disp_pixel.py @@ -99,51 +99,6 @@ # docopt (command line parser) import docopt -def write_envi_hdr(filename, shape, dtype='float32', interleave='bip'): - """ - Crée un fichier ENVI .hdr à partir d'un fichier. - - Parameters: - - filename: nom du fichier sans extension (.hdr sera ajouté) - - shape: tuple (lines, samples, bands) - - dtype: type numpy ('float32', 'int16', ...) - - interleave: 'bsq', 'bil' ou 'bip' - """ - - if len(shape) == 3: - lines, samples, bands = shape - else: - lines, samples = shape - bands = 1 - dtype_map = { - 'uint8': 1, - 'int16': 2, - 'int32': 3, - 'float32': 4, - 'float64': 5, - 'complex64': 6, - 'complex128': 9, - 'uint16': 12, - 'uint32': 13, - 'int64': 14, - 'uint64': 15, - } - - if dtype not in dtype_map: - raise ValueError(f"Unsupported data type for ENVI: {dtype}") - - hdr_content = f"""ENVI -samples = {samples} -lines = {lines} -bands = {bands} -data type = {dtype_map[dtype]} -interleave = {interleave} -byte order = 0 -""" - - with open(filename + '.hdr', 'w') as f: - f.write(hdr_content) - ######################################################################## # Define basis functions ########################################################################