diff --git a/src/ich/__pycache__/algo.cpython-310.pyc b/src/ich/__pycache__/algo.cpython-310.pyc
deleted file mode 100644
index 646093577efcb78c07a364cf7c0abfe40b23e6a8..0000000000000000000000000000000000000000
Binary files a/src/ich/__pycache__/algo.cpython-310.pyc and /dev/null differ
diff --git a/src/ich/__pycache__/args_in.cpython-310.pyc b/src/ich/__pycache__/args_in.cpython-310.pyc
deleted file mode 100644
index 28a51a332587afa620a45f2625518e47917a22bf..0000000000000000000000000000000000000000
Binary files a/src/ich/__pycache__/args_in.cpython-310.pyc and /dev/null differ
diff --git a/src/ich/__pycache__/plotting.cpython-310.pyc b/src/ich/__pycache__/plotting.cpython-310.pyc
deleted file mode 100644
index 18e6f194e092817b5c38a3f9c23a05a84bf490af..0000000000000000000000000000000000000000
Binary files a/src/ich/__pycache__/plotting.cpython-310.pyc and /dev/null differ
diff --git a/src/ich/__pycache__/rst_io.cpython-310.pyc b/src/ich/__pycache__/rst_io.cpython-310.pyc
deleted file mode 100644
index 2518670d4097a85851354dd8f474301d7307a88e..0000000000000000000000000000000000000000
Binary files a/src/ich/__pycache__/rst_io.cpython-310.pyc and /dev/null differ
diff --git a/src/ich/__pycache__/utils.cpython-310.pyc b/src/ich/__pycache__/utils.cpython-310.pyc
deleted file mode 100644
index eb73f28196990210072efb2325d07b429e471dad..0000000000000000000000000000000000000000
Binary files a/src/ich/__pycache__/utils.cpython-310.pyc and /dev/null differ
diff --git a/src/ich/algo.py b/src/ich/algo.py
index b82061756e43ae9f14d4e6bcc4078f57fed5ebf4..24dfd42afef53f21e736abf3187b1a3dfc7467a4 100644
--- a/src/ich/algo.py
+++ b/src/ich/algo.py
@@ -248,6 +248,7 @@ def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_):
cohBlock = cohArray[start_i:end_i, start_j:end_j]
if np.count_nonzero(~np.isnan(cohBlock)) <5:
+ # print('no coh')
count = np.count_nonzero(~np.isnan(lidarBlock))
s_parm = 0
c_parm = 0
@@ -259,22 +260,23 @@ def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_):
return start_i, end_i, start_j, end_j, s_parm, c_parm, rmse_parm, ht, count
# elif (np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 1) or max(lidarBlock.shape)>512:
- elif (np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 1) or max(lidarBlock.shape)>512:
+ elif (np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 3) or max(lidarBlock.shape)>512:
parm = cal_(cohBlock, lidarBlock, htl, htg)
mask = np.where(~np.isnan(lidarBlock), 1, 0)
- if np.all(np.array(parm) == 0):
- parm = parm_.copy()
- del parm_
+ # if np.all(np.array(parm) == 0):
+ # parm = parm_.copy()
+ # del parm_
if lidarBlock.shape[0]*lidarBlock.shape[1]>initial_ws*initial_ws:
# mask[np.shape(mask)[0]//2,np.shape(mask)[1]//2]=1
+ print('blank_blocks',lidarBlock.shape,parm)
np.fill_diagonal(mask, 1)
mask = np.flipud(mask)
np.fill_diagonal(mask, 1)
- s_parm = np.full(lidarBlock.shape, parm[1]) * mask
- c_parm = np.full(lidarBlock.shape, parm[2]) * mask
+ s_parm = np.full(lidarBlock.shape, parm[1])
+ c_parm = np.full(lidarBlock.shape, parm[2])
rmse_parm = np.full(lidarBlock.shape, parm[4])
count = np.count_nonzero(~np.isnan(lidarBlock))
diff --git a/src/ich/args_in.py b/src/ich/args_in.py
index 498c3a1ca9601edde075cc162206e28a5a8a5170..e4891bdf3d071cf02f6c2cadfbc277251b70034c 100644
--- a/src/ich/args_in.py
+++ b/src/ich/args_in.py
@@ -28,7 +28,7 @@ from sklearn.model_selection import train_test_split
from plotting import density_scatter, plot_data, rmse
from algo import arc_sinc_fast,arc_sinc, arc_sinc_, cal_, dynamicWindow
from rst_io import read_bin, write_tif, write_bin
-from utils import blockshaped, unblockshaped, unpad, spatial_intp_lin, split_sample
+from utils import blockshaped, unblockshaped, unpad, spatial_intp_lin, split_sample, spatial_intp_idw
# Initialize EMPatches instance
emp = EMPatches()
@@ -148,6 +148,26 @@ def rvog_inverse(args):
t1 = time.time()
s, c, rmse__, ht_, count = dynamicWindow(cor, lidar_ht_cal, args.window_size, args.htl, args.htg)
+ temp_lidar = blockshaped(lidar_ht_cal, args.window_size, args.window_size)
+ 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)
+ 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 = s*temp_mask
+ c = c*temp_mask
+
+
# Mask out invalid values
ht_[ht_ == 0] = np.nan
s[s == 0] = np.nan
@@ -171,10 +191,18 @@ def rvog_inverse(args):
t2 = time.time()
print("[2/3] Generating height map... ",end=" ")
- # Interpolate calibration parameters
+ # # Interpolate calibration parameters
ci = spatial_intp_lin(c)
si = spatial_intp_lin(s)
+ # # Interpolate calibration parameters
+ # ci = spatial_intp_idw(c)
+ # si = spatial_intp_idw(s)
+
+ # # Interpolate calibration parameters
+ # ci = spatial_intp_kriging(c)
+ # si = spatial_intp_kriging(s)
+
# Clip extreme values
si[si > 1.6] = 1.6
ci[ci > 16] = 16
diff --git a/src/ich/utils.py b/src/ich/utils.py
index 337b4cb4308da0ba0de44956247e05e0735ccaf4..5a6882df9bfc490558e226694c5fba5fcc19513b 100644
--- a/src/ich/utils.py
+++ b/src/ich/utils.py
@@ -21,7 +21,7 @@ from scipy.interpolate import griddata
from tqdm import tqdm
from scipy.interpolate import interpn
from sklearn.model_selection import train_test_split
-
+# from pykrige.ok import OrdinaryKriging
def blockshaped(arr, nrows, ncols):
"""
@@ -73,6 +73,52 @@ def spatial_intp_lin(array):
return intp
+
+# def spatial_intp_kriging(array):
+# x = np.arange(0, array.shape[1])
+# y = np.arange(0, array.shape[0])
+# # Mask invalid values
+# array = np.ma.masked_invalid(array)
+# xx, yy = np.meshgrid(x, y)
+
+# # Get only the valid values
+# x1 = xx[~array.mask]
+# y1 = yy[~array.mask]
+# newarr = array[~array.mask]
+
+# # Perform ordinary kriging
+# OK = OrdinaryKriging(x1, y1, newarr, variogram_model='linear', verbose=False, enable_plotting=False)
+# z, ss = OK.execute('grid', x, y)
+
+# return z
+
+def spatial_intp_idw(array, power=0.5):
+ x = np.arange(0, array.shape[1])
+ y = np.arange(0, array.shape[0])
+ # Mask invalid values
+ array = np.ma.masked_invalid(array)
+ xx, yy = np.meshgrid(x, y)
+
+ # Get only the valid values
+ x1 = xx[~array.mask]
+ y1 = yy[~array.mask]
+ newarr = array[~array.mask]
+
+ # Function to compute IDW
+ def idw(xi, yi, x1, y1, values, power):
+ dist = np.sqrt((x1 - xi) ** 2 + (y1 - yi) ** 2)
+ # Handle case where distance is zero
+ if np.any(dist == 0):
+ return values[dist == 0][0]
+ weights = 1 / (dist ** power)
+ return np.sum(weights * values) / np.sum(weights)
+
+ # Apply IDW for each point in the grid
+ idw_values = np.array([[idw(i, j, x1, y1, newarr, power) for i in x] for j in y])
+
+ return idw_values
+
+
def split_sample(a,frac=0.3):
"""
Returns calibration and validation 2d arrays