dev

LiCSBAS_lib/LiCSBAS_inv_lib.py

@@ -8,6 +8,8 @@
 =========
 Changelog
 =========
+20240930 ML
+ - (finally) found the bug causing nans in inversion of some datasets. Fixed by removing the scipy.sparse functionality. Perhaps just csr_array would do or other tweaking?
 20240423 ML
  - parallelised singular (with correct vel/vconst estimates)
 20231101 Yasser Maghsoudi (and ML), Uni Leeds
@@ -42,15 +44,18 @@
 from astropy.stats import bootstrap
 from astropy.utils import NumpyRNGContext
 #from scipy.sparse.linalg import lsqr as sparselsq
-from scipy.sparse.linalg import lsmr as sparselsq
-from scipy.sparse import csr_array, csc_array  # csr_matrix, csc_matrix # but maybe coo_matrix would be better for G? to be checked
+#from scipy.sparse.linalg import lsmr as sparselsq
+#from scipy.sparse import csr_array, csc_array  # csr_matrix, csc_matrix # but maybe coo_matrix would be better for G? to be checked
 import LiCSBAS_tools_lib as tools_lib
 try:
     from sklearn.linear_model import RANSACRegressor
 except:
     print('not loading RANSAC (optional experimental function)')
 
 
+#debugmode = True
+#print('inversion runs in debug mode - please inform Milan if this works now')
+
 #%%
 def make_sb_matrix(ifgdates):
     """
@@ -197,7 +202,10 @@ def invert_nsbas(unw, G, dt_cum, gamma, n_core, gpu, singular=False, only_sb=Fal
             q = multi.get_context('fork')
             p = q.Pool(n_core)
             if not singular:
-                A = csc_array(Gall)  # or csr?
+                #if debugmode:
+                #A = Gall
+                #else:
+                #    A = csc_array(Gall)  # or csr?
                 _result = p.map(censored_lstsq_slow_para_wrapper, args) #list[n_pt][length]
             else:
                 from functools import partial
@@ -302,13 +310,22 @@ def censored_lstsq_slow_para_wrapper(i):
     if np.mod(i, 100) == 0:
         print('  Running {0:6}/{1:6}th point...'.format(i, unw_tmp.shape[1]), flush=True)
     m = mask[:,i] # drop rows where mask is zero
+    X = np.linalg.lstsq(Gall[m], unw_tmp[m, i], rcond=None)[0]
+    return X
+
+'''
+2024-09-30 - a BUG! the sparselsq ended up on values with e-5 where lstsq return correct values(!)
+reverting back to numpy solution
+
     try:
         #X = np.linalg.lstsq(Gall[m], unw_tmp[m,i], rcond=None)[0]
         X = sparselsq(A[m], unw_tmp[m, i], atol=1e-05, btol=1e-05)[0]
     except:
         X = np.zeros((Gall.shape[1]), dtype=np.float32)*np.nan
+        print('Warning: error during sparselsq - setting to nan')
+        print('')
     return X
-
+'''
 
 #%%
 def invert_nsbas_wls(unw, var, G, dt_cum, gamma, n_core):
@@ -734,18 +751,29 @@ def censored_lstsq_slow(A, B, M):
     X = np.empty((A.shape[1], B.shape[1]))
     # 20231101 update - not tested
     #A = csr_matrix(A) # or csc?
-    A = csc_array(A) # or csr?
+    #if not debugmode:
+    #    A = csc_array(A) # or csr?
+    errs = 0
     for i in range(B.shape[1]):
         if np.mod(i, 100) == 0:
              print('\r  Running {0:6}/{1:6}th point...'.format(i, B.shape[1]), end='', flush=True)
 
         m = M[:,i] # drop rows where mask is zero
+        X[:, i] = np.linalg.lstsq(A[m], B[m, i], rcond=None)[0]
+    return X
+
+'''
+# 20240930 - removing the sparselsq - it just does NOT do any good job! 
         try:
             #X[:,i] = np.linalg.lstsq(A[m], B[m,i], rcond=None)[0]
             # 20231101 update
             X[:, i] = sparselsq(A[m], B[m, i], atol=1e-05, btol=1e-05)[0]
         except:
+            errs = errs+1
             X[:,i] = np.nan
-
     print('')
+    if errs > 0:
+        print('Warning: '+str(errs)+' errors occurred during sparselsq (-> nan increments)')
+        print('')
     return X
+'''

LiCSBAS_lib/LiCSBAS_io_lib.py

@@ -227,6 +227,20 @@ def read_ifg_list(ifg_listfile):
     return ifgdates
 
 
+def read_epochlist(txtfile):
+    f = open(txtfile)
+    line = f.readline()
+    out = []
+    while line:
+        if (line[0] == "2") or (line[0] == "1"):
+            out.append(str(line))
+            line = f.readline()
+        else:
+            line = f.readline()
+            continue
+    return out
+
+
 #%%
 def get_param_par(mlipar, field):
     """

