added inversion rate threshold, and median filter

src/ich/algo.py

@@ -204,17 +204,19 @@ def cal_(temp_cor, temp_gedi, htl, htg):
                         slope, intercept, r_value, p_value, std_err = np.nan, np.nan, np.nan, np.nan, np.nan
                         RMSE = np.nan
                     
-                    result.append([len(x), S_param, C_param, r_value, RMSE])
+                    result.append([len(x), S_param, C_param, r_value, RMSE, np.round((inv_N/nn),2)])
         
         result = np.array(result)
         tempdf = pd.DataFrame(data={'N': result[:, 0], 'S': result[:, 1], 'C': result[:, 2],
-                                    'r': result[:, 3], 'rmse': result[:, 4]})
+                                    'r': result[:, 3], 'rmse': result[:, 4],'inv': result[:, 5]})
         del result
         tempdf.dropna(subset=['rmse'], inplace=True)
         if tempdf.empty:
-            return [0, 0, 0, 0, 0]        
+            return [0, 0, 0, 0, 0, 0]        
         elif nn>6:
             tempdf_coarse = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
+            tempdf_coarse = tempdf_coarse[tempdf_coarse['inv'] > 0.5]
+
             sCoarse = np.round(tempdf_coarse.iloc[0]['S'], 2)
             cCoarse = np.round(tempdf_coarse.iloc[0]['C'], 2)
 
@@ -243,14 +245,14 @@ def cal_(temp_cor, temp_gedi, htl, htg):
             
             result = np.array(result)
             tempdf = pd.DataFrame(data={'N': result[:, 0], 'S': result[:, 1], 'C': result[:, 2],
-                                        'r': result[:, 3], 'rmse': result[:, 4]})
+                                        'r': result[:, 3], 'rmse': result[:, 4], 'inv': result[:, 4]})
             del result
             tempdf.dropna(subset=['rmse'], inplace=True)
             # if nn>6:
             tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
 
             if tempdf.empty and sCoarse==0:
-                return [0, 0, 0, 0, 0]
+                return [0, 0, 0, 0, 0,0]
             
             elif tempdf.empty and sCoarse!=0:
                 list(tempdf_coarse.iloc[0])
@@ -262,7 +264,7 @@ def cal_(temp_cor, temp_gedi, htl, htg):
 
     except Exception as e:
         # print(f"Exception in cal_: {e}")
-        return [0, 0, 0, 0, 0]
+        return [0, 0, 0, 0, 0,0]
 
 def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_):
     rows, cols = cohArray.shape

src/ich/args_in.py

@@ -23,6 +23,10 @@ from scipy.interpolate import interpn
 from empatches import EMPatches
 from sklearn.model_selection import train_test_split
 import concurrent.futures
+
+from scipy.ndimage import median_filter
+
+
 emp = EMPatches()
 
 gdal.UseExceptions()
@@ -145,7 +149,7 @@ def rvog_inverse(args):
             count = np.zeros(np.shape(temp_cor))
 
 
-            parm_ = [0,0,0,0,0]
+            # parm_ = [0,0,0,0,0]
             # print("[1/3] Generating calibration parameters...")
 
             batch_size = 100  # Define your batch size
@@ -167,6 +171,24 @@ def rvog_inverse(args):
             # Unpack results
             s, c, rmse__, count, ht_ = zip(*results)
 
+            temp_mask = np.zeros(temp_lidar.shape)
+            for win in tqdm(range(np.shape(temp_lidar)[0])):
+                mask = temp_lidar[win,:,:].copy()
+                mask[~np.isnan(mask)] = 1 
+                temp_mask[win,:,:] = mask
+
+                if np.all(temp_lidar[win,:,:]==0) or np.all(np.isnan(temp_lidar[win,:,:])):
+                    mask = np.zeros(temp_mask[win,:,:].shape)
+                    # mask[np.shape(mask)[0]//2,np.shape(mask)[1]//2]=1
+                    np.fill_diagonal(mask, 1)
+                    mask = np.flipud(mask)
+                    np.fill_diagonal(mask, 1)
+                    temp_mask[win,:,:] = mask
+
+            temp_mask = unblockshaped(temp_mask, rows,cols)
+
+
+
             s = unblockshaped(np.array(s), rows,cols)
             c = unblockshaped(np.array(c), rows,cols)
             rmse__ = unblockshaped(np.array(rmse__), rows,cols)
@@ -174,6 +196,10 @@ def rvog_inverse(args):
             count = unblockshaped(np.array(count), rows,cols)            
             temp_cor = unblockshaped(temp_cor, rows,cols)
 
+            s = s*temp_mask
+            c = c*temp_mask
+
+
         elif args.algo==2 and args.window_overlap>0:
             temp_cor, indices__ = emp.extract_patches(cor, patchsize=args.window_size, overlap=args.window_overlap)
             temp_lidar, indices__ = emp.extract_patches(lidar_ht_cal, patchsize=args.window_size, overlap=args.window_overlap)
@@ -320,8 +346,6 @@ def rvog_inverse(args):
             s50, c50, _, _, _ = zip(*results)
 
 
-
-
             s20 = unblockshaped(np.array(s20), rows,cols)
             c20 = unblockshaped(np.array(c20), rows,cols)
             s50 = unblockshaped(np.array(s50), rows,cols)
@@ -418,8 +442,13 @@ def rvog_inverse(args):
         t2 = time.time()
         print("[2/3] Generating height map... ",end=" ")
 
+
+
         ci = spatial_intp_lin(c)
         si = spatial_intp_lin(s)
+
+        ci = median_filter(ci, size=args.window_size//2)
+        si = median_filter(si, size=args.window_size//2)
         
         si[si>1.6]=1.6
         ci[ci>16]=16