Skip to content
Snippets Groups Projects
Commit 4bd09d3f authored by Narayanarao Bhogapurapu's avatar Narayanarao Bhogapurapu
Browse files

multi algorithm support

parent b0b7e207
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ich/InSAR_canopy_height.py
+ 25
− 26
View file @ 4bd09d3f
......@@ -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)
src/ich/algo.py
+ 64
− 82
View file @ 4bd09d3f
......@@ -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
src/ich/args_in.py
+ 311
− 236
View file @ 4bd09d3f
......@@ -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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment