修改正射插值模块,改为3次立方插值

Mn
tian jiax 2023-09-18 16:39:59 +08:00
parent 10b484e4bf
commit 477745e863
7 changed files with 109 additions and 70 deletions

View File

@ -1589,6 +1589,18 @@ class IndirectOrthorectification(Orthorectification):
print(exe_cmd)
print(os.system(exe_cmd))
print("==========================================================================")
def calInterpolation_cubic_Wgs84_rc_sar_sigma(self, parameter_path, dem_rc, in_sar, out_sar):
'''
# std::cout << "mode 10";
# std::cout << "SIMOrthoProgram.exe 9 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path";
'''
exe_path = r".\baseTool\x64\Release\SIMOrthoProgram.exe"
exe_cmd = r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path, 9, parameter_path,
dem_rc, in_sar, out_sar)
print(exe_cmd)
print(os.system(exe_cmd))
print("==========================================================================")
def getPowerTif(self,in_ori_path,out_power_path):
'''

View File

@ -486,6 +486,12 @@ class OrthoMain:
pass
def cut_dem(self, dem_merged_path, meta_file_path):
left_up_lon = 0
left_up_lat = 0
def RD_process_handle(self):
# RPC
@ -498,6 +504,7 @@ class OrthoMain:
out_dem_path = self.__workspace_ResampledDEM_path
dem_merged_path=DEMProcess.dem_merged(in_dem_path, meta_file_path, out_dem_path) # 生成TestDEM\mergedDEM_VRT.tif
# self.cut_dem(dem_merged_path, meta_file_path)
# 2、间接定位法求解行列坐标
slc_paths = self.__in_processing_paras["SLC"]
# 2.1 生成映射表
@ -519,11 +526,11 @@ class OrthoMain:
Orthorectification.preCaldem_sar_rc(dem_merged_path,in_slc_path,self.__workspace_Temporary_path,self.__workspace_package_path.replace("\\","\\\\")) # 初步筛选坐标范围
logger.info('progress bar: 40%')
clip_dem_reample_path=os.path.join(self.__workspace_Temporary_path, "SAR_dem.tiff")
infooption=gdal.InfoOptions("-json")
clip_dem_tif_info=gdal.Info(clip_dem_reample_path,options=infooption)
dem_merged_info=gdal.Info(dem_merged_path,options=infooption)
sampling_f=clip_dem_tif_info['size'][0]/dem_merged_info['size'][0]
# clip_dem_reample_path=os.path.join(self.__workspace_Temporary_path, "SAR_dem.tiff")
# infooption=gdal.InfoOptions("-json")
# clip_dem_tif_info=gdal.Info(clip_dem_reample_path,options=infooption)
# dem_merged_info=gdal.Info(dem_merged_path,options=infooption)
# sampling_f=clip_dem_tif_info['size'][0]/dem_merged_info['size'][0]
out_dir_path=self.__workspace_package_path.replace("\\","\\\\")
this_outSpace_path = out_dir_path
@ -565,6 +572,9 @@ class OrthoMain:
# GTC 入射角
GTC_rc_path=os.path.join(self.__workspace_package_path,"ori_sim-ortho.tif")
GTC_out_path=self.__workspace_package_path
parameter_path = os.path.join(self.__workspace_package_path, "orth_para.txt")
dem_rc = os.path.join(self.__workspace_Temporary_path, "dem_rc.tiff")
in_tif_paths = list(glob.glob(os.path.join(slc_paths, '*.tiff')))
for in_tif_path in in_tif_paths:
@ -577,7 +587,8 @@ class OrthoMain:
temp_slc_path=temp_slc_path.replace("_db.tif","-ortho.tif")
#inter_Range2Geo(self,lon_lat_path , data_tiff , grid_path , space)
Orthorectification.inter_Range2Geo(GTC_rc_path,out_power_path,temp_slc_path,Orthorectification.heightspace)
# Orthorectification.inter_Range2Geo(GTC_rc_path,out_power_path,temp_slc_path,Orthorectification.heightspace)
Orthorectification.calInterpolation_cubic_Wgs84_rc_sar_sigma(parameter_path, dem_rc, out_power_path, temp_slc_path)
break
#Orth_Slc.append(temp_slc_path)
# power_list.append(out_power_path)
@ -667,7 +678,7 @@ if __name__ == '__main__':
except Exception:
logger.exception("run-time error!")
finally:
OrthoMain.del_temp_workspace()
# OrthoMain.del_temp_workspace()
pass
end = datetime.datetime.now()
logger.info('running use time: %s ' % (end - start))

View File

@ -327,8 +327,10 @@ class ScatteringMain:
os.remove(this_out_local_inc_angle_rpc_path)
this_out_ori_sim_tiff = os.path.join(out_dir_path, "RD_ori_sim.tif") # out_dir_path + "\\" + "RD_ori_sim.tif"#// 局地入射角
this_in_rpc_lon_lat_path = this_out_ori_sim_tiff
this_in_rpc_lon_lat_path = this_out_ori_sim_tiff
parameter_path = os.path.join(self.__workspace_processing_path, "orth_para.txt")
dem_rc = os.path.join(self.__workspace_preprocessing_path, "dem_rc.tiff")
for in_tif_path in in_tif_paths:
out_tif_path = os.path.join(self.__workspace_preprocessing_path,os.path.splitext(os.path.basename(in_tif_path))[0]) + r"_DB.tif"
@ -344,7 +346,9 @@ class ScatteringMain:
# 移动RPC
#rpc_correction(in_tif_path,rpc_path,out_tif_path,dem_tif_file = None)
Orthorectification.inter_Range2Geo(this_in_rpc_lon_lat_path,out_tif_path,tempout_tif_path,Orthorectification.heightspace)
# Orthorectification.inter_Range2Geo(this_in_rpc_lon_lat_path,out_tif_path,tempout_tif_path,Orthorectification.heightspace)
Orthorectification.calInterpolation_cubic_Wgs84_rc_sar_sigma(parameter_path, dem_rc, out_tif_path,
tempout_tif_path)
#shutil.move(rpc_path,out_tif_path.replace(".tiff",".rpc"))
self.imageHandler.write_quick_view(tempout_tif_path, color_img=False)
else:
@ -384,9 +388,18 @@ class ScatteringMain:
para_dict.update({"imageinfo_ProductLevel": "LEVEL3"})
para_dict.update({"ProductProductionInfo_BandSelection": "1,2"})
para_dict.update({"ProductProductionInfo_AuxiliaryDataDescription": "DEM"})
para_dict.update({"MetaInfo_UnitDes": "DB"}) # 设置单位
CreateProductXml(para_dict, model_path, meta_xml_path).create_standard_xml()
temp_folder = os.path.join(self.__workspace_path, EXE_NAME, 'Output')
out_xml = os.path.join(temp_folder, os.path.basename(meta_xml_path))
if os.path.exists(temp_folder) is False:
os.mkdir(temp_folder)
# CreateProductXml(para_dict, model_path, out_xml).create_standard_xml()
shutil.copy(meta_xml_path, out_xml)
self.make_targz(self.__out_para, self.__workspace_processing_path)
logger.info('process_handle finished!')
logger.info('progress bar: 100%')
return True

View File

@ -1587,6 +1587,18 @@ class IndirectOrthorectification(Orthorectification):
print(exe_cmd)
print(os.system(exe_cmd))
print("==========================================================================")
def calInterpolation_cubic_Wgs84_rc_sar_sigma(self, parameter_path, dem_rc, in_sar, out_sar):
'''
# std::cout << "mode 10";
# std::cout << "SIMOrthoProgram.exe 9 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path";
'''
exe_path = r".\baseTool\x64\Release\SIMOrthoProgram.exe"
exe_cmd = r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path, 9, parameter_path,
dem_rc, in_sar, out_sar)
print(exe_cmd)
print(os.system(exe_cmd))
print("==========================================================================")
def getPowerTif(self,in_ori_path,out_power_path):
'''