batch_LiCSBAS.sh

@@ -7,12 +7,16 @@
 #                                         04: LiCSBAS04op_mask_unw.py (optional)
 #                                         05: LiCSBAS05op_clip_unw.py (optional)
 #  11: LiCSBAS11_check_unw.py
+#  (optional) 120: LiCSBAS120_choose_reference.py   - RECOMMENDED, especially if nullification is used, thus added to cometdev
 #  12: LiCSBAS12_loop_closure.py
 #  13: LiCSBAS13_sb_inv.py
 #  14: LiCSBAS14_vel_std.py
 #  15: LiCSBAS15_mask_ts.py
 #  16: LiCSBAS16_filt_ts.py
 
+#
+# Status of COMET dev version - the experimental functions are turned on with 
+#
 #################
 ### Settings ####
 #################
@@ -64,6 +68,7 @@ p01_get_mli="n" # y/n
 p11_unw_thre=""	# default: 0.3
 p11_coh_thre=""	# default: 0.05
 p11_s_param="n" # y/n
+p120_use="n"  # y/n
 p12_loop_thre=""	# default: 1.5 rad. With --nullify, recommended higher value (as this is an average over the whole scene)
 p12_multi_prime="y"	# y/n. y recommended
 p12_nullify="" # y/n. y recommended
@@ -114,6 +119,7 @@ p05_outGEOCmldir_suffix="" # default: clip
 p05_n_para=$n_para   # default: # of usable CPU
 p11_GEOCmldir=""	# default: $GEOCmldir
 p11_TSdir=""	# default: TS_$GEOCmldir
+p120_ignoreconncomp="n" # y/n
 p12_GEOCmldir=""        # default: $GEOCmldir
 p12_TSdir=""    # default: TS_$GEOCmldir
 p12_n_para=$n_para	# default: # of usable CPU
@@ -137,6 +143,13 @@ p16_nomask="n"	# y/n. default: n
 p16_n_para=$n_para   # default: # of usable CPU
 
 
+# cometdev
+if [ $cometdev -gt 0 ]; then
+    # using --singular, so setting to simple LS instead of WLS
+    p13_inv_alg="LS"
+fi
+
+
 #############################
 ### Run (No need to edit) ###
 #############################
@@ -317,6 +330,22 @@ if [ $start_step -le 11 -a $end_step -ge 11 ];then
 fi
 
 if [ $start_step -le 12 -a $end_step -ge 12 ];then
+  if [ $cometdev -eq 1 ]; then
+    p120_use='y'
+  fi
+  if [ $p120_use == "y" ]; then
+    dirset="-c $GEOCmldir -d $GEOCmldir -t $TSdir "
+    extra=""
+    if [ $p120_ignoreconncomp == "y" ]; then
+        extra="--ignore_comp"
+    fi
+    if [ $check_only == "y" ];then
+      echo "LiCSBAS120_choose_reference.py $dirset "$extra
+    else
+      LiCSBAS120_choose_reference.py $dirset $extra 2>&1 | tee -a $log
+      if [ ${PIPESTATUS[0]} -ne 0 ];then exit 1; fi
+    fi
+  fi
   p12_op=""
   if [ ! -z $p12_GEOCmldir ];then p12_op="$p12_op -d $p12_GEOCmldir"; 
     else p12_op="$p12_op -d $GEOCmldir"; fi
