compute htg from data

src/ich/InSAR_canopy_height.py

@@ -39,8 +39,8 @@ parser = argparse.ArgumentParser(description="Generate canopy height map from In
 
 parser.add_argument("-c", "--correlationFile", dest="corFile", help="Correlation file [0,1]")
 parser.add_argument("-l", "--cal_ht", dest="lidarFile", help="Calibration height file (e.g., LiDAR heights in meters)")
-parser.add_argument("-ll", "--lower_limit", dest="htl", default=0, help="Lower limit of canopy height (m)", type=int)
-parser.add_argument("-ul", "--upper_limit", dest="htg", default=40, help="Upper limit of canopy height (m)", type=int)
+parser.add_argument("-ll", "--lower_limit", dest="htl", default=None, help="Lower limit of canopy height (m)", type=int)
+parser.add_argument("-ul", "--upper_limit", dest="htg", default=None, help="Upper limit of canopy height (m)", type=int)
 parser.add_argument("-w", "--window", dest="window_size", default=10, help="Window size", type=int)
 parser.add_argument("-val", "--validation", dest="validation", default=0, help="Fraction to split for cross-validation", type=float)
 

src/ich/algo.py

@@ -114,7 +114,8 @@ def cal_(cor, gedi, htl, htg):
         tempdf = pd.DataFrame(data={'N': result[:, 0], 'S': result[:, 1], 'C': result[:, 2],
                                     'r': result[:, 3], 'rmse': result[:, 4]})
         tempdf.dropna(subset=['rmse'], inplace=True)
-        tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
+        if nn>6:
+            tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
 
         if tempdf.empty:
             sCoarse = cCoarse = 0
@@ -147,7 +148,8 @@ def cal_(cor, gedi, htl, htg):
             tempdf = pd.DataFrame(data={'N': result[:, 0], 'S': result[:, 1], 'C': result[:, 2],
                                         'r': result[:, 3], 'rmse': result[:, 4]})
             tempdf.dropna(subset=['rmse'], inplace=True)
-            tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
+            if nn>6:
+                tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
 
         if tempdf.empty:
             return [0, 0, 0, 0, 0]
@@ -170,7 +172,7 @@ def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_):
         lidarBlock = lidarArray[start_i:end_i, start_j:end_j]
         cohBlock = cohArray[start_i:end_i, start_j:end_j]
 
-        if np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 3:
+        if np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 1:
             parm = cal_(cohBlock, lidarBlock, htl, htg)
             mask = np.where(~np.isnan(lidarBlock), 1, 0)
 

src/ich/args_in.py

@@ -76,6 +76,15 @@ def rvog_inverse(args):
 
         # Apply masks to filter out invalid or extreme values
         lidar_ht[lidar_ht == 0] = np.nan
+
+        if args.htl==None and args.htg==None:
+            """ Calculating the upper limit of the heights to be estimated from the LiDAR data 
+            Ceiling the upper limit to the closest multiple of 5 for the height at 99.9 percentile """
+            args.htl=0
+            pp99 = np.nanpercentile(lidar_ht,99.9)
+            args.htg = int( 5 * np.ceil( pp99 / 5. ))
+
+
         lidar_ht[lidar_ht < args.htl] = np.nan
         lidar_ht[lidar_ht > args.htg] = np.nan
         cor[cor < 0.1] = np.nan
@@ -116,8 +125,7 @@ def rvog_inverse(args):
             os.mkdir(outDir)
 
         # Create postfix string for output files
-        _postfix_str = '_w' + str(args.window_size) + 'o' + "_" + str(args.htl) + str(args.htg)
-        # _postfix_str = '_w' + str(args.window_size) + 'o' + "%02d" % (args.window_overlap * 10) + "_" + str(args.htl) + str(args.htg)
+        _postfix_str = f"_w{args.window_size}_{args.htl}_{args.htg}"
 
         # Define output file paths
         outcFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_c.tif')
