247 lines
8.7 KiB
Python
247 lines
8.7 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""
|
|
HyTools: Hyperspectral image processing library
|
|
Copyright (C) 2021 University of Wisconsin
|
|
|
|
Authors: Adam Chlus, Zhiwei Ye, Philip Townsend.
|
|
|
|
This program is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, version 3 of the License.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
|
|
BRDF 校正
|
|
"""
|
|
import json
|
|
import ray
|
|
import numpy as np
|
|
import h5py
|
|
from .universal import universal_brdf,apply_universal
|
|
from .flex import flex_brdf,apply_flex,ndvi_stratify, get_kernel_samples, ndvi_bins, get_band_samples
|
|
from ..masks import mask_create
|
|
from ..misc import set_brdf, update_brdf, progbar
|
|
|
|
def apply_brdf_correct(hy_obj,data,dimension,index):
|
|
''' 在内存中应用 BRDF 校正。
|
|
'''
|
|
|
|
if hy_obj.brdf['type'] == 'universal':
|
|
data = apply_universal(hy_obj,data,dimension,index)
|
|
elif hy_obj.brdf['type'] == 'flex':
|
|
data = apply_flex(hy_obj,data,dimension,index)
|
|
elif hy_obj.brdf['type'] == 'local':
|
|
print('Local/class BRDF correction....under development')
|
|
return data
|
|
|
|
def load_brdf_precomputed(hy_obj,brdf_dict):
|
|
with open(brdf_dict['coeff_files'][hy_obj.file_name], 'r') as outfile:
|
|
hy_obj.brdf = json.load(outfile)
|
|
|
|
def set_solar_zn(hy_obj):
|
|
"""设置太阳天顶角归一化值"""
|
|
solar_zn = hy_obj.get_anc('solar_zn')
|
|
solar_zn = np.mean(solar_zn[hy_obj.mask['no_data']])
|
|
hy_obj.brdf['solar_zn_norm_radians'] = float(solar_zn)
|
|
return solar_zn
|
|
|
|
def ndvi_stratify_samples(combine_dict):
|
|
'''创建 NDVI 分箱分层掩膜
|
|
'''
|
|
|
|
ndvi = combine_dict["ndi_samples"]
|
|
class_mask = np.zeros(ndvi.shape)
|
|
|
|
for bin_num in combine_dict['brdf_dict']['bins']:
|
|
start,end = combine_dict['brdf_dict']['bins'][bin_num]
|
|
class_mask[(ndvi > start) & (ndvi <= end)] = bin_num
|
|
|
|
class_mask = class_mask.astype(np.int8)
|
|
combine_dict['ndvi_classes'] = class_mask
|
|
|
|
def get_topo_var_samples_pre(hy_obj):
|
|
'''获取分组地形校正变量，在 ndvi_stratify() 之后运行
|
|
'''
|
|
slope = hy_obj.get_anc('slope')
|
|
cosine_i = hy_obj.cosine_i()
|
|
sample_ind = (hy_obj.ancillary['ndvi_classes'] !=0)
|
|
|
|
return slope[sample_ind], cosine_i[sample_ind]
|
|
|
|
def calc_flex_single_post(combine_data_dict,brdf_dict,load_reflectance_mode):
|
|
|
|
combine_data_dict["brdf_dict"] = brdf_dict
|
|
bad_bands = combine_data_dict['bad_bands']
|
|
# 确定分箱维度并创建类别掩膜
|
|
if brdf_dict['bin_type'] == 'dynamic':
|
|
bins = ndvi_bins(combine_data_dict["ndi_samples"],brdf_dict)
|
|
# 更新分箱数量
|
|
#print(bins)
|
|
combine_data_dict["brdf_dict"]['num_bins']=len(bins) #hy_obj.brdf['num_bins'] = len(bins)
|
|
else:
|
|
bins = brdf_dict['bins']
|
|
|
|
combine_data_dict['brdf_dict']['bins'] = {k:v for (k,v) in enumerate(bins,start=1)}
|
|
|
|
ndvi_stratify_samples(combine_data_dict)
|
|
|
|
coeffs = {}
|
|
good_band_count=0
|
|
for band_num,band in enumerate(bad_bands):
|
|
if ~band:
|
|
coeffs[band_num] = {}
|
|
|
|
if load_reflectance_mode==0:
|
|
band_samples = combine_data_dict["reflectance_samples"][:,good_band_count] #ray.get([a.do.remote(get_band_samples,
|
|
#{'band_num':band_num}) for a in actors])
|
|
else:
|
|
combine_refl = []
|
|
for h5name in combine_data_dict["reflectance_samples"]:
|
|
h5_obj = h5py.File(h5name, "r")
|
|
sub_refl_samples = h5_obj["reflectance_samples"][()][:,good_band_count]
|
|
combine_refl += [sub_refl_samples]
|
|
h5_obj.close()
|
|
band_samples = np.concatenate(combine_refl,axis=0)
|
|
|
|
band_coeffs= []
|
|
for bin_num in combine_data_dict['brdf_dict']['bins']:
|
|
|
|
bin_mask = (combine_data_dict["ndvi_classes"]== bin_num)
|
|
|
|
X = np.concatenate([combine_data_dict["kernels_samples"],np.ones((bin_mask.shape[0],1))],axis=1)[bin_mask] #kernel_samples[:,:3][bin_mask]
|
|
y = band_samples[bin_mask]
|
|
band_coeffs.append(np.linalg.lstsq(X, y,rcond=-1)[0].flatten().tolist())
|
|
coeffs[band_num] = band_coeffs
|
|
progbar(np.sum(~bad_bands[:band_num+1]),np.sum(~bad_bands))
|
|
good_band_count+=1
|
|
|
|
print('\n')
|
|
|
|
combine_data_dict["brdf_dict"]['coeffs'] = coeffs
|
|
|
|
|
|
|
|
def calc_flex_single_pre(hy_obj,brdf_dict):
|
|
''' 获取单个图像的样本，用于未来的 BRDF 系数估计
|
|
'''
|
|
hy_obj.brdf['coeffs'] ={}
|
|
|
|
# 确定分箱维度并创建类别掩膜
|
|
if hy_obj.brdf['bin_type'] == 'dynamic':
|
|
bins = ndvi_bins(hy_obj.ndi()[hy_obj.mask['no_data']],brdf_dict)
|
|
# 更新分箱数量
|
|
hy_obj.brdf['num_bins'] = len(bins)
|
|
else:
|
|
bins = brdf_dict['bins']
|
|
|
|
hy_obj.brdf['bins'] = {k:v for (k,v) in enumerate(bins,start=1)}
|
|
ndvi_stratify(hy_obj)
|
|
kernel_samples= get_kernel_samples(hy_obj)
|
|
|
|
# 循环每个波段
|
|
refl_samples_list = []
|
|
used_band = []
|
|
for band_num,band in enumerate(hy_obj.bad_bands):
|
|
if ~band:
|
|
band_samples = hy_obj.do(get_band_samples, {'band_num':band_num})
|
|
refl_samples_list+=[band_samples[:,None]]
|
|
used_band+=[hy_obj.wavelengths[band_num]]
|
|
|
|
refl_samples = np.concatenate(refl_samples_list,axis=1)
|
|
|
|
slope_samples, cos_i_samples = get_topo_var_samples_pre(hy_obj) # 坡度和余弦i
|
|
|
|
return kernel_samples[:,:2], refl_samples, used_band, slope_samples, cos_i_samples
|
|
|
|
|
|
def calc_brdf_coeffs(actors,config_dict):
|
|
|
|
brdf_dict = config_dict['brdf']
|
|
|
|
if brdf_dict['type'] == 'precomputed':
|
|
print("使用预计算的 BRDF 系数")
|
|
_ = ray.get([a.do.remote(load_brdf_precomputed,
|
|
config_dict['brdf']) for a in actors])
|
|
else:
|
|
# 设置 BRDF 字典
|
|
_ = ray.get([a.do.remote(set_brdf,brdf_dict) for a in actors])
|
|
|
|
# 创建用于计算系数的掩膜
|
|
_ = ray.get([a.gen_mask.remote(mask_create,'calc_brdf',
|
|
brdf_dict['calc_mask']) for a in actors])
|
|
|
|
# 计算平均太阳天顶角
|
|
if isinstance(brdf_dict['solar_zn_type'],str):
|
|
|
|
# 分配每条线的平均太阳天顶角
|
|
solar_zn_samples = ray.get([a.do.remote(set_solar_zn) for a in actors])
|
|
# 计算并分配场景平均太阳天顶角
|
|
if brdf_dict['solar_zn_type'] == 'scene':
|
|
scene_mean = float(np.mean(solar_zn_samples))
|
|
_ = ray.get([a.do.remote(update_brdf,{'key':'solar_zn_norm_radians',
|
|
'value': scene_mean }) for a in actors])
|
|
print("场景平均太阳天顶角 : %s 度" % round(np.degrees(scene_mean),3))
|
|
|
|
elif isinstance(brdf_dict['solar_zn_type'],float):
|
|
_ = ray.get([a.do.remote(update_brdf,{'key':'solar_zn_norm_radians',
|
|
'value': brdf_dict['solar_zn_type']}) for a in actors])
|
|
else:
|
|
print('无法识别的太阳天顶角归一化')
|
|
|
|
print("计算 BRDF 系数")
|
|
if brdf_dict['type']== 'universal':
|
|
universal_brdf(actors,config_dict)
|
|
elif brdf_dict['type'] == 'flex':
|
|
flex_brdf(actors,config_dict)
|
|
elif brdf_dict['type'] == 'local':
|
|
print('本地/类别 BRDF 校正....开发中')
|
|
|
|
_ = ray.get([a.do.remote(lambda x: x.corrections.append('brdf')) for a in actors])
|
|
|
|
def calc_brdf_coeffs_pre(hy_obj,config_dict):
|
|
|
|
brdf_dict = config_dict['brdf']
|
|
|
|
if brdf_dict['type'] == 'precomputed':
|
|
print("使用预计算的 BRDF 系数")
|
|
load_brdf_precomputed(hy_obj,config_dict['brdf'])
|
|
|
|
else:
|
|
# 设置 BRDF 字典
|
|
|
|
set_brdf(hy_obj,brdf_dict)
|
|
set_solar_zn_0 = set_solar_zn(hy_obj)
|
|
|
|
# 创建用于计算系数的掩膜
|
|
hy_obj.gen_mask(mask_create,'calc_brdf',brdf_dict['calc_mask'])
|
|
|
|
kernel_samples, reflectance_samples, used_band, slope_samples, cos_i_samples = calc_flex_single_pre(hy_obj,brdf_dict)
|
|
|
|
hy_obj.corrections.append('brdf')
|
|
|
|
return {
|
|
"set_solar_zn":set_solar_zn_0,
|
|
#"ndvi":hy_obj.ndi(),
|
|
"kernel_samples":kernel_samples,
|
|
"reflectance_samples":reflectance_samples,
|
|
"used_band":used_band,
|
|
"slope_samples":slope_samples,
|
|
"cos_i_samples":cos_i_samples,
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|