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

multi-scale calibration

parent 4bd09d3f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/ich/algo.py
+ 29
− 19
View file @ a81fc8ca
...@@ -73,14 +73,16 @@ def process_st_window(win, temp_cor, temp_lidar, args): ...@@ -73,14 +73,16 @@ def process_st_window(win, temp_cor, temp_lidar, args):
# if np.all(np.array(parm) == 0): # if np.all(np.array(parm) == 0):
# parm = parm_.copy() # parm = parm_.copy()
if np.all(mask == 0) or np.all(np.isnan(mask)): # if np.all(mask == 0) or np.all(np.isnan(mask)):
np.fill_diagonal(mask, 1) # np.fill_diagonal(mask, 1)
mask = np.flipud(mask) # mask = np.flipud(mask)
np.fill_diagonal(mask, 1) # np.fill_diagonal(mask, 1)
mask = np.nan_to_num(mask) # mask = np.nan_to_num(mask)
s = np.full((args.window_size, args.window_size), parm[1]) * 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 # c = np.full((args.window_size, args.window_size), parm[2]) * mask
s = np.full((args.window_size, args.window_size), parm[1])
c = np.full((args.window_size, args.window_size), parm[2])
rmse = np.full((args.window_size, args.window_size), parm[4]) rmse = np.full((args.window_size, args.window_size), parm[4])
count = np.full((args.window_size, args.window_size), count = np.full((args.window_size, args.window_size),
np.count_nonzero(~np.isnan(temp_lidar))) np.count_nonzero(~np.isnan(temp_lidar)))
...@@ -186,6 +188,9 @@ def cal_(temp_cor, temp_gedi, htl, htg): ...@@ -186,6 +188,9 @@ def cal_(temp_cor, temp_gedi, htl, htg):
height[(height >= htg) | (height <= htl)] = np.nan height[(height >= htg) | (height <= htl)] = np.nan
inv_N = np.count_nonzero(~np.isnan(height))
x, y = temp_gedi.flatten(), height.flatten() x, y = temp_gedi.flatten(), height.flatten()
valid = ~np.isnan(x) & ~np.isnan(y) valid = ~np.isnan(x) & ~np.isnan(y)
x, y = x[valid], y[valid] x, y = x[valid], y[valid]
...@@ -206,14 +211,13 @@ def cal_(temp_cor, temp_gedi, htl, htg): ...@@ -206,14 +211,13 @@ def cal_(temp_cor, temp_gedi, htl, htg):
'r': result[:, 3], 'rmse': result[:, 4]}) 'r': result[:, 3], 'rmse': result[:, 4]})
del result del result
tempdf.dropna(subset=['rmse'], inplace=True) tempdf.dropna(subset=['rmse'], inplace=True)
if nn>6:
tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
if tempdf.empty: if tempdf.empty:
sCoarse = cCoarse = 0 return [0, 0, 0, 0, 0]
else: elif nn>6:
sCoarse = np.round(tempdf.iloc[0]['S'], 2) tempdf_coarse = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
cCoarse = np.round(tempdf.iloc[0]['C'], 2) sCoarse = np.round(tempdf_coarse.iloc[0]['S'], 2)
cCoarse = np.round(tempdf_coarse.iloc[0]['C'], 2)
del tempdf del tempdf
result = [] result = []
for S_param in np.arange(sCoarse - 0.1, sCoarse + 0.1, 0.02): for S_param in np.arange(sCoarse - 0.1, sCoarse + 0.1, 0.02):
...@@ -242,11 +246,17 @@ def cal_(temp_cor, temp_gedi, htl, htg): ...@@ -242,11 +246,17 @@ def cal_(temp_cor, temp_gedi, htl, htg):
'r': result[:, 3], 'rmse': result[:, 4]}) 'r': result[:, 3], 'rmse': result[:, 4]})
del result del result
tempdf.dropna(subset=['rmse'], inplace=True) tempdf.dropna(subset=['rmse'], inplace=True)
if nn>6: # if nn>6:
tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True) tempdf = tempdf[tempdf['N'] > 3].sort_values(by=['rmse'], ascending=True)
if tempdf.empty: if tempdf.empty and sCoarse==0:
return [0, 0, 0, 0, 0] return [0, 0, 0, 0, 0]
elif tempdf.empty and sCoarse!=0:
list(tempdf_coarse.iloc[0])
else:
return list(tempdf.iloc[0])
else: else:
return list(tempdf.iloc[0]) return list(tempdf.iloc[0])
......
src/ich/args_in.py
+ 161
− 26
View file @ a81fc8ca
...@@ -185,34 +185,26 @@ def rvog_inverse(args): ...@@ -185,34 +185,26 @@ def rvog_inverse(args):
rmse__=[] rmse__=[]
count = [] count = []
parm_ = [0,0,0,0,0] 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])) batch_size = 100 # Define your batch size
num_windows = np.shape(temp_cor)[0]
results = [None] * num_windows # Preallocate a list to hold results
count.append(np.count_nonzero(~np.isnan(temp_lidar[win,:,:]))) with concurrent.futures.ProcessPoolExecutor() as executor:
gama = temp_cor[win] / parm[1] futures = {}
ht_.append(arc_sinc(gama, parm[2])*mask) 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
parm_ = parm.copy() # Unpack results
s, c, rmse__, count, ht_ = zip(*results)
s = emp.merge_patches(s, indices__,mode='max') s = emp.merge_patches(s, indices__,mode='max')
...@@ -220,7 +212,7 @@ def rvog_inverse(args): ...@@ -220,7 +212,7 @@ def rvog_inverse(args):
rmse__ = emp.merge_patches(rmse__, indices__,mode='max') rmse__ = emp.merge_patches(rmse__, indices__,mode='max')
ht_ = emp.merge_patches(ht_, indices__,mode='max') ht_ = emp.merge_patches(ht_, indices__,mode='max')
temp_cor = emp.merge_patches(temp_cor, indices__,mode='max') temp_cor = emp.merge_patches(temp_cor, indices__,mode='max')
count = emp.merge_patches(count, indices__,mode='mean') count = emp.merge_patches(count, indices__,mode='avg')
elif args.algo==3: elif args.algo==3:
...@@ -256,10 +248,153 @@ def rvog_inverse(args): ...@@ -256,10 +248,153 @@ def rvog_inverse(args):
c = c*temp_mask c = c*temp_mask
temp_cor = cor.copy() temp_cor = cor.copy()
elif args.algo==4:
temp_cor = blockshaped(cor, args.window_size, args.window_size)
temp_lidar = blockshaped(lidar_ht_cal, args.window_size, args.window_size)
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
sw, cw, rmse__, count, ht_ = zip(*results)
del results
print('20')
args.window_size = 20
cor20 = blockshaped(cor, args.window_size, args.window_size)
lidar20 = blockshaped(lidar_ht_cal, args.window_size,args.window_size)
print('20 bolcked')
batch_size = 100 # Define your batch size
num_windows = np.shape(cor20)[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, cor20[win,:,:], lidar20[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
s20, c20, _, _, _ = zip(*results)
del results
print('50')
args.window_size = 50
cor50 = blockshaped(cor, 50, 50)
lidar50 = blockshaped(lidar_ht_cal, 50,50)
batch_size = 100 # Define your batch size
num_windows = np.shape(cor50)[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, cor50[win,:,:], lidar50[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
s50, c50, _, _, _ = zip(*results)
s20 = unblockshaped(np.array(s20), rows,cols)
c20 = unblockshaped(np.array(c20), rows,cols)
s50 = unblockshaped(np.array(s50), rows,cols)
c50 = unblockshaped(np.array(c50), rows,cols)
sw = unblockshaped(np.array(sw), rows,cols)
cw = unblockshaped(np.array(cw), 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)
# # s = np.nanmean(np.dstack([sw,s20,s50])
# s = np.nanmean(np.stack((sw,s20,s50), axis=0), axis=0)
# c = np.nanmean(np.stack((cw,c20,c50), axis=0), axis=0)
temp_lidar = blockshaped(lidar_ht_cal, 10, 10)
sw = blockshaped(sw, 10, 10)
cw = blockshaped(cw, 10, 10)
s20 = blockshaped(s20, 10, 10)
c20 = blockshaped(c20, 10, 10)
s50 = blockshaped(s50, 10, 10)
c50 = blockshaped(c50, 10, 10)
s = np.zeros(temp_lidar.shape)
c = np.zeros(temp_lidar.shape)
temp_mask = np.zeros(temp_lidar.shape)
for win in tqdm(range(np.shape(temp_lidar)[0])):
mask = temp_lidar[win,:,:].copy()
mask[~np.isnan(mask)] = 1
temp_mask[win,:,:] = mask
s[win,:,:] = sw[win,:,:]
c[win,:,:] = cw[win,:,:]
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
s[win,:,:] = np.nanmean(np.stack((sw[win,:,:],s20[win,:,:],s50[win,:,:]), axis=0), axis=0)
c[win,:,:] = np.nanmean(np.stack((cw[win,:,:],c20[win,:,:],c50[win,:,:]), axis=0), axis=0)
# del temp_lidar,sw,s20,s50,cw,c20,c50
temp_mask = unblockshaped(temp_mask, rows,cols)
s = unblockshaped(np.array(s), rows,cols)
c = unblockshaped(np.array(c), rows,cols)
s = s*temp_mask
c = c*temp_mask
# c = np.nanmean([sw,c20,c50])
else: else:
raise ValueError('Invalid algorithm type!') raise ValueError('Invalid algorithm type!')
ht_[ht_==0] = np.nan ht_[ht_==0] = np.nan
s[s==0]=np.nan s[s==0]=np.nan
c[c==0]=np.nan c[c==0]=np.nan
......
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