diff --git a/src/ich/InSAR_canopy_height.py b/src/ich/InSAR_canopy_height.py
index f6eb711c0435b46a9bf0aadb8e4f6dfee39dcf20..468d29de56055731d780e0fc3e1f7bd555fb4ade 100644
--- a/src/ich/InSAR_canopy_height.py
+++ b/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)
diff --git a/src/ich/algo.py b/src/ich/algo.py
index 72dc89e385c7ba91050f50ca34f92e8f6672c734..fc448c6110695c2374e2e2437d96fb7babe88522 100644
--- a/src/ich/algo.py
+++ b/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)
diff --git a/src/ich/args_in.py b/src/ich/args_in.py
index 4c38560d6e4e08dd6cf57aada9cf9e47e57a809f..bb11d9506d898a8bcc5568f9bd681d4d526a12db 100644
--- a/src/ich/args_in.py
+++ b/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)
diff --git a/src/ich/plotting.py b/src/ich/plotting.py
index a1f6d0270a3c159a7420bb4ead3b0da81912155b..a90c065481189899ee5324f7a369377ca3da07f0 100644
--- a/src/ich/plotting.py
+++ b/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