# -*- coding: UTF-8 -*- """ @Project ：MicroProduct @File ：LeafIndexMain.PY @Function ：主函数 @Author ：SHJ @Date ：2022/11/07 @Version ：2.0.0 修改历史： [修改序列] [修改日期] [修改者] [修改内容] 1 2022-6-27 李明明 1.增加配置文件config.ini; 2.修复快速图全黑的问题; 3.内部处理使用地理坐标系(4326)； """ from osgeo import gdalconst from tool.algorithm.algtools.PreProcess import PreProcess as pp # 此行放在下面会报错，最好放在上面 from tool.algorithm.xml.AlgXmlHandle import ManageAlgXML, CheckSource, InitPara # 导入xml文件读取与检查文件 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 from tool.algorithm.algtools.filter.lee_Filter import Filter from LeafIndexXmlInfo import CreateDict, CreateStadardXmlFile from tool.csv.csvHandle import csvHandle import os import shutil import logging import datetime import glob import numpy as np import scipy.spatial.transform # 用于解决打包错误 import scipy.spatial.transform._rotation_groups # 用于解决打包错误 import scipy.special.cython_special # 用于解决打包错误 import pyproj._compat from scipy.interpolate import griddata import sys import multiprocessing from tool.file.fileHandle import fileHandle from sample_process import read_sample_csv,combine_sample_attr,ReprojectImages2,read_tiff,check_sample,split_sample_list from tool.LAI.LAIProcess import train_WMCmodel,test_WMCModel,process_tiff cover_id_list = [] threshold_of_ndvi_min = 0 threshold_of_ndvi_max = 0 multiprocessing_num = int(cf.get('multiprocessing_num')) leaf_index_value_min = float(cf.get('product_value_min')) leaf_index_value_max = float(cf.get('product_value_max')) tar = r'-' + cf.get('tar') productLevel = cf.get('productLevel') FILTER_SISE = int(cf.get('filter_sise')) if cf.get('debug') == 'True': DEBUG = True else: DEBUG = False file = fileHandle(DEBUG) EXE_NAME = cf.get('exe_name') LogHandler.init_log_handler('run_log\\' + EXE_NAME) logger = logging.getLogger("mylog") csvh = csvHandle() # env_str = os.path.split(os.path.realpath(__file__))[0] env_str =os.path.dirname(os.path.abspath(sys.argv[0])) os.environ['PROJ_LIB'] = env_str class LeafIndexMain: """ 叶面积指数主函数 """ def __init__(self, alg_xml_path): self.alg_xml_path = alg_xml_path self.imageHandler = ImageHandler() self.BlockProcess = BlockProcess() self.__alg_xml_handler = ManageAlgXML(alg_xml_path) self.__check_handler = CheckSource(self.__alg_xml_handler) self.__workspace_path, self.__workspace_soimois_path,self.__workspace_soimois_path = None, None, None self.__input_paras, self.__processing_paras, self.__processing_paras = {}, {}, {} self.__out_para = None # 参考影像路径,坐标系 self.__ref_img_path, self.__proj = '', '' # 宽/列数,高/行数 self.__cols, self.__rows = 0, 0 # 影像投影变换矩阵 self.__geo = [0, 0, 0, 0, 0, 0] # 极化影像的名称 self.__sar_tif_name = None 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 self.__workspace_path = self.__alg_xml_handler.get_workspace_path() 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") SrcImageName = os.path.split(self.__input_paras["SAR"]['ParaValue'])[1].split('.tar.gz')[0] result_name = SrcImageName + tar +".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'] if Polarization == "HV": self.__sar_tif_name = "HV" elif Polarization == "HH": self.__sar_tif_name = "HH" elif Polarization == "VV": self.__sar_tif_name = "VV" elif Polarization == 'empty' or Polarization == 'Empty' or Polarization == 'EMPTY': self.__sar_tif_name = 'empty' else: raise ValueError("input Para 'Polarization' is not 'HV'、'HH'、'VV'or 'empty'!") self.__processing_paras = InitPara.init_processing_paras(self.__input_paras) self.__processing_paras.update(self.get_tar_gz_inf(self.__processing_paras["sar_path0"])) logger.info('check_source success!') return True def get_tar_gz_inf(self, tar_gz_path): para_dic = {} name = os.path.split(tar_gz_path)[1].rstrip('.tar.gz') file_dir = os.path.join(self.__workspace_preprocessing_path, name + '\\') file.de_targz(tar_gz_path, file_dir) # 元文件字典 para_dic.update(InitPara.get_meta_dic_new(InitPara.get_meta_paths(file_dir, name), name)) # tif路径字典 pol_dic = InitPara.get_polarization_mode(InitPara.get_tif_paths(file_dir, name)) flag_list = [0, 0] if self.__sar_tif_name == 'empty': for key, in_tif_path in pol_dic.items(): # 获取极化类型 if 'HV' == key and flag_list[0] == 0: para_dic.update({'HV': in_tif_path}) flag_list[0] = 1 elif 'HH' == key and flag_list[0] == 0: para_dic.update({'HH': in_tif_path}) flag_list[0] = 1 elif 'VV' == key and flag_list[0] == 0: para_dic.update({'VV': in_tif_path}) flag_list[0] = 1 elif 'LocalIncidenceAngle' == key: para_dic.update({'LocalIncidenceAngle': in_tif_path}) flag_list[1] = 1 else: continue else: for key, in_tif_path in pol_dic.items(): # 获取极化类型 if self.__sar_tif_name == key and flag_list[0] == 0: para_dic.update({self.__sar_tif_name: in_tif_path}) flag_list[0] = 1 elif 'LocalIncidenceAngle' == key: para_dic.update({'LocalIncidenceAngle': in_tif_path}) flag_list[1] = 1 if flag_list != [1, 1]: raise Exception('There are not ' + self.__sar_tif_name + '、LocalIncidenceAngle or meta.xml in path:', tar_gz_path) return para_dic def __create_work_space(self): """ 删除原有工作区文件夹,创建新工作区文件夹 """ self.__workspace_preprocessing_path = self.__workspace_path + EXE_NAME + "\\Temporary\\preprocessing""\\" self.__workspace_preprocessed_path = self.__workspace_path + EXE_NAME + "\\Temporary\\preprocessed""\\" self.__workspace_processing_path = self.__workspace_path + EXE_NAME + "\\Temporary\\processing""\\" self.__workspace_maskcai_image_path = self.__workspace_path + EXE_NAME + "\\Temporary\\maskcai_image""\\" self.__workspace_maskcai_SoilMoi_path = self.__workspace_path + EXE_NAME + "\\Temporary\\maskcai_SoilMoi""\\" self.__workspace_maskcai_localang_path = self.__workspace_path + EXE_NAME + "\\Temporary\\maskcai_localang\\" self.__workspace_cai_sartif_path = self.__workspace_path + EXE_NAME + "\\Temporary\\cai_sartif""\\" self.__workspace_cai_SoilMoi_path = self.__workspace_path + EXE_NAME + "\\Temporary\\cai_SoilMoi""\\" self.__workspace_cai_localang_path = self.__workspace_path + EXE_NAME + "\\Temporary\\cai_localang""\\" self.__workspace_Leafindex_path = self.__workspace_path + EXE_NAME + "\\Temporary\\Leafindex""\\" self.__product_dic = self.__workspace_processing_path + "product\\" path_list = [self.__workspace_preprocessing_path, self.__workspace_preprocessed_path, self.__workspace_processing_path, self.__workspace_maskcai_image_path, self.__workspace_maskcai_SoilMoi_path, self.__workspace_maskcai_localang_path, self.__workspace_cai_sartif_path, self.__workspace_cai_SoilMoi_path, self.__workspace_cai_localang_path, self.__workspace_Leafindex_path] file.creat_dirs(path_list) logger.info('create new workspace success!') def del_temp_workspace(self): """ 临时工作区 """ if DEBUG is True: return path = self.__workspace_path + EXE_NAME + r'\Temporary' if os.path.exists(path): file.del_folder(path) def preprocess_handle(self): """ 预处理 """ para_names = [self.__sar_tif_name, 'LocalIncidenceAngle', "NDVI", "surface_coverage", 'soilMeasured'] ref_img_name = self.__sar_tif_name p = pp() self.__preprocessed_paras = p.preprocessing(para_names, ref_img_name, self.__processing_paras, self.__workspace_preprocessing_path, self.__workspace_preprocessed_path) self.__ref_img_path, self.__cols, self.__rows, self.__proj, self.__geo = p.get_ref_inf( self.__preprocessed_paras[ref_img_name]) logger.info('preprocess_handle success!') logger.info('progress bar: 30%') def cal_roi(self): """ 创建掩膜，并裁剪成统一范围 """ # 1、利用叠掩区域生成局部入射角的Mask processing_path = self.__workspace_processing_path # 利用角度为nan生成Mask pp.check_LocalIncidenceAngle(self.__preprocessed_paras['LocalIncidenceAngle'], self.__preprocessed_paras['LocalIncidenceAngle']) angle_nan_mask_path = processing_path + 'angle_nan_mask.tif' roi.trans_tif2mask(angle_nan_mask_path, self.__preprocessed_paras['LocalIncidenceAngle'], np.nan) # 利用影像的有效范围生成MASK tif_nan_mask_path = processing_path + 'tif_mask_nan.tif' roi.trans_tif2mask(tif_nan_mask_path, self.__preprocessed_paras[self.__sar_tif_name], np.nan) # tif_zero_mask_path = processing_path + 'tif_mask_zero.tif' # roi.trans_tif2mask(tif_zero_mask_path, self.__preprocessed_paras[self.__sar_tif_name], 0,0) # 2、利用cover计算植被覆盖范围 cover_mask_path = self.__workspace_processing_path + "SurfaceCover_mask.tif" if self.__processing_paras['CoveringIDs'] == 'empty': cover_data = ImageHandler.get_data(self.__preprocessed_paras["surface_coverage"]) cover_data[np.where(np.isnan(cover_data))] = 0 cover_id_list = list(np.unique(cover_data)) else: cover_id_list = list(self.__processing_paras['CoveringIDs'].split(';')) cover_id_list = [int(num) for num in cover_id_list] # cover_id_list = list(self.__processing_paras['CoveringIDs'].split(';')) # cover_id_list = [int(num) for num in cover_id_list] roi.trans_cover2mask(cover_mask_path, self.__preprocessed_paras["surface_coverage"], cover_id_list) # 3、利用NDVI计算裸土范围该指数的输出值在 -1.0 和 1.0 之间，大部分表示植被量， # 负值主要根据云、水和雪而生成 # 接近零的值则主要根据岩石和裸土而生成。 # 较低的（小于等于 0.1）NDVI 值表示岩石、沙石或雪覆盖的贫瘠区域。 # 中等值（0.2 至 0.3）表示灌木丛和草地 # 较高的值（0.6 至 0.8）表示温带雨林和热带雨林。 ndvi_mask_path = self.__workspace_processing_path + "Ndvi_mask.tif" ndvi_scope = list(self.__processing_paras['NDVIScope'].split(';')) threshold_of_ndvi_min = float(ndvi_scope[0]) threshold_of_ndvi_max = float(ndvi_scope[1]) roi.trans_tif2mask(ndvi_mask_path, self.__preprocessed_paras["NDVI"], threshold_of_ndvi_min, threshold_of_ndvi_max) logger.info('create masks success!') # 4、利用覆盖范围和裸土范围 生成总的MASK bare_land_mask_path = self.__workspace_processing_path + "bare_land_mask.tif" roi.combine_mask(bare_land_mask_path, ndvi_mask_path, cover_mask_path) roi.combine_mask(bare_land_mask_path, bare_land_mask_path, tif_nan_mask_path) # roi.combine_mask(bare_land_mask_path, bare_land_mask_path, tif_zero_mask_path) roi.combine_mask(bare_land_mask_path, bare_land_mask_path, angle_nan_mask_path) logger.info('combine_mask success!') # 5、使用ROI区域进行掩膜 roi.cal_roi(self.__workspace_maskcai_localang_path + "LocalIncidenceAngle_preproed.tif", self.__preprocessed_paras["LocalIncidenceAngle"], bare_land_mask_path, background_value=1) logger.info('create ROI three image success!') return bare_land_mask_path def cal_leaf_index(self, vm_file, image_file, angle_file, v2_path, num): """ 叶面积指数计算公式 V1:植被层特征参数，V1设定为1 v2：植被层特征参数，V2为叶面积指数 A,B,C,D 经验系数 vm 土壤体积含水量 """ # 滤波 # image_filter_path = v2_path.replace("LeafAreaIndexProduct","lee_filer_img") # f.lee_filter(image_file, image_filter_path, FILTER_SISE) a = self.__processing_paras["A"] b = self.__processing_paras["B"] c = self.__processing_paras["C"] d = self.__processing_paras["D"] # k = self.__processing_paras["K"] # m = self.__processing_paras["M"] v1 = 1 vm_value = self.imageHandler.get_band_array(vm_file, 1) image_value = self.imageHandler.get_band_array(image_file, 1) # Filter.lee_filter_array(image_value, image_value, FILTER_SISE) angle = self.imageHandler.get_band_array(angle_file, 1) cos_angle_value = np.cos(angle) # 开始计算叶面积指数 back_soil = c * vm_value - d t = (image_value - a * v1 * cos_angle_value) / (back_soil - a * v1 * cos_angle_value) where_0 = np.where(t <= 0) t[where_0] = 1 array_v2 = -1 * (np.log(t)) * cos_angle_value / (2 * b) array_v2[where_0] = 0 self.imageHandler.write_img(v2_path, self.__proj, self.__geo, array_v2) logger.info('block:%s cal LeafIndex success!', num) def create_soil_moisture_tif(self): """ 创建土壤介电常数影像 """ # 读取实测值，从经纬度坐标系转为图像坐标系 measured_data_img = csvh.trans_measuredata(csvh.readcsv(self.__processing_paras['MeasuredData']), self.__ref_img_path) # 插值法生成土壤水分二维影像 grid_x, grid_y = np.mgrid[0:self.__rows:1, 0:self.__cols:1] points = np.array(measured_data_img)[:, 0:2] values = np.array(measured_data_img)[:, 2].T grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest') self.__processing_paras['SoilMoisture'] = os.path.join(self.__workspace_maskcai_SoilMoi_path, 'soil_moisture.tif') self.imageHandler.write_img(self.__processing_paras['SoilMoisture'], self.__proj, self.__geo, grid_z0) logger.info('create soil_moisture image success!') logger.info('progress bar: 40%') def cal_empirical_parameters(self): work_path = self.__workspace_path + EXE_NAME + "\\Temporary\\empirical""\\" # b. 结果工作 result_dir_path = self.__workspace_path + EXE_NAME + "\\Temporary\\empirical_result""\\" path_list = [work_path, result_dir_path] file.creat_dirs(path_list) # 1. 后向散射系数 dB sigma_path = self.__workspace_maskcai_image_path + self.__sar_tif_name+'.tif' # 2. 局地入射角 incident_angle_path = self.__workspace_maskcai_localang_path + "LocalIncidenceAngle_preproed.tif" # 3. 样本csv地址 lai_csv_path = self.__processing_paras['laiMeasuredData'] # 4. NDVI影像地址 -- 修正模型 NDVI_tiff_path = self.__preprocessed_paras["NDVI"] # 5. 土壤含水量影像地址 soil_water_tiff_path = self.__preprocessed_paras['soilMeasured'] # 6. 土壤含水量样本地址 soil_water_csv_path = r"" # 7. 选择土壤含水量影像 soil_water = 'tiff' # 8. 输出图片 train_err_image_path = os.path.join(result_dir_path, "train_image.png") NDVI_min = -1 # 完全裸土对应的 NDVI 值 NDVI_max = 1 # 完全植被覆盖对应的 NDVI 值 # 临时变量 soil_tiff_resample_path = os.path.join(work_path, "soil_water.tiff") # 与 后向散射系数同样分辨率的 土壤水分影像 NDVI_tiff_resample_path = os.path.join(work_path, 'NDVI.tiff') # 与 后向散射系数产品同样分辨率的 NDVI影像 incident_angle_resample_path = os.path.join(work_path, "localincangle.tiff") # 读取数据 lai_sample = read_sample_csv(lai_csv_path) # 读取样本数据 sigma_tiff = read_tiff(sigma_path) # 读取后向散射系数 incident_angle = read_tiff(incident_angle_path) # 读取局地入射角 # 对于土壤水分、NDVI做重采样 ReprojectImages2(soil_water_tiff_path, sigma_path, soil_tiff_resample_path, resampleAlg=gdalconst.GRA_Bilinear) ReprojectImages2(NDVI_tiff_path, sigma_path, NDVI_tiff_resample_path, resampleAlg=gdalconst.GRA_Bilinear) ReprojectImages2(incident_angle_path, sigma_path, incident_angle_resample_path, resampleAlg=gdalconst.GRA_NearestNeighbour) soil_water_tiff = read_tiff(soil_tiff_resample_path) # 读取土壤含水量影像 NDVI_tiff = read_tiff(NDVI_tiff_resample_path) # 引入NDVI incident_angle = read_tiff(incident_angle_resample_path) # 读取局地入射角 # 处理归一化植被指数 F_VEG = (NDVI_tiff['data'] - NDVI_min) / (NDVI_max - NDVI_min) # 处理得到植被覆盖度 soil_water_tiff['data'] = soil_water_tiff['data'] / 100.0 # 转换为百分比 incident_angle['data'] = incident_angle['data'] * np.pi / 180.0 # 转换为弧度值 sigma_tiff['data'] = np.power(10, (sigma_tiff['data'] / 10)) # 转换为线性值 # float32 转 float64 soil_water_tiff['data'] = soil_water_tiff['data'].astype(np.float64) incident_angle['data'] = incident_angle['data'].astype(np.float64) sigma_tiff['data'] = sigma_tiff['data'].astype(np.float64) # 将土壤水分与lai样本之间进行关联 lai_water_sample = [] # ['日期', '样方编号', '经度', '纬度', 'LAI','土壤含水量'] if soil_water == 'tiff': lai_water_sample = combine_sample_attr(lai_sample, soil_water_tiff) pass else: # 这个暂时没有考虑 pass # 将入射角、后向散射系数与lai样本之间进行关联 lai_water_inc_list = combine_sample_attr(lai_water_sample, incident_angle) # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角'] lai_waiter_inc_sigma_list = combine_sample_attr(lai_water_inc_list, sigma_tiff) # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数'] # lai_waiter_inc_sigma_NDVI_list=combine_sample_attr(lai_waiter_inc_sigma_list,NDVI_tiff) # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数','NDVI'] lai_waiter_inc_sigma_list = check_sample( lai_waiter_inc_sigma_list) # 清理样本 ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数'] # lai_waiter_inc_sigma_NDVI_list=check_sample(lai_waiter_inc_sigma_NDVI_list) # ['日期','样方编号','经度','纬度','叶面积指数',"后向散射系数",'土壤含水量','入射角','后向散射系数','NDVI'] # 数据集筛选 lai_waiter_inc_sigma_list_result = [] # 筛选保留的数据集 logger.info("保留得数据集如下") for i in range(len(lai_waiter_inc_sigma_list)): if i in []: continue logger.info(str(lai_waiter_inc_sigma_list[i])) lai_waiter_inc_sigma_list_result.append(lai_waiter_inc_sigma_list[i]) lai_waiter_inc_sigma_list = lai_waiter_inc_sigma_list_result # [sample_train,sample_test]=split_sample_list(lai_waiter_inc_sigma_list,0.6) # step 1 切分数据集 [sample_train, sample_test] = [lai_waiter_inc_sigma_list[:], lai_waiter_inc_sigma_list[:]] # step 1 切分数据集 logger.info("训练模型") a = self.__processing_paras["A"] b = self.__processing_paras["B"] c = self.__processing_paras["C"] d = self.__processing_paras["D"] params_X0 = [a, b, c, d, 0.771, -0.028] params_arr = train_WMCmodel(sample_train, params_X0, train_err_image_path, False) logging.info("模型初值:\t{}".format(str(params_X0))) logging.info("训练得到的模型系数:\t{}".format(str(params_arr))) self.__processing_paras.update({"A": params_arr[0]}) self.__processing_paras.update({"B": params_arr[1]}) self.__processing_paras.update({"C": params_arr[2]}) self.__processing_paras.update({"D": params_arr[3]}) def block_process(self,start): """ 生成叶面积指数产品 """ # 生成土壤水分影像 # self.create_soil_moisture_tif() if os.path.exists(self.__preprocessed_paras['soilMeasured']): soil_new = os.path.join(self.__workspace_maskcai_SoilMoi_path, 'soil_moisture.tif') shutil.copy(self.__preprocessed_paras['soilMeasured'], soil_new) lee_path = os.path.join(self.__workspace_preprocessed_path, self.__sar_tif_name + '.tif') Filter().lee_process_sar(self.__preprocessed_paras[self.__sar_tif_name], lee_path, 3, 0.25) shutil.copyfile(lee_path, self.__workspace_maskcai_image_path + self.__sar_tif_name+'.tif') # shutil.copyfile(self.__preprocessed_paras[self.__sar_tif_name], self.__workspace_maskcai_image_path + self.__sar_tif_name+'.tif') logger.info('progress bar: 50%') # 模型训练得到经验系数 self.cal_empirical_parameters() block_size = self.BlockProcess.get_block_size(self.__rows, self.__cols) 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!') # 计算每小块的叶面积指数 in_tif_paths = list(glob.glob(os.path.join(self.__workspace_cai_sartif_path, '*.tif'))) num = 0 # 开启多进程处理 processes_num = min([len(in_tif_paths), multiprocessing_num]) # Filter().lee_filter_multiprocess(in_tif_paths, in_tif_paths, FILTER_SISE, processes_num) pool = multiprocessing.Pool(processes=processes_num) pl = [] for name in in_tif_paths: suffix = '_' + name.split('_')[-4] + "_" + name.split('_')[-3] + "_" + name.split('_')[-2] + "_" + \ name.split('_')[-1] tif_path = name vm_path = self.__workspace_cai_SoilMoi_path + "soil_moisture" + suffix angle_path = self.__workspace_cai_localang_path + "LocalIncidenceAngle_preproed" + suffix v2_save_path = self.__workspace_Leafindex_path + "LeafAreaIndexProduct" + suffix pl.append(pool.apply_async(self.cal_leaf_index, (vm_path, tif_path, angle_path, v2_save_path, num))) logger.info('total:%s, block:%s calculating leaf_index!', len(in_tif_paths), num) num = num + 1 pool.close() pool.join() logger.info('all img cal LeafIndex success!') logger.info('progress bar: 80%') ref_tif_path = self.__ref_img_path date_type = 'float32' 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() SrcImageName = os.path.split(self.__input_paras["SAR"]['ParaValue'])[1].split('.tar.gz')[0] + tar + '.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!') proj, geos, data = self.imageHandler.read_img(product_path) data[data < leaf_index_value_min] = leaf_index_value_min data[data > leaf_index_value_max] = leaf_index_value_max self.imageHandler.write_img(product_path, proj, geos, data) # 生成快视图 self.imageHandler.write_quick_view(product_path) xml_path = "./model_meta.xml" tem_folder = self.__workspace_path + EXE_NAME + r"\Temporary""\\" 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() # 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) model_path = "./product.xml" meta_xml_path = os.path.join(self.__product_dic, SrcImageName + tar + ".meta.xml") para_dict = CreateMetaDict(image_path, self.__processing_paras['Origin_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": productLevel}) 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() 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) # 文件夹打包 file.make_targz(self.__out_para, self.__product_dic) logger.info('process_handle success!') logger.info('progress bar: 100%') def process_handle(self,start): self.cal_roi() self.block_process(start) return True if __name__ == '__main__': multiprocessing.freeze_support() # 解决多进程打包错误 start = datetime.datetime.now() try: if len(sys.argv) < 2: xml_path = 'LeafAreaIndex.xml' else: xml_path = sys.argv[1] main_handler = LeafIndexMain(xml_path) if main_handler.check_source() is False: raise Exception('check_source() failed!') if main_handler.preprocess_handle() is False: raise Exception('preprocess_handle() failed!') if main_handler.process_handle(start) is False: raise Exception('process_handle() failed!') logger.info('successful production of ' + EXE_NAME + ' products!') except Exception: logger.exception('run-time error!') finally: main_handler.del_temp_workspace() end = datetime.datetime.now() msg = 'running use time: %s ' % (end - start) logger.info(msg)