diff --git a/ich.sh b/ich.sh index e15ec9e8f16cd51bb15448e5862aae0738c5b425..1afd3ae7bb5ce270b05c242a40ab40d168e4e49d 100644 --- a/ich.sh +++ b/ich.sh @@ -3,19 +3,20 @@ set -xuo pipefail basedir=$(dirname "$(readlink -f "$0")") -# ich_py="conda run --no-capture-output -n ich ${basedir}/src/ich/InSAR_canopy_height.py" -# ich_py="conda run --no-capture-output ${basedir}/src/ich/InSAR_canopy_height.py" ich_py="python ${basedir}/src/ich/InSAR_canopy_height.py" -correlationFile="$(ls input/*_corr.tif)" -cal_ht="$(ls input/*_lidar.tif)" +# Ensure there is at least one file of each type +correlationFile=$(ls input/*_corr.tif 2>/dev/null || { echo "No correlation files found"; exit 1; }) +cal_ht=$(ls input/*_lidar.tif 2>/dev/null || { echo "No lidar files found"; exit 1; }) +# Initialize options array options=() -[[ "${1:--}" != "-" ]] && options=("${options[@]}" --lower_limit "${1:--}") -[[ "${2:--}" != "-" ]] && options=("${options[@]}" --upper_limit "${2:--}") -[[ "${3:--}" != "-" ]] && options=("${options[@]}" --window "${3:--}") -[[ "${4:--}" != "-" ]] && options=("${options[@]}" --validation "${4:--}") -[[ "${4:--}" != "-" ]] && options=("${options[@]}" --overlap "${5:--}") +# Add options based on provided arguments +[[ "${1:-}" != "-" ]] && options+=("--lower_limit" "${1}") +[[ "${2:-}" != "-" ]] && options+=("--upper_limit" "${2}") +[[ "${3:-}" != "-" ]] && options+=("--window" "${3}") +[[ "${4:-}" != "-" ]] && options+=("--validation" "${4}") -${ich_py} --correlationFile "${correlationFile}" --cal_ht "${cal_ht}" "${options[@]}" +# Run the Python script with the determined options +${ich_py} --correlationFile "${correlationFile}" --cal_ht "${cal_ht}" "${options[@]}" \ No newline at end of file diff --git a/src/ich/InSAR_canopy_height.py b/src/ich/InSAR_canopy_height.py index 586a4f764b3de06c720aaebd537084701421d83e..f6eb711c0435b46a9bf0aadb8e4f6dfee39dcf20 100644 --- a/src/ich/InSAR_canopy_height.py +++ b/src/ich/InSAR_canopy_height.py @@ -5,48 +5,50 @@ Created on Thu Dec 27 14:40:29 2023 @author: Narayanarao """ - -import argparse,os,sys,warnings,time - -import numpy as np, pandas as pd +import argparse +import os +import sys +import warnings +import numpy as np +import pandas as pd from osgeo import gdal from scipy import interpolate import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.colors import Normalize - from scipy.stats import linregress from skimage.util.shape import view_as_blocks from mpl_toolkits.axes_grid1 import make_axes_locatable from scipy.interpolate import griddata from tqdm import tqdm from scipy.interpolate import interpn +from plotting import density_scatter, plot_data, rmse +from algo import arc_sinc,arc_sinc_vectorized, cal_ +from rst_io import read_bin, write_tif, write_bin +from utils import blockshaped, unblockshaped, unpad, spatial_intp_lin +from args_in import rvog_inverse + +# Configure environment and warnings gdal.UseExceptions() warnings.filterwarnings('ignore') warnings.filterwarnings('error') np.seterr(divide='ignore', invalid='ignore') +# Command-line argument parsing +parser = argparse.ArgumentParser(description="Generate canopy height map from InSAR coherence data.") -from plotting import density_scatter,plot_data,rmse -from algo import arc_sinc,arc_sinc_,cal_ -from rst_io import read_bin, write_tif, write_bin -from utils import blockshaped,unblockshaped,unpad,spatial_intp_lin -from args_in import rvog_inverse -########################################################################## -parser = argparse.ArgumentParser() - -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 (m)") -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("-w", "--window",dest = "window_size", default = 10, help="Size", type=int) -parser.add_argument("-val", "--validation",dest = "validation", default = 0, help="fraction to split cross validation", type=float) -parser.add_argument("-ol", "--overlap",dest = "window_overlap", default = 0, help="window overlap fraction", type=float) - +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("-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) args = parser.parse_args() - if __name__ == "__main__": - - rvog_inverse(args) + try: + rvog_inverse(args) + except Exception as e: + print(f"An error occurred: {e}") + sys.exit(1) diff --git a/src/ich/algo.py b/src/ich/algo.py index 775651a9f26e2919d1ac6aafbd47a1dc66413d33..72dc89e385c7ba91050f50ca34f92e8f6672c734 100644 --- a/src/ich/algo.py +++ b/src/ich/algo.py @@ -21,275 +21,211 @@ from mpl_toolkits.axes_grid1 import make_axes_locatable from scipy.interpolate import griddata from tqdm import tqdm from scipy.interpolate import interpn +from sklearn.metrics import mean_squared_error as mse +from concurrent.futures import ProcessPoolExecutor, as_completed +import os -def rmse(predictions, targets): - return np.sqrt(((predictions - targets) ** 2).mean()) - -def arc_sinc_(x, c_param): - # Get rid of extreme values by set all values where x > 1 equal to 1, and x < 0 equal to 0 - if x >1: - x=1 - elif x<0: - x=0 - # x[(x > 1)] = 1 - # x[(x < 0)] = 0 - - # Create array of increments between 0 and pi of size pi/100 +def arc_sinc(x, c_param): + x = np.clip(x, 0, 1) # Ensure x is between 0 and 1 XX = np.linspace(0, np.pi, num=100, endpoint=True) - - # Set the first value of XX to eps to avoid division by zero issues -> Paul's suggestion - XX[0] = np.spacing(1) - - # Calculate sinc for XX and save it to YY - ## YY = sinc(XX / math.pi) + XX[0] = np.spacing(1) # Avoid division by zero issues YY = np.sin(XX) / XX - - # Reset the first value of XX to zero and the first value of YY to the corresponding output XX[0] = 0 YY[0] = 1 - - # Set the last value of YY to 0 to avoid NaN issues YY[-1] = 0 - # Flip XX and YY left to right XX = XX[::-1] YY = YY[::-1] - # print(XX.shape,YY.shape) - # Run interpolation - # XX and YY are your original values, x is the query values, and y is the interpolated values that correspond to x - interp_func = interpolate.interp1d(YY, XX * c_param, kind='slinear') + interp_func = interpolate.interp1d(YY, XX * c_param, kind='slinear', bounds_error=False, fill_value=0) y = interp_func(x) - - # Set all values in y less than 0 equal to 0 - y[(y < 0)] = 0 - # return y - return y - -def arc_sinc(x, c_param): - # Get rid of extreme values by set all values where x > 1 equal to 1, and x < 0 equal to 0 - x[(x > 1)] = 1 - x[(x < 0)] = 0 - - # Create array of increments between 0 and pi of size pi/100 + return np.clip(y, 0, None) # Ensure no negative heights +def arc_sinc_vectorized(x, c_param): + # Ensure x and c_param are numpy arrays + x = np.asarray(x) + c_param = np.asarray(c_param) + + # Clip x to the range [0, 1] + x = np.clip(x, 0, 1) + + # Create array of increments between 0 and pi XX = np.linspace(0, np.pi, num=100, endpoint=True) - - # Set the first value of XX to eps to avoid division by zero issues -> Paul's suggestion + + # Set the first value of XX to eps to avoid division by zero issues XX[0] = np.spacing(1) - + # Calculate sinc for XX and save it to YY - ## YY = sinc(XX / math.pi) YY = np.sin(XX) / XX - - # Reset the first value of XX to zero and the first value of YY to the corresponding output + + # Reset the first value of XX and YY XX[0] = 0 YY[0] = 1 - - # Set the last value of YY to 0 to avoid NaN issues YY[-1] = 0 - # Flip XX and YY left to right + + # Reverse arrays for interpolation XX = XX[::-1] YY = YY[::-1] - # Run interpolation - # XX and YY are your original values, x is the query values, and y is the interpolated values that correspond to x - interp_func = interpolate.interp1d(YY, XX * c_param, kind='slinear') - y = interp_func(x) - - # Set all values in y less than 0 equal to 0 - y[(y < 0)] = 0 - # return y - return y -def dynamicWindow(cohArray, lidarArray, initial_ws, htl,htg): - rows, cols = cohArray.shape - c_parm = np.zeros((rows, cols)) # Initialize output array to store results - s_parm = np.zeros((rows, cols)) - rmse_parm = np.zeros((rows, cols)) - count = np.zeros((rows, cols)) - ht_ = np.zeros((rows, cols)) - - result = [] - i = 0 - while i < rows: - # print("%0.1f /100 "%((i/rows)*100),end="\r") - j = 0 - while j < cols: - # Reset to the initial block size for each block - S = initial_ws - # start_i = max(i - S // 2, 0) - # end_i = min(i + S // 2 + 1, rows) - # start_j = max(j - S // 2, 0) - # end_j = min(j + S // 2 + 1, cols) - - start_i = max(i - S // 2, 0) - end_i = min(i + S // 2, rows) - start_j = max(j - S // 2, 0) - end_j = min(j + S // 2, cols) - - while True: - lidarBlock = lidarArray[start_i:end_i, start_j:end_j] - cohBlock = cohArray[start_i:end_i, start_j:end_j] - - # Check if the block if there are 3 or more lidar samples - if np.isfinite(lidarBlock).any() and (np.count_nonzero(~np.isnan(lidarBlock))>3): - parm = cal_(cohBlock,lidarBlock, htl,htg) - mask = lidarBlock.copy() - mask[~np.isnan(mask)] = 1 - # mask[mask!=1] = 0 - - if np.all(np.array(parm) == 0): - parm=parm_.copy() - if lidarBlock.shape[0]*lidarBlock.shape[1]>initial_ws*initial_ws: - # 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) - - s_parm[start_i:end_i, start_j:end_j] = np.full(lidarBlock.shape, parm[1])*mask - c_parm[start_i:end_i, start_j:end_j] = np.full(lidarBlock.shape, parm[2])*mask - rmse_parm[start_i:end_i, start_j:end_j] = np.full(lidarBlock.shape, parm[4]) - count[start_i:end_i, start_j:end_j] = np.full(lidarBlock.shape, np.count_nonzero(~np.isnan(lidarBlock))) - - gama = cohBlock / parm[1] - ht_[start_i:end_i, start_j:end_j] = arc_sinc(gama, parm[2])*mask - parm_ = parm.copy() - # print(f"Block at ({i}, {j}) with size {lidarBlock.shape} : Result = {parm}") - - break - else: - # Increase the block size equally in all four directions - S += 2 - - # Update block boundaries while staying within the array limits - # start_i = max(i - S // 2, 0) - # end_i = min(i + S // 2 + 1, rows) - # start_j = max(j - S // 2, 0) - # end_j = min(j + S // 2 + 1, cols) - - start_i = max(i - S // 2, 0) - end_i = min(i + S // 2, rows) - start_j = max(j - S // 2, 0) - end_j = min(j + S // 2, cols) - - # If block size exceeds array dimensions, stop expanding - if start_i == 0 and end_i == rows and start_j == 0 and end_j == cols: - print(f"Unable to find finite block at position ({i}, {j}).") - break - - # Move to the next block position in the row - j += initial_ws + + # Initialize result array + y = np.zeros_like(x) + + # Handle the interpolation for each unique c_param value if needed + for unique_c in np.unique(c_param): + # Interpolation function for the current unique c_param + interp_func = interpolate.interp1d(YY, XX * unique_c, kind='slinear', bounds_error=False, fill_value=0) - # Move to the next block position in the column - i += initial_ws + # Apply the interpolation + mask = (c_param == unique_c) + y[mask] = interp_func(x[mask]) - return s_parm, c_parm, rmse_parm, ht_, count + # Clip the result to ensure non-negative values + return np.clip(y, 0, None) -def cal_(cor,gedi,htl,htg): - # scale = 1 +def cal_(cor, gedi, htl, htg): try: temp_cor = cor.copy() temp_gedi = gedi.copy() nn = np.count_nonzero(~np.isnan(temp_cor)) - result=[] - for S_param in np.arange(.2,1.6,0.1): - for C_param in np.arange(2,16,1): - # t2 = time.time() - gama = temp_cor / S_param - gama[gama>1] = 1 - gama[gama<0] = 0 + result = [] + + # Define ranges for parameters + S_range = np.arange(0.2, 1.6, 0.1) + C_range = np.arange(2, 16, 1) + + for S_param in S_range: + for C_param in C_range: + gama = np.clip(temp_cor / S_param, 0, 1) height = arc_sinc(gama, C_param) - height[np.isnan(gama)] = np.nan - height[height>=htg] = np.nan - height[height<=htl] = np.nan - del gama - - x, y = temp_gedi.flatten(),height.flatten() - del height - - xy = np.vstack([x,y]).T - del x,y - xy = xy[~np.isnan(xy).any(axis=1)] - x,y = xy[:,0], xy[:,1] - del xy - try: - slope, intercept, r_value, p_value, std_err = linregress(x,y) - RMSE = rmse(x,y) - # rd = rmsd(x,y) - # nrd = nrmsd(x,y) - # ubrd = ubrmsd(x,y) - except: - slope, intercept, r_value, p_value, std_err = np.nan,np.nan,np.nan,np.nan,np.nan - RMSE = np.nan - # rd = np.nan - # nrd = np.nan - # ubrd = np.nan - result.append([np.shape(x)[0],S_param,C_param,r_value,RMSE]) - - result = np.array(result) - 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.drop(tempdf[tempdf.N < nn*0.1].index, inplace=True) - tempdf.drop(tempdf[tempdf.N <= 3].index, inplace=True) - tempdf = tempdf.sort_values(by=['rmse'], ascending=True) - - # print(tid,', %.1f, N=%d,nn=%d, %0.2f, %.2f, %0.2f, %0.2f'%(u['class'].iloc[ii], - # tempdf.iloc[0]['N'],nn, tempdf.iloc[0]['S'],tempdf.iloc[0]['C'], - # tempdf.iloc[0]['r'],tempdf.iloc[0]['rmse'])) - - sCoarse = np.round(tempdf.iloc[0]['S'],2) - cCoarse = np.round(tempdf.iloc[0]['C'],2) + height[(height >= htg) | (height <= htl)] = np.nan + x, y = temp_gedi.flatten(), height.flatten() + valid = ~np.isnan(x) & ~np.isnan(y) + x, y = x[valid], y[valid] + + if len(x) > 0: + try: + slope, intercept, r_value, p_value, std_err = linregress(x, y) + RMSE = np.sqrt(np.mean((x - y) ** 2)) + except Exception as e: + # print(f"Linear regression error: {e}") + 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=[] - for S_param in np.arange(sCoarse-0.1,sCoarse+.1,0.02): - for C_param in np.arange(cCoarse-1,cCoarse+1,.2): - # t2 = time.time() - gama = temp_cor / S_param - gama[gama>1] = 1 - gama[gama<0] = 0 - height = arc_sinc(gama, C_param) - height[np.isnan(gama)] = np.nan - height[height>=htg] = np.nan - height[height<=htl] = np.nan - del gama - - x, y = temp_gedi.flatten(),height.flatten() - del height - - xy = np.vstack([x,y]).T - del x,y - xy = xy[~np.isnan(xy).any(axis=1)] - x,y = xy[:,0], xy[:,1] - del xy - try: - slope, intercept, r_value, p_value, std_err = linregress(x,y) - RMSE = rmse(x,y) - # rd = rmsd(x,y) - # nrd = nrmsd(x,y) - # ubrd = ubrmsd(x,y) - except: - slope, intercept, r_value, p_value, std_err = np.nan,np.nan,np.nan,np.nan,np.nan - RMSE = np.nan - # rd = np.nan - # nrd = np.nan - # ubrd = np.nan - - - # print('N: %d, S %0.2f, C %0.2f r = %0.2f, RMSE = %0.2f, ubRMSD = %0.2f, t=%0.2f'%(np.shape(x)[0], - # S_param,C_param,r_value,RMSE,ubrd, - # (time.time()-t2))) - result.append([int(np.shape(x)[0]),np.round(S_param,2),np.round(C_param,2),np.round(r_value,3),np.round(RMSE,3)]) result = np.array(result) - 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.sort_values(by=['rmse'], ascending=True) - # tempdf.drop(tempdf[tempdf.N < nn*0.1].index, inplace=True) - tempdf.drop(tempdf[tempdf.N <= 3].index, inplace=True) + 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 tempdf.empty: - tempdf = pd.DataFrame(data={'N':0,'S':0,'C':0, 'r':0,'rmse':0},index=[0]) - return list(tempdf.iloc[0]) + sCoarse = cCoarse = 0 + else: + sCoarse = np.round(tempdf.iloc[0]['S'], 2) + cCoarse = np.round(tempdf.iloc[0]['C'], 2) + + result = [] + for S_param in np.arange(sCoarse - 0.1, sCoarse + 0.1, 0.02): + for C_param in np.arange(cCoarse - 1, cCoarse + 1, 0.2): + gama = np.clip(temp_cor / S_param, 0, 1) + height = arc_sinc(gama, C_param) + height[(height >= htg) | (height <= htl)] = np.nan + x, y = temp_gedi.flatten(), height.flatten() + valid = ~np.isnan(x) & ~np.isnan(y) + x, y = x[valid], y[valid] + + if len(x) > 0: + try: + slope, intercept, r_value, p_value, std_err = linregress(x, y) + RMSE = np.sqrt(np.mean((x - y) ** 2)) + except Exception as e: + # print(f"Linear regression error: {e}") + 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), np.round(S_param, 2), np.round(C_param, 2), np.round(r_value, 3), np.round(RMSE, 3)]) + + result = np.array(result) + 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 tempdf.empty: + return [0, 0, 0, 0, 0] else: return list(tempdf.iloc[0]) except Exception as e: - tempdf = pd.DataFrame(data={'N':0,'S':0,'C':0, 'r':0,'rmse':0},index=[0]) - return list(tempdf.iloc[0]) + # print(f"Exception in cal_: {e}") + return [0, 0, 0, 0, 0] + +def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_): + rows, cols = cohArray.shape + S = initial_ws + start_i = max(i - S // 2, 0) + end_i = min(i + S // 2, rows) + start_j = max(j - S // 2, 0) + end_j = min(j + S // 2, cols) + + while True: + 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: + parm = cal_(cohBlock, lidarBlock, htl, htg) + mask = np.where(~np.isnan(lidarBlock), 1, 0) + + if np.all(np.array(parm) == 0): + parm = parm_.copy() + if lidarBlock.shape[0]*lidarBlock.shape[1]>initial_ws*initial_ws: + # 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) + + s_parm = np.full(lidarBlock.shape, parm[1]) * mask + c_parm = np.full(lidarBlock.shape, parm[2]) * mask + rmse_parm = np.full(lidarBlock.shape, parm[4]) + count = np.count_nonzero(~np.isnan(lidarBlock)) + + gama = cohBlock / parm[1] + ht = arc_sinc(gama, parm[2]) * mask + + return start_i, end_i, start_j, end_j, s_parm, c_parm, rmse_parm, ht, count + else: + S += 2 + start_i = max(i - S // 2, 0) + end_i = min(i + S // 2, rows) + start_j = max(j - S // 2, 0) + end_j = min(j + S // 2, cols) + + if start_i == 0 and end_i == rows and start_j == 0 and end_j == cols: + print(f"Unable to find finite block at position ({i}, {j}).") + return start_i, end_i, start_j, end_j, np.zeros_like(lidarBlock), np.zeros_like(lidarBlock), 0, np.zeros_like(lidarBlock), 0 + +def dynamicWindow(cohArray, lidarArray, initial_ws, htl, htg): + rows, cols = cohArray.shape + c_parm = np.zeros((rows, cols)) + s_parm = np.zeros((rows, cols)) + rmse_parm = np.zeros((rows, cols)) + count = np.zeros((rows, cols)) + ht_ = np.zeros((rows, cols)) + + parm_ = [0, 0, 0, 0, 0] + + num_workers = min(16, os.cpu_count()) + + with ProcessPoolExecutor(max_workers=num_workers) as executor: + futures = [] + for i in range(0, rows, initial_ws): + for j in range(0, cols, initial_ws): + futures.append(executor.submit(process_block, i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_)) + + for future in as_completed(futures): + start_i, end_i, start_j, end_j, s_p, c_p, r_p, ht, cnt = future.result() + s_parm[start_i:end_i, start_j:end_j] = s_p + c_parm[start_i:end_i, start_j:end_j] = c_p + rmse_parm[start_i:end_i, start_j:end_j] = r_p + ht_[start_i:end_i, start_j:end_j] = ht + count[start_i:end_i, start_j:end_j] = cnt + + return s_parm, c_parm, rmse_parm, ht_, count diff --git a/src/ich/args_in.py b/src/ich/args_in.py index 5c6aa30a11d7443a78a4c084859c8e7dd5767b21..69b07cd58d9ba9add38bbbe96505e92657d00304 100644 --- a/src/ich/args_in.py +++ b/src/ich/args_in.py @@ -5,15 +5,18 @@ Created on Thu Dec 27 14:40:29 2023 @author: Narayanarao """ -import argparse,os,sys,warnings,time - -import numpy as np, pandas as pd +import argparse +import os +import sys +import warnings +import time +import numpy as np +import pandas as pd from osgeo import gdal from scipy import interpolate import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.colors import Normalize - from scipy.stats import linregress from skimage.util.shape import view_as_blocks from mpl_toolkits.axes_grid1 import make_axes_locatable @@ -22,160 +25,124 @@ from tqdm import tqdm from scipy.interpolate import interpn from empatches import EMPatches from sklearn.model_selection import train_test_split +from plotting import density_scatter, plot_data, rmse +from algo import arc_sinc, arc_sinc_, arc_sinc_vectorized, cal_, dynamicWindow +from rst_io import read_bin, write_tif, write_bin +from utils import blockshaped, unblockshaped, unpad, spatial_intp_lin, split_sample + +# Initialize EMPatches instance emp = EMPatches() +# Configure GDAL and NumPy warnings gdal.UseExceptions() warnings.filterwarnings('ignore') warnings.filterwarnings('error') np.seterr(divide='ignore', invalid='ignore') - -from plotting import density_scatter,plot_data,rmse -from algo import 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 ########################################################################## def rvog_inverse(args): - + """ + Main function for processing InSAR canopy height estimation. + + Args: + args (argparse.Namespace): Command-line arguments containing file paths and parameters. + """ try: - t0 = time.time() - print( "-corFile {} -lidarFile {} -htl {} -htg {} -window_size {} -validation {} -window_overlap {}".format( - args.corFile, - args.lidarFile, - args.htl, - args.htg, - args.window_size, - args.validation, - args.window_overlap - + t0 = time.time() # Record start time + + # Print the arguments for debugging + print("-corFile {} -lidarFile {} -htl {} -htg {} -window_size {} -validation {}".format( + args.corFile, + args.lidarFile, + args.htl, + args.htg, + args.window_size, + args.validation )) - lidar_ht = read_bin(args.lidarFile) + # Read input files + lidar_ht = read_bin(args.lidarFile) # LiDAR height data corFile = args.corFile - cor = read_bin(corFile) - + cor = read_bin(corFile) # Correlation file - rPad,cPad = args.window_size-np.shape(lidar_ht)[0]%args.window_size,args.window_size-np.shape(lidar_ht)[1]%args.window_size + # Calculate padding needed to make dimensions divisible by window size + rPad, cPad = args.window_size - np.shape(lidar_ht)[0] % args.window_size, args.window_size - np.shape(lidar_ht)[1] % args.window_size pad_width = [(0, rPad), (0, cPad)] - if rPad!=0 or cPad!=0: - lidar_ht = np.pad(lidar_ht, pad_width, mode='constant',constant_values=0) - cor = np.pad(cor, pad_width, mode='constant',constant_values=0) + # Pad arrays to handle edge cases + if rPad != 0 or cPad != 0: + lidar_ht = np.pad(lidar_ht, pad_width, mode='constant', constant_values=0) + cor = np.pad(cor, pad_width, mode='constant', constant_values=0) + # Apply masks to filter out invalid or extreme values + lidar_ht[lidar_ht == 0] = np.nan + lidar_ht[lidar_ht < args.htl] = np.nan + lidar_ht[lidar_ht > args.htg] = np.nan + cor[cor < 0.1] = np.nan - lidar_ht[lidar_ht==0] = np.nan - lidar_ht[lidar_ht<args.htl] = np.nan - lidar_ht[lidar_ht>args.htg] = np.nan - cor[cor<0.1]=np.nan - rows,cols = np.shape(cor)[0],np.shape(cor)[1] - ###################################################### + rows, cols = np.shape(cor)[0], np.shape(cor)[1] - if args.validation>0 and args.validation<1: + ###################################################### + # Handle cross-validation if specified + if args.validation > 0 and args.validation < 1: temp_lidar = blockshaped(lidar_ht, args.window_size, args.window_size) lidar_ht_cal = np.zeros(np.shape(temp_lidar)) lidar_ht_val = np.zeros(np.shape(temp_lidar)) - - for win in range(np.shape(temp_lidar)[0]): - if np.count_nonzero(~np.isnan(temp_lidar[win,:,:])) > 10: - lidar_ht_cal[win,:,:],lidar_ht_val[win,:,:] = split_sample(temp_lidar[win,:,:],args.validation) + for win in range(np.shape(temp_lidar)[0]): + if np.count_nonzero(~np.isnan(temp_lidar[win, :, :])) > 10: + lidar_ht_cal[win, :, :], lidar_ht_val[win, :, :] = split_sample(temp_lidar[win, :, :], args.validation) else: - lidar_ht_cal[win,:,:] = temp_lidar[win,:,:] + lidar_ht_cal[win, :, :] = temp_lidar[win, :, :] + + lidar_ht_cal = unblockshaped(lidar_ht_cal, rows, cols) + lidar_ht_val = unblockshaped(lidar_ht_val, rows, cols) - lidar_ht_cal = unblockshaped(lidar_ht_cal, rows,cols) - lidar_ht_val = unblockshaped(lidar_ht_val, rows,cols) - - lidar_ht_cal[lidar_ht_cal==0]=np.nan - lidar_ht_val[lidar_ht_val==0]=np.nan + lidar_ht_cal[lidar_ht_cal == 0] = np.nan + lidar_ht_val[lidar_ht_val == 0] = np.nan - elif args.validation==0: + elif args.validation == 0: lidar_ht_cal = lidar_ht.copy() lidar_ht_val = lidar_ht.copy() else: print("\n Invalid validation fraction!!\n Please enter validation fraction between 0 and 1; \n -val=[0,1)") sys.exit(1) - # del lidar_ht ######################################################### + # Set up output directories and filenames os.chdir(os.path.dirname(corFile)) outDir = r'../output' if not os.path.exists(outDir): os.mkdir(outDir) - _postfix_str = '_w'+str(args.window_size)+'o'+"%02d" %(args.window_overlap*10)+"_"+str(args.htl)+str(args.htg) - - outcFile = os.path.join(outDir,os.path.basename(corFile).split('.tif')[0]+_postfix_str+'_c.tif') - outsFile = os.path.join(outDir,os.path.basename(corFile).split('.tif')[0]+_postfix_str+'_s.tif') - outciFile = os.path.join(outDir,os.path.basename(corFile).split('.tif')[0]+_postfix_str+'_ci.tif') - outsiFile = os.path.join(outDir,os.path.basename(corFile).split('.tif')[0]+_postfix_str+'_si.tif') - outrmseFile = os.path.join(outDir,os.path.basename(corFile).split('.tif')[0]+_postfix_str+'_rmse.tif') - - 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.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) - - print("[1/3] Generating calibration parameters...") - s=[] - c=[] - ht_=[] - rmse__=[] - count = [] - parm_ = [0,0,0,0,0] - for win in tqdm(range(np.shape(temp_cor)[0])): - # for win in (range(np.shape(temp_cor)[0])): - parm = cal_(temp_cor[win],temp_lidar[win],args.htl,args.htg) - - mask = temp_lidar[win].copy() - mask[~np.isnan(mask)] = 1 - - if np.all(np.array(parm) == 0): - parm=parm_.copy() - if np.all(mask==0) or np.all(np.isnan(mask)): - # 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) - - - s.append(np.full((args.window_size,args.window_size), parm[1])*mask) - c.append(np.full((args.window_size,args.window_size), parm[2])*mask) - # s.append(np.full((args.window_size,args.window_size), parm[1])) - # c.append(np.full((args.window_size,args.window_size), parm[2])) - - rmse__.append(np.full((args.window_size,args.window_size), parm[4])) - - count.append(np.count_nonzero(~np.isnan(temp_lidar[win,:,:]))) - gama = temp_cor[win] / parm[1] - ht_.append(arc_sinc(gama, parm[2])*mask) - - parm_ = parm.copy() - - - s = emp.merge_patches(s, indices__,mode='max') - c = emp.merge_patches(c, indices__,mode='max') - rmse__ = emp.merge_patches(rmse__, indices__,mode='max') - ht_ = emp.merge_patches(ht_, indices__,mode='max') - temp_cor = emp.merge_patches(temp_cor, indices__,mode='max') - count = emp.merge_patches(count, indices__,mode='mean') - - - else: - print("[1/3] Generating calibration parameters...") - s,c,rmse__,ht_, count = dynamicWindow(cor,lidar_ht_cal, args.window_size,args.htl,args.htg) - - ht_[ht_==0] = np.nan - s[s==0]=np.nan - c[c==0]=np.nan - rmse__[rmse__==0]=np.nan - - if rPad!=0 or cPad!=0: + # 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) + + # Define output file paths + outcFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_c.tif') + outsFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_s.tif') + outciFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_ci.tif') + outsiFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_si.tif') + outrmseFile = os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_rmse.tif') + + 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') + + # Generate calibration parameters + print("[1/3] Generating calibration parameters...") + s, c, rmse__, ht_, count = dynamicWindow(cor, lidar_ht_cal, args.window_size, args.htl, args.htg) + + # Mask out invalid values + ht_[ht_ == 0] = np.nan + s[s == 0] = np.nan + c[c == 0] = np.nan + rmse__[rmse__ == 0] = np.nan + + # Remove padding from arrays if applied + if rPad != 0 or cPad != 0: lidar_ht_cal = unpad(lidar_ht_cal, pad_width) c = unpad(c, pad_width) s = unpad(s, pad_width) @@ -184,57 +151,52 @@ def rvog_inverse(args): count = unpad(count, pad_width) temp_cor = unpad(cor, pad_width) - print(lidar_ht_cal.shape,ht_.shape) - print(np.nanmean(lidar_ht_cal),np.nanstd(lidar_ht_cal)) - print(np.nanmean(ht_),np.nanstd(ht_)) - - plot_data(lidar_ht_cal,ht_,args.htl,args.htg,'Lidar height(m)','InSAR inverted (m)',pltName) - + # Plot calibration results + plot_data(lidar_ht_cal, ht_, args.htl, args.htg, 'Lidar height (m)', 'InSAR inverted (m)', pltName) + # Interpolate calibration parameters ci = spatial_intp_lin(c) si = spatial_intp_lin(s) - - si[si>1.6]=1.6 - ci[ci>16]=16 - ci[ci<0]=0 - si[si<0]=0 - + + # Clip extreme values + si[si > 1.6] = 1.6 + ci[ci > 16] = 16 + ci[ci < 0] = 0 + si[si < 0] = 0 + + # Calculate gamma values and estimate height gama = (temp_cor / si) gama = gama.flatten() - gama[gama>1]=1 - gama[gama<0]=0 + gama[gama > 1] = 1 + gama[gama < 0] = 0 c_temp = ci.flatten() - c_temp[c_temp<0] = np.nan - height=[] - + c_temp[c_temp < 0] = np.nan + height = [] + print("[2/3] Generating height map...") - for ii in tqdm(range(np.size(c_temp))): - # for ii in (range(np.size(c_temp))): - height.append(arc_sinc_(gama[ii], c_temp[ii])) - + height = arc_sinc_vectorized(gama, c_temp) height = np.array(height).reshape(np.shape(temp_cor)) - height[height<0] = 0 - - if args.validation>0 and args.validation<1: + height[height < 0] = 0 + + # If validation is enabled, save validation results + if args.validation > 0 and args.validation < 1: lidar_ht_val = unpad(lidar_ht_val, pad_width) - write_tif(os.path.join(outDir,os.path.basename(corFile).split('.tif')[0]+_postfix_str+'_val.tif'), - 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) + write_tif(os.path.join(outDir, os.path.basename(corFile).split('.tif')[0] + _postfix_str + '_val.tif'), + 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) + # Write all output files print("[3/3] Writing output...") - 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)) - write_tif(outrmseFile,rmse__,os.path.basename(args.corFile)) - - 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)) + 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)) + write_tif(outrmseFile, rmse__, os.path.basename(args.corFile)) + 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("Completed in %0.2f sec"%(time.time()-t0)) + print("Completed in %0.2f sec" % (time.time() - t0)) except Exception as e: - print ("Task failed due to ", e) + print("Task failed due to", e) sys.exit(1)