View File

@ -15,6 +15,7 @@ from tool.algorithm.xml.AlgXmlHandle import ManageAlgXML, CheckSource, InitPara
from tool.algorithm.algtools.logHandler import LogHandler
from tool.algorithm.image.ImageHandle import ImageHandler
from tool.algorithm.block.blockprocess import BlockProcess
from tool.algorithm.xml.CreateMetaDict import CreateMetaDict, CreateProductXml
from tool.config.ConfigeHandle import Config as cf
from tool.algorithm.xml.CreatMetafile import CreateMetafile
from tool.algorithm.algtools.ROIAlg import ROIAlg as roi
@ -91,7 +92,12 @@ class LeafIndexMain:
self.__create_work_space()
self.__task_id = self.__alg_xml_handler.get_task_id()
self.__input_paras = self.__alg_xml_handler.get_input_paras() # 获取输入文件夹中的数据名、类型、地址
self.__out_para = os.path.join(self.__workspace_path,EXE_NAME, 'Output', r"LeafAreaIndexProduct.tar.gz")
# self.__out_para = os.path.join(self.__workspace_path,EXE_NAME, 'Output', r"LeafAreaIndexProduct.tar.gz")
SrcImageName = os.path.split(self.__input_paras["SAR"]['ParaValue'])[1].split('.tar.gz')[0]
result_name = SrcImageName + "-lef.tar.gz"
# out_name = os.path.splitext(os.path.splitext(os.path.basename(self.__input_paras['SLC']['ParaValue']))[0])[0]
# self.__out_para = os.path.join(self.__workspace_path, EXE_NAME, 'Output', "BackScatteringProduct.tar.gz")
self.__out_para = os.path.join(self.__workspace_path, EXE_NAME, 'Output', result_name)
self.__alg_xml_handler.write_out_para("LeafAreaIndexProduct", self.__out_para) # 写入输出参数
Polarization = self.__input_paras['Polarization']['ParaValue']
@ -315,9 +321,9 @@ class LeafIndexMain:
block_size = self.BlockProcess.get_block_size(self.__rows, self.__cols)
self.BlockProcess.cut(self.__workspace_maskcai_image_path, self.__workspace_cai_sartif_path, ['tif', 'tiff'], 'tif', block_size)
self.BlockProcess.cut(self.__workspace_maskcai_localang_path, self.__workspace_cai_localang_path, ['tif', 'tiff'], 'tif', block_size)
self.BlockProcess.cut(self.__workspace_maskcai_SoilMoi_path, self.__workspace_cai_SoilMoi_path, ['tif', 'tiff'], 'tif', block_size)
self.BlockProcess.cut_new(self.__workspace_maskcai_image_path, self.__workspace_cai_sartif_path, ['tif', 'tiff'], 'tif', block_size)
self.BlockProcess.cut_new(self.__workspace_maskcai_localang_path, self.__workspace_cai_localang_path, ['tif', 'tiff'], 'tif', block_size)
self.BlockProcess.cut_new(self.__workspace_maskcai_SoilMoi_path, self.__workspace_cai_SoilMoi_path, ['tif', 'tiff'], 'tif', block_size)
logger.info('mask data image success!')
# 计算每小块的叶面积指数
@ -349,14 +355,16 @@ class LeafIndexMain:
logger.info('progress bar: 80%')
ref_tif_path = self.__ref_img_path
date_type = 'float32'
self.BlockProcess.combine(self.__workspace_Leafindex_path, self.__cols, self.__rows, self.__workspace_processing_path + "combine", 'tif', ['tif'], date_type)
self.BlockProcess.combine_new(self.__workspace_Leafindex_path, self.__cols, self.__rows, self.__workspace_processing_path + "combine", 'tif', ['tif'], date_type)
# 添加地理信息
out_path = self.__workspace_processing_path + "combine\\LeafAreaIndexProduct.tif"
self.BlockProcess.assign_spatial_reference_byfile(ref_tif_path, out_path)
# 截取roi区域
bare_land_mask_path = self.cal_roi()
product_path = self.__product_dic + "LeafAreaIndexProduct.tif"
SrcImageName = os.path.split(self.__input_paras["SAR"]['ParaValue'])[1].split('.tar.gz')[0] + '-lef.tif'
# product_path = self.__product_dic + "LeafAreaIndexProduct.tif"
product_path = os.path.join(self.__product_dic, SrcImageName)
roi.cal_roi(product_path, out_path, bare_land_mask_path, background_value=np.nan)
logger.info('block_process finished!')
@ -374,28 +382,41 @@ class LeafIndexMain:
image_path = product_path
out_path1 = os.path.join(tem_folder, "trans_geo_projcs.tif")
out_path2 = os.path.join(tem_folder, "trans_projcs_geo.tif")
par_dict = CreateDict(image_path, out_path1, out_path2).calu_nature(start)
model_xml_path = os.path.join(tem_folder, "creat_standard.meta.xml") # 输出xml路径
id_min = 0
id_max = 0
for id in cover_id_list:
if id < id_min:
id_min = id
if id > id_max:
id_max = id
set_threshold = [id_max, id_min, threshold_of_ndvi_min, threshold_of_ndvi_max]
CreateStadardXmlFile(xml_path, self.alg_xml_path, par_dict, set_threshold, model_xml_path).create_standard_xml()
# par_dict = CreateDict(image_path, out_path1, out_path2).calu_nature(start)
# model_xml_path = os.path.join(tem_folder, "creat_standard.meta.xml") # 输出xml路径
#
# id_min = 0
# id_max = 0
# for id in cover_id_list:
# if id < id_min:
# id_min = id
# if id > id_max:
# id_max = id
# set_threshold = [id_max, id_min, threshold_of_ndvi_min, threshold_of_ndvi_max]
#
# CreateStadardXmlFile(xml_path, self.alg_xml_path, par_dict, set_threshold, model_xml_path).create_standard_xml()
#
SrcImagePath = self.__input_paras["SAR"]['ParaValue']
paths = SrcImagePath.split(';')
SrcImageName=os.path.split(paths[0])[1].split('.tar.gz')[0]
if len(paths) >= 2:
for i in range(1, len(paths)):
SrcImageName = SrcImageName+";"+os.path.split(paths[i])[1].split('.tar.gz')[0]
meta_xml_path = self.__product_dic+EXE_NAME+"Product.meta.xml"
CreateMetafile(self.__processing_paras['META'], self.alg_xml_path, model_xml_path, meta_xml_path).process(SrcImageName)
# if len(paths) >= 2:
# for i in range(1, len(paths)):
# SrcImageName = SrcImageName+";"+os.path.split(paths[i])[1].split('.tar.gz')[0]
# meta_xml_path = self.__product_dic+EXE_NAME+"Product.meta.xml"
# CreateMetafile(self.__processing_paras['META'], self.alg_xml_path, model_xml_path, meta_xml_path).process(SrcImageName)
model_path = "./product.xml"
meta_xml_path = os.path.join(self.__workspace_processing_path, SrcImageName + "-lef.meta.xml")
para_dict = CreateMetaDict(image_path, self.__processing_paras['META'], self.__workspace_processing_path,
out_path1, out_path2).calu_nature()
para_dict.update({"imageinfo_ProductName": "叶面积指数"})
para_dict.update({"imageinfo_ProductIdentifier": "LeafAreaIndex"})
para_dict.update({"imageinfo_ProductLevel": "LEVEL3"})
para_dict.update({"ProductProductionInfo_BandSelection": "1"})
para_dict.update({"ProductProductionInfo_AuxiliaryDataDescription": "MeasuredData,NDVI,LandCover"})
para_dict.update({"MetaInfo_Unit": "None"}) # 设置单位
CreateProductXml(para_dict, model_path, meta_xml_path).create_standard_xml()
# 文件夹打包
file.make_targz(self.__out_para, self.__product_dic)