@@ -328,14 +357,14 @@ if [ $start_step -le 12 -a $end_step -ge 12 ];then
   if [ ! -z $p12_rm_ifg_list ];then p12_op="$p12_op --rm_ifg_list $p12_rm_ifg_list"; fi
   if [ ! -z $p12_n_para ];then p12_op="$p12_op --n_para $p12_n_para";
   elif [ ! -z $n_para ];then p12_op="$p12_op --n_para $n_para";fi
-  if [ $check_only == "y" ];then
-    echo "LiCSBAS12_loop_closure.py $p12_op"
-  else
-    if [ $cometdev -eq 1 ]; then
+  if [ $cometdev -eq 1 ]; then
      extra='--nullify'
     else
      extra=''
-    fi
+  fi
+  if [ $check_only == "y" ];then
+    echo "LiCSBAS12_loop_closure.py $p12_op "$extra
+  else
     LiCSBAS12_loop_closure.py $extra $p12_op 2>&1 | tee -a $log
     if [ ${PIPESTATUS[0]} -ne 0 ];then exit 1; fi
   fi
@@ -367,7 +396,7 @@ if [ $start_step -le 13 -a $end_step -ge 13 ];then
      if [ -z $p13_n_unw_r_thre ];then
        extra=$extra' --n_unw_r_thre 0.4'
      fi
-     #extra='--singular --nopngs'
+     extra=$extra' --singular' # --nopngs'
     else
      extra=''
     fi

bin/LiCSBAS11_check_unw.py

@@ -360,9 +360,10 @@ def main(argv=None):
             rasorg = os.path.join(ifgdir, ifgdates[i], rasname)
 
             if not os.path.exists(rasorg):
-                print('\nERROR: No browse image {} available!\n'
-                      .format(rasorg), file=sys.stderr)
-                return 2
+                print('WARNING: No browse image {} available!\n'.format(rasorg))
+                print('assuming there is an error and skipping this ifg')
+                bad_ifgdates.append(ifgdates[i])
+                continue
 
         ### Identify bad ifgs and link ras
         if rate_cov[i] < unw_cov_thre or coh_avg_ifg[i] < coh_thre or \

bin/LiCSBAS120_choose_reference.py

@@ -38,6 +38,13 @@
 =====
 LiCSBAS120_choose_reference.py [-h] [-f FRAME_DIR] [-g UNW_DIR] [-t TS_DIR] [-w WIN] [-r [0-1]] [--w_unw [0-1]] [--w_coh [0-1]] [--w_con [0-1]] [--w_hgt [0-1]] [--refx [0-1]] [--refy [0-1]]
 """
+
+#%% Change log
+'''
+20240928 ML
+ - check for existence of conn. components, avoiding if not 
+'''
+
 from LiCSBAS_meta import *
 import numpy as np
 import matplotlib.pyplot as plt
@@ -330,8 +337,10 @@ def discard_ifg_with_all_nans_at_ref():
     ### Check ref exist in unw. If not, list as noref_ifg
     noref_ifg = []
     for ifgd in ifgdates:
-
-        unwfile = os.path.join(ifgdir, ifgd, ifgd+'.unw')
+        # add about ori:
+        unwfile = os.path.join(ifgdir, ifgd, ifgd+'.unw.ori')
+        if not os.path.exists(unwfile):
+            unwfile = os.path.join(ifgdir, ifgd, ifgd+'.unw')
         unw_data = io_lib.read_img(unwfile, length, width)
         unw_ref = unw_data[refy1:refy2, refx1:refx2]
 
@@ -421,6 +430,17 @@ def plot_networks():
     plot_lib.plot_strong_weak_cuts_network(retained_ifgs, bperp, weak_links, edge_cuts, node_cuts, pngfile, plot_weak=True)
 
 
+def check_components_existence():
+    missingcomp = False
+    for ifgd in ifgdates:
+        confile = os.path.join(ccdir, ifgd, ifgd + '.conncomp')
+        if not os.path.exists(confile):
+            missingcomp = True
+            break
+    if missingcomp:
+        print('At least one ifg has missing conncomp file - not using conncomps')
+        args.ignore_comp=True
+
 
 def main():
     global retained_ifgs
@@ -430,6 +450,8 @@ def main():
     read_length_width()
     decide_reference_window_size()
     get_ifgdates()
+    if not args.ignore_comp:
+        check_components_existence()
 
     calc_block_sum_of_unw_coh_component_size()
     calc_height_std()