567 lines
29 KiB
Python
567 lines
29 KiB
Python
|
# -*- 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'))
|
|||
|
|
|
|||
|
|
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')
|
|||
|
|
tar = '-' + cf.get('target')
|
|||
|
|
productLevel = cf.get('productLevel')
|
|||
|
|
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.__workspace_preprocessed_path)
|
|||
|
|
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", "Covering", '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["Covering"])
|
|||
|
|
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]
|
|||
|
|
roi.trans_cover2mask(cover_mask_path, self.__preprocessed_paras["Covering"], 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)
|
|||
|
|
print("模型初值:" + str(params_X0))
|
|||
|
|
print("训练得到的模型系数:" + 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)
|
|||
|
|
|
|||
|
|
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")
|
|||
|
|
|
|||
|
|
SrcImagePath = self.__input_paras["SAR"]['ParaValue']
|
|||
|
|
paths = SrcImagePath.split(';')
|
|||
|
|
SrcImageName=os.path.split(paths[0])[1].split('.tar.gz')[0]
|
|||
|
|
|
|||
|
|
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)
|