View File

@ -256,36 +256,6 @@ class BlockProcess:
end_y) + '.' + out_type)
# print(out_dir_images)
# cut_factor_row = int(np.ceil(image.shape[0] / out_size))
# cut_factor_clo = int(np.ceil(image.shape[1] / out_size))
# for i in range(cut_factor_row):
# for j in range(cut_factor_clo):
#
# if i == cut_factor_row - 1:
# i = image.shape[0] / out_size - 1
# else:
# pass
#
# if j == cut_factor_clo - 1:
# j = image.shape[1] / out_size - 1
# else:
# pass
#
# start_x = int(np.rint(i * out_size))
# start_y = int(np.rint(j * out_size))
# end_x = int(np.rint((i + 1) * out_size))
# end_y = int(np.rint((j + 1) * out_size))
# out_dir_images = os.path.join(out_dir, img_name + '_' + str(start_x) + '_' + str(end_x) + '_' + str(start_y) + '_' + str(
# end_y) + '.' + out_type)
# + '/' + img_name \
# + '_' + str(start_x) + '_' + str(end_x) + '_' + str(start_y) + '_' + str(
# end_y) + '.' + out_type
# temp_image = image[start_x:end_x, start_y:end_y]
# out_image = Image.fromarray(temp_data)
# out_image = Image.fromarray(temp_image)
# out_image.save(out_dir_images)
data = ImageHandler.get_data(each_dir)
if ImageHandler.get_bands(each_dir) > 1:
# temp_data = data[:,start_x:end_x, start_y:end_y]
@ -445,12 +415,12 @@ class BlockProcess:
# if __name__ == '__main__':
# bp = BlockProcess()
# # # cut
# data_dir = r"D:\BaiduNetdiskDownload\HF\cut"
# out_dir = r"D:\BaiduNetdiskDownload\HF\cut_out"
# data_dir = r"D:\micro\WorkSpace\LandCover\Temporary\processing\feature_tif\cut"
# out_dir = r"D:\micro\WorkSpace\LandCover\Temporary\processing\feature_tif\combine"
# file_type = ['tif']
# out_type = 'tif'
# cut_size = 512
#
# cut_size = 1024
# #
# bp.cut_new(data_dir, out_dir, file_type, out_type, cut_size)
# # # combine
# # data_dir=r"D:\Workspace\SoilMoisture\Temporary\test"

View File

@ -235,7 +235,7 @@ class MachineLeaning:
@staticmethod
def sel_optimal_feature_set(X_train, Y_train, threshold=0.01):
"""
筛选最优特征组合
筛选最优特征组合(极度随机树
"""
model = ExtraTreesClassifier()
max = np.max(Y_train)