microproduct/vegetationHeight-L-SAR/VegetationHeightMain.py

935 lines
49 KiB
Python
Raw Permalink Blame History

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

"""
@Project microproduct
@File OnePlantHeightMain.PY
@Function :主函数
@Author LMM
@Date 2021/10/19 14:39
@Version 1.0.0
"""
from VegetationHeightPrePro import PreProcess as Pp # 此行放在下面会报错,最好放在上面
from tool.algorithm.algtools.MetaDataHandler import Calibration
from tool.algorithm.xml.AlgXmlHandle import ManageAlgXML, CheckSource
from tool.algorithm.polsarpro.AHVToPolsarpro import AHVToPolsarpro
# from AlgXmlHandle import ManageAlgXML, CheckSource # 导入xml文件读取与检查文件
from VegetationHeightAlg import PlantHeightAlg, ROIAlg # 此处的ROIAlg可以对多波段数据进行掩膜与tool文件中的ROIAlg不同
from VegetationHeightOrthoAlg import IndirectOrthorectification as Inor
from tool.algorithm.image.ImageHandle import ImageHandler
from tool.algorithm.algtools.logHandler import LogHandler
from AHVToPolSarPro import AHVToPolSarProS2
from tool.algorithm.xml.CreatMetafile import CreateMetafile
from VegetationHeightXmlInfo import CreateDict, CreateStadardXmlFile
from Orthorectification import fine_registration,Orthorectification_RD, header_file_read_angle
import numpy as np
import os
import shutil
import logging
import datetime
import sys
import tarfile
import glob
import gc
import math
import sys
DEBUG = False
EXE_NAME = 'VegetationHeight'
LogHandler.init_log_handler('run_log\\' + EXE_NAME)
logger = logging.getLogger("mylog")
env_str = os.path.split(os.path.realpath(__file__))[0]
os.environ['PROJ_LIB'] = env_str
class VegetationHeightMain:
"""
植被高度主函数
"""
def __init__(self, alg_xml_path):
self.alg_xml_path = alg_xml_path
self.imageHandler = ImageHandler()
self.AHVToPolSarProS2 = AHVToPolSarProS2()
self.PlantHeightAlg = PlantHeightAlg()
self.ROIAlg = ROIAlg()
self.__alg_xml_handler = ManageAlgXML(alg_xml_path)
self.__check_handler = CheckSource(self.__alg_xml_handler)
self.__workspace_path = None
self.__workspace_soimois_path = None
self.__workspace_tem_dir_path = None
self.__task_id = None
self.__input_paras = {}
# self.__output_paras = {}
self.__processing_paras = {}
self.__preprocessed_paras = {}
self.__ahv = []
self.preprocessed_path = None
self.dem_pro = None
self.dem_geo = None
self.out_phi_tif = None
self.out_phi = None
self.out_kz_tif = None
self.out_kz = None
self.tif_angle_mask_path = None
self.__out_para = None
self.__feature_name_list = []
# 参考影像路径
self.__ref_img_path = ''
# 宽/列数
self.__cols = 0
# 高/行数
self.__rows = 0
# 坐标系
self.__proj = ''
# 影像投影变换矩阵
self.__geo = [0, 0, 0, 0, 0, 0]
def check_source(self):
"""
检查算法相关的配置文件,图像,辅助文件是否齐全
"""
env_str = os.getcwd()
logger.info("sysdir: %s", env_str)
if self.__check_handler.check_alg_xml() is False:
return False
if self.__check_handler.check_run_env() is False:
return False
input_para_names = ["MasterSar", "AuxiliarySar", "DEM"]
if self.__check_handler.check_input_paras(input_para_names) is False:
return False
self.__workspace_path = self.__alg_xml_handler.get_workspace_path()
self.__workspace_tem_dir_path = os.path.join(self.__workspace_path, EXE_NAME, "Temporary") # 获取根文件夹地址
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.__output_paras = self.__alg_xml_handler.get_output_paras() # 获取输出文件夹中的数据名、类型、地址
self.__processing_paras = self.__init_processing_paras(self.__input_paras) # 输出{文件名:地址}
# 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', r"VegetationHeightProduct.tar.gz")
self.__alg_xml_handler.write_out_para("VegetationHeight", self.__out_para) #写入输出参数
logger.info('check_source success!')
logger.info('progress bar :10')
return True
def __init_processing_paras(self, names):
""""
判断是文件,则获取其完整路径
判断是数值,则输出其值
将获取的路径和值都添加到 processing_paras 中
"""
processing_paras = {}
for name in names:
para = self.__input_paras[name]
if para is None:
logger.error(name + "is None!")
return False
if para['ParaType'] == 'File':
para_path = para['ParaValue']
if para['DataType'] == 'tif':
processing_paras.update({name: para_path})
elif para['DataType'] == "hgt":
processing_paras.update({name: para_path})
elif para['DataType'] == "xml":
processing_paras.update({name: para_path})
elif para['DataType'] == "nc":
processing_paras.update({name: para_path})
elif para['DataType'] == 'tar.gz':
tar_gz_dic = self.__dec_tar_gz(para_path, self.__workspace_unpack_path, para['ParaName'])
processing_paras.update(tar_gz_dic)
elif para['ParaType'] == 'Value':
if para['DataType'] == 'float':
value = float(para['ParaValue'])
elif para['DataType'] == 'int':
value = int(para['ParaValue'])
else: # 默认string
value = para['ParaValue']
processing_paras.update({name: value})
return processing_paras
def __dec_tar_gz(self, tar_gz_path, out_dir, para_name):
"""
输出以键值对{名称:路径}表示压缩包中头文件、全极化
:param tar_gz_path:
:param out_dir:
:param para_name:
:return:
"""
if os.path.exists(tar_gz_path) is False:
raise Exception('can not find file:', tar_gz_path)
# 创建文件夹
name = os.path.split(tar_gz_path)[1].split('.tar.gz')[0]
# name = os.path.split(tar_gz_path)[1].rstrip('.tar.gz')
file_dir = os.path.join(out_dir, name + '\\')
if os.path.exists(file_dir) is False:
os.makedirs(file_dir)
# 解压
t = tarfile.open(tar_gz_path)
t.extractall(path=file_dir)
para_dic = {}
hh_flag, hv_flag, vh_flag, vv_flag, angle_flag, ori_sim_flag = 0, 0, 0, 0, 0,0
# 获取头文件
if os.path.exists(file_dir + name + '\\'):
in_meta_paths = list(glob.glob(os.path.join(file_dir + name, '*.meta.xml')))
in_incidence_paths = list(glob.glob(os.path.join(file_dir + name, '*.incidence.xml')))
else:
in_meta_paths = list(glob.glob(os.path.join(file_dir, '*.meta.xml')))
in_incidence_paths =list(glob.glob(os.path.join(file_dir, '*.incidence.xml')))
if in_meta_paths is None or in_meta_paths == [] or in_incidence_paths is None or in_incidence_paths == [] :
raise Exception('.meta.xml or .incidence.xml data is missing in the tar_gz')
else:
name01 = "meta_" + para_name
para_dic.update({name01: in_meta_paths})
name02 = "meta_folder_" + para_name
para_dic.update({name02: file_dir}) # 记录解压的文件夹
name03 = "ladar_angle_" + para_name
para_dic.update({name03: in_incidence_paths})
self.image_meta_xml = in_meta_paths
# 获取文件夹内的文件
if os.path.exists(file_dir + name + '\\'):
in_tif_paths = list(glob.glob(os.path.join(file_dir + name + '\\', '*.tif')))
for temp in list(glob.glob(os.path.join(file_dir + name + '\\', '*.tiff'))):
in_tif_paths.append(temp)
else:
in_tif_paths = list(glob.glob(os.path.join(file_dir, '*.tif')))
for temp in list(glob.glob(os.path.join(file_dir, '*.tiff'))):
in_tif_paths.append(temp)
for in_tif_path in in_tif_paths:
# 获取极化类型
if 'hh' in os.path.basename(in_tif_path) or 'HH' in os.path.basename(in_tif_path):
name = 'HH_'+para_name
para_dic.update({name: in_tif_path})
hh_flag = 1
elif 'hv' in os.path.basename(in_tif_path) or 'HV' in os.path.basename(in_tif_path):
name = 'HV_' + para_name
para_dic.update({name: in_tif_path})
hv_flag = 1
elif 'vh' in os.path.basename(in_tif_path) or 'VH' in os.path.basename(in_tif_path):
name = 'VH_' + para_name
para_dic.update({name: in_tif_path})
vh_flag = 1
elif 'vv' in os.path.basename(in_tif_path) or 'VV' in os.path.basename(in_tif_path):
name = 'VV_' + para_name
para_dic.update({name: in_tif_path})
vv_flag = 1
elif 'IncidenceAngle' in os.path.basename(in_tif_path) or 'orth_sar_Angle' in os.path.basename(in_tif_path)\
or "local_incidence" in os.path.basename(in_tif_path) or "orth_sar_incidence" in os.path.basename(in_tif_path):
name = 'IncidenceAngle_' + para_name
para_dic.update({name: in_tif_path})
angle_flag = 1
elif 'ori_sim' in os.path.basename(in_tif_path) or 'Ori_sim' in os.path.basename(in_tif_path):
name = 'ori_sim_' + para_name
para_dic.update({name: in_tif_path})
ori_sim_flag = 1
else:
continue
if hh_flag != 1 or hv_flag != 1 or vh_flag != 1 or vv_flag != 1 or angle_flag != 1 or ori_sim_flag != 1 :
raise Exception('can not found files: HH、HV、VH、VV、IncidenceAngle、ori_sim in path:', tar_gz_path)
return para_dic
def __create_work_space(self):
"""
删除原有工作区文件夹,创建新工作区文件夹
"""
self.__workspace_preprocessing_path = self.__workspace_tem_dir_path + r"\preprocessing""\\"
self.__workspace_preprocessed_path = self.__workspace_tem_dir_path + r"\preprocessed""\\"
self.__workspace_preprocessed2_path = self.__workspace_tem_dir_path + r"\preprocessed2""\\"
self.__workspace_cut_path = self.__workspace_tem_dir_path + r"\cut""\\"
self.__workspace_Output_path = self.__workspace_path + r"VegetationHeight\Output""\\"
self.__workspace_master_slc_path = self.__workspace_tem_dir_path + r"\master_slc""\\"
self.__workspace_slave_slc_path = self.__workspace_tem_dir_path + r"\slave_slc""\\"
self.__workspace_unpack_path = self.__workspace_tem_dir_path + r"\unpack""\\"
self.__workspace_kz_phi_path = self.__workspace_tem_dir_path + r"\kz_phi""\\"
self.__workspace_roi_path = self.__workspace_tem_dir_path + r"\roi""\\"
path_list = [self.__workspace_preprocessing_path, self.__workspace_preprocessed_path,self.__workspace_preprocessed2_path,
self.__workspace_unpack_path, self.__workspace_cut_path, self.__workspace_master_slc_path,
self.__workspace_slave_slc_path, self.__workspace_Output_path, self.__workspace_kz_phi_path,
self.__workspace_roi_path]
for path in path_list:
if os.path.exists(path):
if DEBUG is True:
continue
self.del_folder(path)
os.makedirs(path)
else:
os.makedirs(path)
logger.info('create new workspace success!')
@staticmethod
def force_del_file(file_path):
"""
强制删除文件
"""
if os.path.isdir(file_path):
for main_dir, subdir, file_name_list in os.walk(file_path):
for filename in file_name_list:
apath = main_dir + filename
# noinspection PyBroadException
try:
os.remove(apath)
except Exception as error: # 使用windows命令行强制删除
os.system("del /f /q %s" % apath)
elif os.path.isfile(file_path) is True:
# noinspection PyBroadException
try:
os.remove(file_path)
except Exception as error: # 使用windows命令行强制删除
os.system("del /f /q %s" % file_path)
def del_file(self, path_data):
"""
只删除文件,不删除文件夹
"""
for i in os.listdir(path_data): # os.listdir(path_data)#返回一个列表,里面是当前目录下面的所有东西的相对路径
file_data = path_data + "\\" + i # 当前文件夹的下面的所有东西的绝对路径
if os.path.isfile(file_data) is True: # os.path.isfile判断是否为文件,如果是文件,就删除.如果是文件夹.递归给del_file.
os.remove(file_data)
return True
else:
self.del_file(file_data)
@staticmethod
def del_folder(path_data):
"""
删除整个文件夹
"""
if os.path.isdir(path_data):
shutil.rmtree(path_data)
@staticmethod
def del_folder_list(path_data_list):
"""
删除文件夹列表
"""
for path_data in path_data_list:
if os.path.isdir(path_data):
shutil.rmtree(path_data)
def del_temp_workspace(self):
"""
临时工作区
"""
if DEBUG is True:
return
path = self.__workspace_path + EXE_NAME + r'\Temporary'
if os.path.exists(path):
self.del_folder(path)
def preprocess_handle(self):
"""
预处理:检查坐标系、裁剪图像、重采样
"""
para_names = ['HH_MasterSar', 'HV_MasterSar', 'VH_MasterSar', 'VV_MasterSar', 'HH_AuxiliarySar',
'HV_AuxiliarySar', 'VH_AuxiliarySar', 'VV_AuxiliarySar', 'IncidenceAngle_MasterSar',
'IncidenceAngle_AuxiliarySar', 'DEM',"ori_sim_MasterSar","ori_sim_AuxiliarySar"]
# 检查图像坐标系
self.check_img_projection(para_names)
logger.info('progress bar :20')
try:
intersect_shp_path = self.__workspace_preprocessing_path + "IntersectPolygon.shp"
# 计算相交区域,将shp保存在 preprocessing 中
para_names_01 = ['VH_MasterSar', 'VH_AuxiliarySar','HH_MasterSar', 'HH_AuxiliarySar','HV_MasterSar',
'HV_AuxiliarySar','VV_MasterSar', 'VV_AuxiliarySar', 'IncidenceAngle_MasterSar',
'IncidenceAngle_AuxiliarySar', 'DEM', "ori_sim_MasterSar","ori_sim_AuxiliarySar"]
if self.cal_intersect_shp(intersect_shp_path, para_names_01) is False:
raise Exception('create intersect shp fail!')
logger.info('create intersect shp success!')
except:
intersect_shp_path=None
pass
# 裁剪图像存在cut文件夹中
cutted_img_path = self.cut_img(para_names, intersect_shp_path)
if cutted_img_path is None:
raise Exception('cut images success!')
logger.info('cut images success!')
self.del_folder(self.__workspace_preprocessing_path)
logger.info('progress bar :30')
# 重采样:重采样到跟主影像(MasterSar)的HH极化相同的分辨率然后保存到临时目录
self.preprocessed_path = self.resampling_img(para_names, cutted_img_path, cutted_img_path["HH_MasterSar"])
# 清除预处理缓存文件
logger.info('resampling imgs sucess!')
logger.info('preprocess_handle finished!')
self.del_folder(self.__workspace_cut_path)
logger.info('progress bar :40')
# 使用掩膜裁剪影像
try:
self.tif_angle_mask_path = self.create_roi()
except:
self.tif_angle_mask_path=None
pass
shutil.copy(self.preprocessed_path['HH_MasterSar'],os.path.join(self.__workspace_tem_dir_path,"HH_MasterSar_preprocessed.tif"))
logger.info('progress bar :50')
return True
def check_img_projection(self, para_names):
"""
读取每一张图像,检查图像坐标系,如果是投影坐标,则转换为地理坐标
:param para_names:需要检查的参数名称
"""
if len(para_names) == 0:
return False
for name in para_names:
proj = self.imageHandler.get_projection(self.__processing_paras[name])
keyword = proj.split("[", 2)[0]
if keyword == "PROJCS":
# 地理坐标系
para_dir = os.path.split(self.__processing_paras[name])
out_para = self.__workspace_preprocessing_path + \
para_dir[1].split(".", 1)[0] + "_EPSG32649.tif"
Pp.trans_projcs2geogcs(out_para, self.__processing_paras[name]) # 投影转地理
self.__processing_paras[name] = out_para
pass
elif keyword == "GEOGCS":
# 投影坐标系
pass
def cal_intersect_shp(self, shp_path, para_names):
"""
:param shp_path:相交区域矢量文件保存区域
:param para_names:判断相交影像的名称
:return: True or False
"""
if len(para_names) == 0:
return False
scopes = ()
for name in para_names:
scope_tuple = (self.imageHandler.get_scope(self.__processing_paras[name]),)
scopes += scope_tuple
intersect_polygon = Pp.intersect_polygon(scopes)
if intersect_polygon is None:
logger.error('image range does not overlap!')
return False
if Pp.write_polygon_shp(shp_path, intersect_polygon) is False:
return False
return True
def cut_img(self, para_names, shp_path):
"""
使用矢量数据裁剪影像
:param para_names:需要检查的参数名称
:param shp_path裁剪的shp文件
"""
if len(para_names) == 0:
return None
cutted_img_path = {}
# noinspection PyBroadException
try:
if not shp_path is None:
for name in para_names:
input_path = self.__processing_paras[name]
output_path = self.__workspace_cut_path + name + "_cut.tif"
Pp.cut_img(output_path, input_path, shp_path)
cutted_img_path.update({name: output_path})
logger.info("cut %s success!", name)
else:
for name in para_names:
input_path = self.__processing_paras[name]
output_path = self.__workspace_cut_path + name + "_cut.tif"
shutil.copy(input_path,output_path)
cutted_img_path.update({name: output_path})
logger.info("cut %s success!", name)
except BaseException:
try:
for name in para_names:
input_path = self.__processing_paras[name]
output_path = self.__workspace_cut_path + name + "_cut.tif"
shutil.copy(input_path,output_path)
cutted_img_path.update({name: output_path})
logger.info("cut %s success!", name)
except BaseException:
logger.error('cut_img failed!')
logger.error('cut_img failed!')
return None
return cutted_img_path
def resampling_img(self, para_names, img_paths, refer_img_path):
"""
以主影像为参考,对影像重采样
:param para_names:需要检查的参数名称
:param img_paths待重采样影像路径
:param refer_img_path参考影像路径
"""
if len(para_names) == 0 or len(img_paths) == 0:
return
prepro_imgs_path = {}
for name in para_names:
img_path = img_paths[name]
output_para = self.__workspace_preprocessed_path + name + "_preprocessed.tif"
Pp.resampling_by_scale(img_path, output_para, refer_img_path)
prepro_imgs_path.update({name: output_para})
logger.info("resampling %s success!", name)
return prepro_imgs_path
def check_LocalIncidenceAngle(self):
"""
将角度的无效值设置为nan把角度值转为弧度值
"""
# 主影像雷达入射角
in_tif_path_01 = self.preprocessed_path['IncidenceAngle_MasterSar']
angle_01 = self.imageHandler.get_data(in_tif_path_01)
if np.mean(angle_01) > np.pi:
angle_01 = np.deg2rad(angle_01) # 角度转弧度
angle_01[np.where(angle_01 >= 0.5 * np.pi)] = np.nan
angle_01[np.where(angle_01 < 0)] = np.nan
out_tif_01_path = os.path.splitext(in_tif_path_01)[0] + r"_check.tif"
self.imageHandler.write_img(out_tif_01_path, self.__proj, self.__geo, angle_01)
if os.path.isfile(in_tif_path_01) is True:
os.remove(in_tif_path_01)
self.preprocessed_path['IncidenceAngle_MasterSar'] = out_tif_01_path
# 辅影像雷达入射角
in_tif_path_02 = self.preprocessed_path['IncidenceAngle_AuxiliarySar']
angle_02 = self.imageHandler.get_data(in_tif_path_02)
if np.mean(angle_02) > np.pi:
angle_02 = np.deg2rad(angle_02) # 角度转弧度
angle_02[np.where(angle_02 >= 0.5 * np.pi)] = np.nan
angle_02[np.where(angle_02 < 0)] = np.nan
out_tif_path_02 = os.path.splitext(in_tif_path_02)[0] + r"_check.tif"
self.imageHandler.write_img(out_tif_path_02, self.__proj, self.__geo, angle_02)
if os.path.isfile(in_tif_path_02) is True:
os.remove(in_tif_path_02)
self.preprocessed_path['IncidenceAngle_AuxiliarySar'] = out_tif_path_02
def create_roi(self):
"""
计算 ROI
:return:
"""
# 利用影像的有效范围生成MASK
# tif_mask_0_path = self.__workspace_roi_path + "tif_mask_0.tif"
# ROIAlg.trans_tif2mask(tif_mask_0_path, self.preprocessed_path["HH_MasterSar"], 0, 0) # 掩膜影像中的0值
tif_mask_nan_path = self.__workspace_roi_path + "tif_mask_nan.tif"
ROIAlg.trans_tif2mask(tif_mask_nan_path, self.preprocessed_path["HV_MasterSar"], np.nan) # 掩膜影像中的空值
tif_auxil_nan_path = self.__workspace_roi_path + "tif_mask_nan.tif"
ROIAlg.trans_tif2mask(tif_auxil_nan_path, self.preprocessed_path["HV_AuxiliarySar"], np.nan) # 掩膜影像中的空值
self.check_LocalIncidenceAngle()
# 利用雷达入射角的有效范围生成MASK
incidence_angle_mask_path = self.__workspace_roi_path + "incidence_angle_mask.tif"
ROIAlg.trans_tif2mask(incidence_angle_mask_path, self.preprocessed_path["IncidenceAngle_MasterSar"], np.nan)
incidence_angle_Auxiliary_path = self.__workspace_roi_path + "incidence_angle_auxiliary.tif"
ROIAlg.trans_tif2mask(incidence_angle_Auxiliary_path, self.preprocessed_path["IncidenceAngle_AuxiliarySar"], np.nan)
# 利用行列号的范围生成MASK
ori_sim_MasterSar_path= self.__workspace_roi_path + "ori_sim_MasterSar_path.tif"
ROIAlg.trans_tif2mask(ori_sim_MasterSar_path, self.preprocessed_path["ori_sim_MasterSar"], 0, 99999999)
ori_sim_AuxiliarySar_path= self.__workspace_roi_path + "ori_sim_AuxiliarySar_path.tif"
ROIAlg.trans_tif2mask(ori_sim_AuxiliarySar_path, self.preprocessed_path["ori_sim_AuxiliarySar"], 0, 99999999)
tif_angle_mask_path = self.__workspace_roi_path + "tif_angle_mask.tif"
# ROIAlg.combine_mask(tif_angle_mask_path, tif_mask_0_path, incidence_angle_mask_path)
ROIAlg.combine_mask(tif_angle_mask_path, tif_mask_nan_path, incidence_angle_mask_path)
ROIAlg.combine_mask(tif_angle_mask_path, tif_angle_mask_path, incidence_angle_Auxiliary_path)
ROIAlg.combine_mask(tif_angle_mask_path, tif_angle_mask_path, ori_sim_MasterSar_path)
ROIAlg.combine_mask(tif_angle_mask_path, tif_angle_mask_path, ori_sim_AuxiliarySar_path)
ROIAlg.combine_mask(tif_angle_mask_path, tif_angle_mask_path, tif_auxil_nan_path)
logger.info('combine_mask success!')
logger.info('create ROI image success!')
# 计算roi区域
# for maindir, subdir, file_name_list in os.walk(self.__workspace_preprocessed_path):
# for filename in file_name_list:
# apath = maindir + filename
# ROIAlg.cal_roi(apath, apath, tif_angle_mask_path, background_value=1)
# logger.info('ROI cai image success!')
return tif_angle_mask_path
@staticmethod
def make_targz(output_filename, source_dir):
"""
一次性打包整个根目录。空子目录会被打包。
如果只打包不压缩,将"w:gz"参数改为"w:""w"即可。
:param output_filename:输出压缩包的完整路径eg:'E:\test.tar.gz'
:param source_dir:需要打包的跟目录eg: 'E:\testFfile\'打包文件夹里面的所有文件,'E:\testFfile'打包文件夹
"""
with tarfile.open(output_filename, "w:gz") as tar:
tar.add(source_dir, arcname=os.path.basename(source_dir))
def calc_fe_kz(self):
current_path = os.path.dirname(os.path.abspath(__file__))
head_config = os.path.join(current_path, "config.yaml")
master_meta_path_01 = self.__workspace_unpack_path + "MasterSar" # 主影像文件夹\
master_info= Inor(head_config).parasInfomation(master_meta_path_01)
Auxiliary_meta_path_02 = self.__workspace_unpack_path + "AuxiliarySar"
aux_info= Inor(head_config).parasInfomation(Auxiliary_meta_path_02)
FE_kx_exe= os.path.join(current_path,"PolSARproSim_FE_Kz.exe")
"""
PolSARproSim_FE_Kz out_dir Nlig Ncol dy (m) f0 (GHz) teta0 (deg) H0 (m) Bh (m) Bv (m)
Nlig 行数 Ncol 列数 dy 列分辨率(距离向分辨率) f0 频率 tata0 入射角 H0 卫星高度 Bh B_h Bv B_v
"""
"""PolSARproSim_FE_Kz.exe D:\forestHeight\Result\flat_remove 6764 5809 4.73307896 5.404999242769673 39.9528122 795449.5 102.36594559438615 -9.164148740434687"""
out_dir=self.__workspace_kz_phi_path
Nlig=master_info['ImageInformation']['height']
Ncol=master_info['ImageInformation']['width']
dy= master_info['ImageInformation']['ImageWidthSpace']
f0=1/master_info['ImageInformation']['lambda']
teta0=(master_info['ImageInformation']['incidenceAngle']["NearRange"]+master_info['ImageInformation']['incidenceAngle']["FarRange"])/2
H0=master_info['plate']['height']
# Bh=Bcosa
# Bv=Bsina
B=(master_info['plate']["position"][1]-aux_info['plate']["position"][1])**2+(master_info['plate']["position"][1]-aux_info['plate']["position"][2])**2+(master_info['plate']["position"][1]-aux_info['plate']["position"][3])**2
B=math.sqrt(B)
Bv=aux_info['plate']["height"]-master_info['plate']["height"]
Bh=math.sqrt(B**2-Bv**2)
cmd="{} {} {} {} {} {} {} {} {} {}".format(FE_kx_exe,out_dir,Nlig ,Ncol, dy ,f0, teta0, H0 ,Bh,Bv)
logger.info("calc kz {} :{} ".format(os.system(cmd),cmd))
self.out_kz = self.__workspace_kz_phi_path + "kz.bin" # degree rad
self.out_phi = self.__workspace_kz_phi_path + "flat_earth.bin" # degree rad
return master_info['ImageInformation']['lambda'],master_info['ImageInformation']['height'] , master_info['ImageInformation']['width'] ,aux_info['ImageInformation']['lambda'],aux_info['ImageInformation']['height'] , aux_info['ImageInformation']['width']
return True
def calc_phi(self, row, col):
"""计算phi参数"""
Mas_ori_sim = self.preprocessed_path["ori_sim_MasterSar"] # ori_sim文件路径
dem_resampled_path=Mas_ori_sim
self.dem_pro = ImageHandler().get_projection(dem_resampled_path) # 这里以dem数据的pro、geo作为后续输出tif的pro、geo
self.dem_geo = ImageHandler().get_geotransform(dem_resampled_path)
# #
try:
head_config = os.path.join(os.getcwd(), "config.ymal")
master_meta_path_01 = self.__workspace_unpack_path + "MasterSar" # 主影像文件夹
dem_resampled_path = self.preprocessed_path['DEM']
r1, lamda_01,height_01, width_01 = Inor(head_config).IndirectOrthorectification(master_meta_path_01, dem_resampled_path)
Auxiliary_meta_path_02 = self.__workspace_unpack_path + "AuxiliarySar" # 辅影像文件夹
r2, lamda_02,height_02, width_02= Inor(head_config).IndirectOrthorectification(Auxiliary_meta_path_02, dem_resampled_path)
phi_land = 4 * math.pi * (r1 - r2) * lamda_01
phi_land[np.where(phi_land is np.nan)] = 0
self.out_phi_tif = self.__workspace_kz_phi_path + "phi.tif"
ImageHandler().write_img(self.out_phi_tif, self.dem_pro, self.dem_geo, phi_land)
# no_file1=np.ones((int(row/6),int(col/2)),dtype=float)
no_file1 = np.ones((int(row / 1), int(col / 1)), dtype=float)
no_file_file1 = self.__workspace_kz_phi_path +"no_file_file1.tif"
ImageHandler().write_img(no_file_file1, self.dem_pro, self.dem_geo, no_file1, no_data='null')
array=Pp.resampling_array(self.out_phi_tif, no_file_file1)
self.out_phi = self.__workspace_kz_phi_path + "phi.bin"
self.AHVToPolSarProS2.array_bin(self.out_phi, array)
logger.info("calc phi finish ")
logger.info("phi write into tif finish ")
return lamda_01,height_01, width_01,lamda_02,height_02, width_02
except:
return False
pass
def calc_kz(self,row,col,lamda_01,height_01, width_01,lamda_02,height_02, width_02,multilook_row=1,multilook_col=1):
"""计算kz参数"""
# 方法2 通过正射计算的雷达入射角进行计算
try:
"""
#################################
head_config = os.path.join(os.getcwd(), "config.ymal")
master_meta_path_01 = self.__workspace_unpack_path + "MasterSar" # 主影像文件夹
dem_resampled_path = self.preprocessed_path['DEM']
r1, lamda_01,height_01, width_01 = Inor(head_config).IndirectOrthorectification(master_meta_path_01, dem_resampled_path)
Auxiliary_meta_path_02 = self.__workspace_unpack_path + "AuxiliarySar" # 辅影像文件夹
r2, lamda_02,height_02, width_02 = Inor(head_config).IndirectOrthorectification(Auxiliary_meta_path_02, dem_resampled_path)
# 2、将斜距数组写成tif文件
out_r1 = self.__workspace_kz_phi_path + "r1.tif"
out_r2 = self.__workspace_kz_phi_path + "r2.tif"
self.dem_pro = ImageHandler().get_projection(dem_resampled_path) # 这里以dem数据的pro、geo作为后续输出tif的pro、geo
self.dem_geo = ImageHandler().get_geotransform(dem_resampled_path)
ImageHandler().write_img(out_r1, self.dem_pro, self.dem_geo, r1)
ImageHandler().write_img(out_r2, self.dem_pro, self.dem_geo, r2)
master_incident_file = self.preprocessed_path['IncidenceAngle_MasterSar']
auxiliary_incident_file = self.preprocessed_path['IncidenceAngle_AuxiliarySar']
kz_array = Inor(head_config).calc_kz_array(master_incident_file, auxiliary_incident_file,
out_r1, out_r2, lamda_01)
##################################
"""
# 方法3 通过头文件的生成雷达入射角进行计算
###################################
m_incidence_file=self.__processing_paras["ladar_angle_MasterSar"][0] # _incidence.xml文件路径
a_incidence_file=self.__processing_paras["ladar_angle_AuxiliarySar"][0] # _incidence.xml文件路径
Mas_ori_sim = self.preprocessed_path["ori_sim_MasterSar"] # ori_sim文件路径
Aux_ori_sim = self.preprocessed_path["ori_sim_AuxiliarySar"]
m_all_angle_array=header_file_read_angle().read_all_angle(m_incidence_file,height_01,width_01)
m_angle_array = header_file_read_angle().get_orisim_angle(Mas_ori_sim, m_all_angle_array, height_01,width_01)
a_all_angle_array = header_file_read_angle().read_all_angle(a_incidence_file, height_02, width_02)
a_angle_array = header_file_read_angle().get_orisim_angle(Aux_ori_sim, a_all_angle_array, height_02, width_02)
# chazhi = m_angle_array - a_angle_array
v1 = m_angle_array / 180 * math.pi # 转成弧度
v2 = a_angle_array / 180 * math.pi
mean = (v1 + v2) / 2
kz_array = 4 * math.pi * (abs(v2 - v1)) / (np.sin(mean) * lamda_01)
self.mask_m_angle_array=m_angle_array
self.mask_a_angle_array = a_angle_array
####################################
self.out_kz_tif = self.__workspace_kz_phi_path + "kz.tif"
ImageHandler().write_img(self.out_kz_tif, self.dem_pro, self.dem_geo, kz_array)
if multilook_row!=1 and multilook_col!=1:
# no_file=np.ones((int(row/6),int(col/2)),dtype=float)
no_file = np.ones((int(row / multilook_row), int(col / multilook_col)), dtype=float)
no_file_file2 = self.__workspace_kz_phi_path +"no_file_file2.tif"
ImageHandler().write_img(no_file_file2, self.dem_pro, self.dem_geo, no_file, no_data='null')
array=Pp.resampling_array(self.out_kz_tif, no_file_file2)
array[np.where(array is np.nan)] = 0
self.out_kz = self.__workspace_kz_phi_path + "kz.bin"
self.AHVToPolSarProS2.array_bin(self.out_kz, array)
else:
self.out_kz = self.__workspace_kz_phi_path + "kz.bin"
self.AHVToPolSarProS2.array_bin(self.out_kz, kz_array)
logger.info("calc kz finish ")
except:
return False
def prod_handle(self):
"""
植被高度计算
"""
# 0 FE KZ
# 1、slc->s2
master_param = ['HH_MasterSar', 'HV_MasterSar', 'VH_MasterSar', 'VV_MasterSar']
auxiliary_param = ['HH_AuxiliarySar', 'HV_AuxiliarySar', 'VH_AuxiliarySar', 'VV_AuxiliarySar']
for para_m in master_param:
shutil.move(self.preprocessed_path[para_m], self.__workspace_master_slc_path)
for para_a in auxiliary_param:
shutil.move(self.preprocessed_path[para_a], self.__workspace_slave_slc_path)
atp = AHVToPolsarpro()
# 主影像->s2
master_s2 = self.__workspace_preprocessed2_path + "master_s2""\\"
os.makedirs(master_s2) # 输出文件夹
# self.AHVToPolSarProS2.api_ahv_to_polsarpro_s2(master_s2, self.__workspace_master_slc_path) # 全极化影像转S2矩阵
polarization = ['HH', 'HV', 'VH', 'VV']
calibration = Calibration.get_Calibration_coefficient(self.__processing_paras['meta_MasterSar'][0], polarization)
master_tif_path = atp.calibration(calibration, in_ahv_dir=self.__workspace_master_slc_path)
self.AHVToPolSarProS2.api_ahv_to_polsarpro_s2(master_s2, master_tif_path) # 全极化影像转S2矩阵
# 辅影像->s2
slave_s2 = self.__workspace_preprocessed2_path + "slave_s2""\\"
os.makedirs(slave_s2) # 输出文件夹
# self.AHVToPolSarProS2.api_ahv_to_polsarpro_s2(slave_s2, self.__workspace_slave_slc_path) # 全极化影像转S2矩阵
polarization = ['HH', 'HV', 'VH', 'VV']
calibration = Calibration.get_Calibration_coefficient(self.__processing_paras['meta_AuxiliarySar'][0], polarization)
aux_tif_path = atp.calibration(calibration, in_ahv_dir=self.__workspace_slave_slc_path)
self.AHVToPolSarProS2.api_ahv_to_polsarpro_s2(slave_s2, aux_tif_path) # 全极化影像转S2矩阵
logger.info("slc->s2 finish ")
try:
s11_bin = os.path.join(master_s2, 's11.bin')
arrayt11 = AHVToPolSarProS2().read_complex_bin_to_array(s11_bin)
rows, cols=arrayt11.shape[1],arrayt11.shape[2]
lamda_01,height_01, width_01,lamda_02,height_02, width_02=VegetationHeightMain.calc_phi(rows, cols)
except:
lamda_01,height_01, width_01,lamda_02,height_02, width_02=self.calc_fe_kz()
# self.force_del_file(self.__workspace_preprocessed_path)
# 2、去地平效应 s2->flat_earth_removal->s2
current_path = os.path.dirname(os.path.abspath(__file__))
# flat_earth_removal_master_slave_path = os.path.join(current_path, "flat_earth_removal_MasterSlave.exe")
flat_earth_removal_master_slave_path = os.path.join(current_path, "flat_earth_removal_Slave.exe")
master_fer =master_s2 ## self.__workspace_preprocessed2_path + "master_fer""\\"
#os.makedirs(master_fer)
slave_fer =slave_s2 #self.__workspace_preprocessed2_path + "slave_fer""\\"
#os.makedirs(slave_fer)
"""
# Soft/bin/data_process_dual/flat_earth_removal_Slave.exe -ids \x22$FlatEarthSlaveDirInput\x22
# -ods \x22$FlatEarthSlaveDirOutput\x22
# -fe \x22$FlatEarthFile\x22
# -cf $FlatEarthConjugate
# -fmt $FlatEarthFormat
# -ieee $FlatEarthIEEE
# -idf $DataFormatActive"
strcpy(UsageHelp,"\nflat_earth_removal_Slave.exe\n");
strcat(UsageHelp,"\nParameters:\n");
strcat(UsageHelp," (string) -ids input slave directory\n");
strcat(UsageHelp," (string) -ods output slave directory\n");
strcat(UsageHelp," (string) -fe input flat Earth file\n");
strcat(UsageHelp," (string) -fmt output format (cmplx / realdeg / realrad)\n");
strcat(UsageHelp," (int) -cf conjugate flag (1/0)\n");
strcat(UsageHelp," (int) -ieee ieee flag (1/0)\n");
strcat(UsageHelp," (string) -idf input data format (SPP, S2)\n");
strcat(UsageHelp,"\nOptional Parameters:\n");
strcat(UsageHelp," (noarg) -help displays this message\n");
"""
cmd="{} -ids {} -ods {} -fe {} -fmt cmplx -cf 0 -ieee 0 -idf s2".format(flat_earth_removal_master_slave_path,slave_s2,slave_fer,self.out_phi)
logger.info("remove plant flat {}: {} ".format(os.system(cmd),cmd))
#self.PlantHeightAlg.flat_earth_removal(master_s2, slave_s2, flat_earth_removal_master_slave_path,
# self.out_phi, master_fer, slave_fer, *(1, 0))
logger.info("flat earth removal finish")
# 3、s2+s2->t6
data_convert_s2_t6_path = os.path.join(current_path, "data_convert_MLK_S2_T6.exe")
master_slave_t6 = self.__workspace_preprocessed2_path + "master_slave_t6""\\"
os.makedirs(master_slave_t6)
self.PlantHeightAlg.s2_to_t6(master_fer, slave_fer, data_convert_s2_t6_path, master_slave_t6)
logger.info("s2 data convert t6 finish")
logger.info('progress bar :60')
t6_T11_bin = os.path.join(master_slave_t6, 'T11.bin')
coh_arrayt = AHVToPolSarProS2().read_none_complex_bin_to_array(t6_T11_bin)
rows, cols=coh_arrayt.shape[0],coh_arrayt.shape[1]
# 4、T6->boxcar_filter->T6
logger.info('start computing the filter...')
boxcar_filter_tool_path = os.path.join(current_path, "boxcar_filter_T6.exe")
# boxcar_filter_tool_path = os.path.join(current_path, "lee_refined_filter_T6.exe")
master_slave_t6_box = self.__workspace_preprocessed2_path + "master_slave_t6_box""\\"
os.makedirs(master_slave_t6_box)
self.PlantHeightAlg.polsar_boxcar_filter(master_slave_t6_box, boxcar_filter_tool_path,
master_slave_t6, *(3, 3, 0, 0))
# PlantHeightAlg().polsar_lee_filter(master_slave_t6_box, boxcar_filter_tool_path,
# master_slave_t6, *(3, 3, 0, 0))
logger.info("T6 lee_refined_filter finish")
logger.info('progress bar :85')
# 5、 T6->coherence_estimation->T6
coherence_estimation_path = os.path.join(current_path, "complex_coherence_estimation.exe")
# coherence_opt_estimation_path = os.path.join(current_path, "complex_coherence_opt_estimation_T6.exe")
master_slave_t6_box_coh = master_slave_t6_box
cor_list=["HH", "HV", "VV", "HHpVV", "HHmVV", "HVpVH", "LL", "LR", "RR", "HHVV"]
"""
complex_coherence_estimation.exe -id D:\forestHeight\Result\master_slave_t6_box -od D:\forestHeight\Result\master_slave_t6_box_coh1 -iodf T6 -type HHpVV -nwr 7 -nwc 7 -ofr 0 -ofc 0 -fnr 6065 -fnc 3716
"""
for cor_type in cor_list:
cmd="{} -id {} -od {} -iodf T6 -type {} -nwr 7 -nwc 7 -ofr 0 -ofc 0 -fnr {} -fnc {}".format(coherence_estimation_path,master_slave_t6_box,master_slave_t6_box_coh,cor_type,height_01, width_01)
logger.info("{} {}".format(os.system(cmd),cmd))
coherence_estimation_opt_path = os.path.join(current_path, "complex_coherence_opt_estimation.exe")
cmd="{} -id {} -od {} -iodf T6 -nwr 7 -nwc 7 -ofr 0 -ofc 0 -fnr {} -fnc {}".format(coherence_estimation_opt_path,master_slave_t6_box,master_slave_t6_box_coh,height_01, width_01)
logger.info("{} {}".format(os.system(cmd),cmd))
#self.PlantHeightAlg.complex_coherence_estimation_T6(master_slave_t6_box,
# coherence_estimation_path,
# master_slave_t6_box_coh, * (7, 7, 7, 7))
logger.info("T6->coherence estimation->T6 finish")
logger.info('progress bar :95')
#VegetationHeightMain.calc_kz(rows, cols,lamda_01,height_01, width_01,lamda_02,height_02, width_02, 1, 1) # 参数不为(1,1),须注释掉下面的掩膜部分
# 6、 T6->RVOG->产品
height_estimation_inversion_rvog_path = os.path.join(current_path,
"forest_height_estimation.exe")
"""
forest_height_estimation.exe -id D:\forestHeight\Result\master_slave_t6_box_coh1 -od D:\forestHeight\Result\master_slave_t6_box_coh_rvog -kz D:\forestHeight\Result\flat_remove\kz.bin -avg 0 -nr 6065 -nc 3716
"""
cmd="{} -id {} -od {} -kz {} -avg 0 -nr {} -nc {}".format(height_estimation_inversion_rvog_path,master_slave_t6_box_coh,master_slave_t6_box_coh,self.out_kz,height_01, width_01)
logger.info("{} {}".format(os.system(cmd),cmd))
#height_estimation_inversion_rvog_path = os.path.join(current_path,
# "height_estimation_inversion_procedure_RVOG.exe")
#file_gamma_high_file = master_slave_t6_box_coh + "cmplx_coh_HV.bin"
#file_gamma_low_file = master_slave_t6_box_coh + "cmplx_coh_HHmVV.bin"
# file_gamma_high_file = master_slave_t6_box_coh + "cmplx_coh_avg_HV.bin"
# file_gamma_low_file = master_slave_t6_box_coh + "cmplx_coh_avg_HHmVV.bin"
master_slave_t6_box_coh_rvog = self.__workspace_preprocessed2_path + "master_slave_t6_box_coh_rvog""\\"
os.makedirs(master_slave_t6_box_coh_rvog)
"""
forest_height_estimation.exe -id D:\forestHeight\Result\master_slave_t6_box_coh1 -od D:\forestHeight\Result\master_slave_t6_box_coh_rvog -kz D:\forestHeight\Result\flat_remove\kz.bin -avg 0 -nr 6065 -nc 3716
"""
shutil.copy(os.path.join(master_slave_t6_box_coh,"height.bin"),os.path.join(master_slave_t6_box_coh_rvog,"RVOG_heights.bin"))
#self.PlantHeightAlg.height_estimation_RVOG(master_slave_t6_box_coh, height_estimation_inversion_rvog_path,
# file_gamma_high_file, file_gamma_low_file,
# master_slave_t6_box_coh_rvog, self.out_kz, *(11, 0.4))
logger.info("T6->RVOG->产品 finish")
# 7、生成快视图并生成压缩包
vegetation_height_tif = self.__workspace_preprocessed2_path + "vegetation_height_tif""\\"
os.makedirs(vegetation_height_tif)
bin_path = master_slave_t6_box_coh_rvog + "RVOG_heights.bin"
array =np.fromfile(bin_path,dtype=np.float32).reshape(height_01, width_01) #AHVToPolSarProS2().read_none_complex_bin_to_array(bin_path)
# 对结果进行掩膜
#############
if self.tif_angle_mask_path is None:
array = array
else:
roi_img = self.imageHandler.get_band_array(self.tif_angle_mask_path)
array = array * roi_img # 掩膜
#self.mask_m_angle_array[np.where(self.mask_m_angle_array!=0)]=1
#self.mask_a_angle_array[np.where(self.mask_a_angle_array!=0)]=1
#self.mask_m_angle_array[np.where(self.mask_m_angle_array==0)]=math.nan
#self.mask_a_angle_array[np.where(self.mask_a_angle_array==0)]=math.nan
#array=array*self.mask_m_angle_array
#array=array*self.mask_a_angle_array
##############
self.dem_pro = ImageHandler().get_projection(self._VegetationHeightMain__input_paras["DEM"]['ParaValue']) # 这里以dem数据的pro、geo作为后续输出tif的pro、geo
self.dem_geo = ImageHandler().get_geotransform(self._VegetationHeightMain__input_paras["DEM"]['ParaValue'])
im_pro = self.dem_pro
im_geo = self.dem_geo
plant_height_product = vegetation_height_tif + "plant_heights_32649.tif"
self.imageHandler.write_img(plant_height_product, im_pro, im_geo, array, -999) # 生成tif
# self.force_del_file(plant_height)
# self.force_del_file(output_para)
# 快视图
self.imageHandler.write_quick_view(plant_height_product) # 生成快视图
xml_path = "./model_meta.xml"
tem_folder=self.__workspace_path + EXE_NAME + r"\Temporary""\\"
image_path= plant_height_product
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()
model_xml_path = os.path.join(tem_folder, "creat_standard.meta.xml") # 输出xml路径
CreateStadardXmlFile(xml_path, self.alg_xml_path, par_dict, model_xml_path).create_standard_xml()
SrcImageName = os.path.split(self.__input_paras["MasterSar"]['ParaValue'])[1].split('.tar.gz')[0]
# 生成压缩包
meta_xml_path = vegetation_height_tif+EXE_NAME+"Product.mate.xml"
CreateMetafile(self.image_meta_xml[0], self.alg_xml_path, model_xml_path, meta_xml_path).process(SrcImageName)
self.make_targz(self.__out_para, vegetation_height_tif)
logger.info("write quick view and .tar.gz finish")
self.del_folder(self.__workspace_master_slc_path)
self.del_folder(self.__workspace_slave_slc_path)
self.del_folder(self.__workspace_kz_phi_path)
self.del_folder(self.__workspace_preprocessed2_path)
self.del_folder(self.__workspace_roi_path)
logger.info('progress bar :100')
return True
if __name__ == "__main__":
start = datetime.datetime.now()
# noinspection PyBroadException
try:
if len(sys.argv) < 2:
xml_path = 'VegetationHeight.xml'
else:
xml_path = sys.argv[1]
VegetationHeightMain = VegetationHeightMain(xml_path)
if VegetationHeightMain.check_source() is False:
raise Exception('check_source() failed!')
if VegetationHeightMain.preprocess_handle() is False:
raise Exception('preprocess_handle() failed!')
#if VegetationHeightMain.calc_phi() is False:
# raise Exception('calc_kz_phi() failed!')
#if VegetationHeightMain.calc_kz() is False:
# raise Exception('calc_kz_phi() failed!')
if VegetationHeightMain.prod_handle() is False:
raise Exception('prod_handle() failed!')
logger.info('successful production of ' + EXE_NAME + ' products!')
except Exception:
logger.exception("run-time error!")
finally:
VegetationHeightMain.del_temp_workspace()
end = datetime.datetime.now()
msg = 'running use time: %s ' % (end - start)
logger.info(msg)