diff --git a/src/ich/InSAR_canopy_height.py b/src/ich/InSAR_canopy_height.py index 66989fbd605a76703f8dfaa5d7e1d418578b6856..a700df55e96a0bad55951954d9e3e888c75a2447 100644 --- a/src/ich/InSAR_canopy_height.py +++ b/src/ich/InSAR_canopy_height.py @@ -5,50 +5,49 @@ Created on Thu Dec 27 14:40:29 2023 @author: Narayanarao """ -import argparse -import os -import sys -import warnings -import numpy as np -import pandas as pd + +import argparse,os,sys,warnings,time + +import numpy as np, 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.") -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=0, 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) +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("-al", "--algorithm",dest = "algo", default = 1, help="Algorithm Type", type=int) +parser.add_argument("-ol", "--overlap",dest = "window_overlap", default = 0, help="window overlap fraction", type=float) + args = parser.parse_args() + if __name__ == "__main__": - try: - rvog_inverse(args) - except Exception as e: - print(f"An error occurred: {e}") - sys.exit(1) + + rvog_inverse(args) diff --git a/src/ich/algo.py b/src/ich/algo.py index 0b299b8347f29c25398901c1fd88319b0c9d6498..ed8874965d911d95e20dc8ed3a8d44d60ecadae1 100644 --- a/src/ich/algo.py +++ b/src/ich/algo.py @@ -23,12 +23,22 @@ 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 concurrent.futures import os from concurrent.futures import ThreadPoolExecutor from scipy.interpolate import PchipInterpolator import gc import tracemalloc +from utils import blockshaped,unblockshaped + + + + + + + + # Precompute XX and YY XX = np.linspace(0, np.pi, num=100, endpoint=True) XX[0] = np.spacing(1) # Avoid division by zero @@ -51,6 +61,38 @@ def arc_sinc_fast(x, c_param): # Clip final result to the range [0, inf), works for both 1D and 2D arrays return np.maximum(y, 0) +################################################# +################################################# + +def process_st_window(win, temp_cor, temp_lidar, args): + parm = cal_(temp_cor, temp_lidar, args.htl, args.htg) + + mask = temp_lidar.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)): + np.fill_diagonal(mask, 1) + mask = np.flipud(mask) + np.fill_diagonal(mask, 1) + + mask = np.nan_to_num(mask) + s = np.full((args.window_size, args.window_size), parm[1]) * mask + c = np.full((args.window_size, args.window_size), parm[2]) * mask + rmse = np.full((args.window_size, args.window_size), parm[4]) + count = np.full((args.window_size, args.window_size), + np.count_nonzero(~np.isnan(temp_lidar))) + + gama = temp_cor / parm[1] + ht = arc_sinc(gama, parm[2]) * mask + + return s, c, rmse, count, ht +################################################# +################################################## +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 @@ -91,62 +133,39 @@ def arc_sinc_(x, c_param): # return y return y - 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) - XX[0] = np.spacing(1) # Avoid division by zero issues - YY = np.sin(XX) / XX - XX[0] = 0 - YY[0] = 1 - YY[-1] = 0 - XX = XX[::-1] - YY = YY[::-1] - interp_func = interpolate.interp1d(YY, XX * c_param, kind='slinear', bounds_error=False, fill_value=0) - y = interp_func(x) - return np.clip(y, 0, None) # Ensure no negative heights - + # 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 -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 + # Create array of increments between 0 and pi of size pi/100 XX = np.linspace(0, np.pi, num=100, endpoint=True) - - # Set the first value of XX to eps to avoid division by zero issues + + # 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) YY = np.sin(XX) / XX - - # Reset the first value of XX and YY + + # Reset the first value of XX to zero and the first value of YY to the corresponding output XX[0] = 0 YY[0] = 1 - YY[-1] = 0 - # Reverse arrays for interpolation + # 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] - - # 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) - - # Apply the interpolation - mask = (c_param == unique_c) - y[mask] = interp_func(x[mask]) - - # Clip the result to ensure non-negative values - return np.clip(y, 0, None) + # 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 cal_(temp_cor, temp_gedi, htl, htg): try: @@ -380,40 +399,3 @@ def process_batch(batch_futures, cohArray, lidarArray, initial_ws, htl, htg, par # Invoke garbage collection gc.collect() return results - - -## blocks (windows) are processed in parallel -# 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 = os.cpu_count()-1 - -# futures = [] -# with ProcessPoolExecutor(max_workers=num_workers) as executor: -# 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_)) - -# # Initialize the progress bar with the total number of futures -# with tqdm(total=len(futures)) as pbar: -# completed_jobs = 0 -# 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 - -# completed_jobs += 1 -# if completed_jobs % 100 == 0: # Update every 100 jobs -# pbar.update(100) - -# return s_parm, c_parm, rmse_parm, ht_, count diff --git a/src/ich/args_in.py b/src/ich/args_in.py index 81ea4a74b9c36b2f0932e5cfea0d6ff196312675..a77b52fc70e580dae7dd407c1c1ebca2a1a9be00 100644 --- a/src/ich/args_in.py +++ b/src/ich/args_in.py @@ -5,266 +5,341 @@ Created on Thu Dec 27 14:40:29 2023 @author: Narayanarao """ -import argparse -import os -import sys -import warnings -import time -import numpy as np -import pandas as pd +import argparse,os,sys,warnings,time + +import numpy as np, 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 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 empatches import EMPatches 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, spatial_intp_idw - -# Initialize EMPatches instance +import concurrent.futures 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_ +from algo import arc_sinc_fast,arc_sinc, arc_sinc_, cal_, dynamicWindow +from algo import process_st_window +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() # Record start time - - # Print the arguments for debugging - print("-corFile {} -lidarFile {} -htl {} -htg {} -window_size {} -validation {}".format( + + try: + t0 = time.time() + print( "-corFile {} -lidarFile {} -htl {} -htg {} -window_size {} -validation {} -algo {} -window_overlap {}".format( args.corFile, args.lidarFile, args.htl, args.htg, args.window_size, - args.validation - )) - - # Read input files - lidar_ht = read_bin(args.lidarFile) # LiDAR height data - corFile = args.corFile - cor = read_bin(corFile) # Correlation file - - # 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)] - - # 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 - - if args.htl==0 and args.htg==0: - """ 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 - - rows, cols = np.shape(cor)[0], np.shape(cor)[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) - else: - 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[lidar_ht_cal == 0] = np.nan - lidar_ht_val[lidar_ht_val == 0] = np.nan - - 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)") + args.validation, + args.algo, + args.window_overlap + + )) + + lidar_ht = read_bin(args.lidarFile) + corFile = args.corFile + cor = read_bin(corFile) + + + 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) + + + 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] + ###################################################### + + 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) + else: + 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[lidar_ht_cal==0]=np.nan + lidar_ht_val[lidar_ht_val==0]=np.nan + + 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 + ######################################################### + 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.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... ") + t1 = time.time() + + + if args.algo==1: + temp_cor = blockshaped(cor, args.window_size, args.window_size) + temp_lidar = blockshaped(lidar_ht_cal, args.window_size, args.window_size) + + s = np.zeros(np.shape(temp_cor)) + c = np.zeros(np.shape(temp_cor)) + ht_ = np.zeros(np.shape(temp_cor)) + rmse__ = np.zeros(np.shape(temp_cor)) + count = np.zeros(np.shape(temp_cor)) + + + parm_ = [0,0,0,0,0] + # print("[1/3] Generating calibration parameters...") + + batch_size = 100 # Define your batch size + num_windows = np.shape(temp_cor)[0] + results = [None] * num_windows # Preallocate a list to hold results + + with concurrent.futures.ProcessPoolExecutor() as executor: + futures = {} + for start in range(0, num_windows, batch_size): + end = min(start + batch_size, num_windows) + for win in range(start, end): + futures[executor.submit(process_st_window, win, temp_cor[win,:,:], temp_lidar[win,:,:], args)] = win + + for future in tqdm(concurrent.futures.as_completed(futures),total=len(futures)): + win_index = futures[future] # Get the corresponding index + result = future.result() + results[win_index] = result # Store result directly at the corresponding index + + # Unpack results + s, c, rmse__, count, ht_ = zip(*results) + + s = unblockshaped(np.array(s), rows,cols) + c = unblockshaped(np.array(c), rows,cols) + rmse__ = unblockshaped(np.array(rmse__), rows,cols) + ht_ = unblockshaped(np.array(ht_), rows,cols) + count = unblockshaped(np.array(count), rows,cols) + temp_cor = unblockshaped(temp_cor, rows,cols) + + 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) + + # 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') + + + elif args.algo==3: + 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) + + """ uncomment below to get cal coefficents uniformly across all windows """ + # for win in tqdm(range(np.shape(temp_lidar)[0])): + # 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 + + + 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 = s*temp_mask + c = c*temp_mask + + temp_cor = cor.copy() + + else: + raise ValueError('Invalid algorithm type!') + + 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: + lidar_ht_cal = unpad(lidar_ht_cal, pad_width) + temp_cor = unpad(temp_cor, pad_width) + s = unpad(s, pad_width) + c = unpad(c, pad_width) + ht_ = unpad(ht_, pad_width) + rmse__ = unpad(rmse__, pad_width) + count = unpad(count, 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=" ") + + 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 + + # Calculate gamma values and estimate height + gama = (temp_cor / si) + gama = gama.flatten() + gama[gama > 1] = 1 + gama[gama < 0] = 0 + c_temp = ci.flatten() + c_temp[c_temp < 0] = np.nan + height = [] + + # height = arc_sinc_vectorized(gama, c_temp) + # height = arc_sinc_vectorized_parallel(gama, c_temp) + height = arc_sinc_fast(gama, c_temp) + + # 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 = np.array(height).reshape(np.shape(temp_cor)) + height[height < 0] = 0 + print("%0.2f sec"%(time.time()-t2)) + + + 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) + + t3 = time.time() + # Write all output files + 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)) + 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("%0.2f sec"%(time.time()-t3)) + print("\n Total computing time : %0.2f sec" % (time.time() - t0)) + + + + except Exception as e: + print ("Task failed due to ", e) sys.exit(1) - - ######################################################### - # Set up output directories and filenames - os.chdir(os.path.dirname(corFile)) - outDir = r'../output' - if not os.path.exists(outDir): - os.mkdir(outDir) - - # Create postfix string for output files - _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') - 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') - - 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... ",end=" ") - 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) - - """ uncomment below to get cal coefficents uniformly across all windows """ - # for win in tqdm(range(np.shape(temp_lidar)[0])): - # 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 - - - 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 = s*temp_mask - c = c*temp_mask - - - # 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) - ht_ = unpad(ht_, pad_width) - rmse__ = unpad(rmse__, pad_width) - 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) - 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 - 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 - c_temp = ci.flatten() - c_temp[c_temp < 0] = np.nan - height = [] - print('...') - # height = arc_sinc_vectorized(gama, c_temp) - # height = arc_sinc_vectorized_parallel(gama, c_temp) - height = arc_sinc_fast(gama, c_temp) - - # 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 = 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: - 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) - - - t3 = time.time() - # Write all output files - 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)) - 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("%0.2f sec"%(time.time()-t3)) - - print("\n Total computing time : %0.2f sec" % (time.time() - t0)) - - # except Exception as e: - # print("Task failed due to", e) - # sys.exit(1)