invert_ramp_topo_unw.py flatsim option

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=<path> [--ref_zone=<jstart,jend,istart,iend> ] [--int_path=<path>] [--out_path=<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,19 +1600,25 @@ 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)
         rms_map = np.ones((ds.RasterYSize,ds.RasterXSize))
         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()

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=<file> --int_path=<path> [--prefix=<value>] [--suffix=<value>] [--rlook=<value>]
+
+Options:
+-h --help           Show this screen.
+--int_list=<file>   Text file containing list of interferograms dates in two colums, $data1 $date2
+--int_path=<path>   Path to input interferograms directory
+--prefix=<value>    Prefix name $prefix$date1_$date2$suffix_$rlookrlks.tiff [default: '']
+--suffix=<value>    Suffix name $prefix$date1_$date2$suffix_$rlookrlks.tiff [default: '']
+--rlook=<value>     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()

NSBAS-playground/utils/nsb_plot_interferograms_network.py

@@ -22,17 +22,16 @@
 nsb_plot_interferograms_network
 
 Usage:
-  nsb_plot_interferograms_network [-o <out_path>] <pair_file> <baseline_file>
+  nsb_plot_interferograms_network [-o <out_path>] <pair_file> <baseline_file> [-w <y/n>]
   nsb_plot_interferograms_network -h | --help
 
 Options:
   -o OUT_PATH  Write figure to file.
   -h --help    Show this screen.
+  -w <y/n>     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,16 @@
 if arguments["-o"]:
    mpl.use('pdf')
 
+if arguments["-w"] == None:
+    w = 'n'
+else :
+    w = 'y'
+
 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,28 +67,54 @@
 
 # Load graph from <pair_file> and <baseline_file>
 graph = pydot.Dot("interferogram_network", graph_type="digraph")
+weight = []
 for line in open(arguments["<baseline_file>"], "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["<pair_file>"], "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.PuBu
+
 # Draw nodes
 x, y = [], []
 for n in graph.get_nodes():
     x.append(date2num(datetime.strptime(n.get_label(), "%Y%m%d")))
     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 +132,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")