@@ -129,10 +137,15 @@ def rvog_inverse(args):
         outhtFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_ht.tif')
         outcntFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_count.tif')
         pltName = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '.png')
-        pltvalName = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_val.png')
+        
+        if args.validation > 0 and args.validation < 1:
+            pltName = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + f"_cal_{int(100-args.validation*100)}.png")
+        
+        pltvalName = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + f"_val_{int(args.validation*100)}.png")
 
         # Generate calibration parameters
-        print("[1/3] Generating calibration parameters...")
+        print("[1/3] Generating calibration parameters... ",end=" ")
+        t1 = time.time()
         s, c, rmse__, ht_, count = dynamicWindow(cor, lidar_ht_cal, args.window_size, args.htl, args.htg)
 
         # Mask out invalid values
@@ -151,8 +164,12 @@ def rvog_inverse(args):
             count = unpad(count, pad_width)
             temp_cor = unpad(cor, pad_width)
 
+        print("%0.2f sec"%(time.time()-t1))
         # Plot calibration results
         plot_data(lidar_ht_cal, ht_, args.htl, args.htg, 'Lidar height (m)', 'InSAR inverted (m)', pltName)
+        
+        t2 = time.time()
+        print("[2/3] Generating height map... ",end=" ")
 
         # Interpolate calibration parameters
         ci = spatial_intp_lin(c)
@@ -173,10 +190,10 @@ def rvog_inverse(args):
         c_temp[c_temp < 0] = np.nan
         height = []
 
-        print("[2/3] Generating height map...")
         height = arc_sinc_vectorized(gama, c_temp)
         height = np.array(height).reshape(np.shape(temp_cor))
         height[height < 0] = 0
+        print("%0.2f sec"%(time.time()-t2))
 
         # If validation is enabled, save validation results
         if args.validation > 0 and args.validation < 1:
@@ -185,8 +202,10 @@ def rvog_inverse(args):
                       lidar_ht_val, os.path.basename(args.corFile))
             plot_data(lidar_ht_val, height, args.htl, args.htg, 'Lidar height (m)', 'InSAR inverted (m)', pltvalName)
 
+        
+        t3 = time.time()
         # Write all output files
-        print("[3/3] Writing output...")
+        print("[3/3] Writing output... ",end=" ")
         write_tif(outhtFile, height, os.path.basename(args.corFile))
         write_tif(outcFile, c, os.path.basename(args.corFile))
         write_tif(outsFile, s, os.path.basename(args.corFile))
@@ -194,8 +213,9 @@ def rvog_inverse(args):
         write_tif(outciFile, ci, os.path.basename(args.corFile))
         write_tif(outsiFile, si, os.path.basename(args.corFile))
         write_tif(outcntFile, count, os.path.basename(args.corFile))
+        print("%0.2f sec"%(time.time()-t3))
 
-        print("Completed in %0.2f sec" % (time.time() - t0))
+        print("\n      Total computing time : %0.2f sec" % (time.time() - t0))
 
     except Exception as e:
         print("Task failed due to", e)

src/ich/plotting.py

@@ -105,4 +105,7 @@ def plot_data(x,y,htl,htg,xlab,ylab,pltName):
     plt.tight_layout()
     if pltName:
         plt.savefig(pltName)
-    print('%s, r = %0.2f, RMSE = %0.2f N = %d'%(pltName,r_value,RMSE,np.size(x)))
\ No newline at end of file
+    if '_val' in pltName:
+        print('%s %s; \t r = %0.2f, RMSE = %0.2f N = %d'%("      Scatterplot saved to: ",pltName,r_value,RMSE,np.size(x)))
+    else:
+        print('%s %s; \t r = %0.2f, RMSE = %0.2f N = %d'%("      Scatterplot saved to: ",pltName,r_value,RMSE,np.size(x)))
\ No newline at end of file