1080 lines
53 KiB
Python
1080 lines
53 KiB
Python
"""
|
||
@Project ：microproduct
|
||
@File ：AtmosphericDelayMain.py
|
||
@Function ：主程序，执行入口
|
||
@Author ：LMM
|
||
@Date ：2021/09/3 14:39
|
||
@Version ：1.0.0
|
||
"""
|
||
import pyproj._compat #解决打包错误
|
||
from osgeo import gdal
|
||
import glob
|
||
from pykrige import OrdinaryKriging # 强制引入包
|
||
from tool.algorithm.algtools.PreProcess import PreProcess as pp
|
||
from tool.algorithm.transforml1a.transHandle import TransImgL1A
|
||
from tool.algorithm.xml.AlgXmlHandle import ManageAlgXML, CheckSource # 导入xml文件读取与检查文件
|
||
from tool.algorithm.xml.CreatMetafile import CreateMetafile
|
||
from AtmosphericDelayAuxData import NcHandle,ReadImage
|
||
from AtmosphericDelayAlg import AtmosphericDelay as Ady
|
||
from tool.algorithm.algtools.ROIAlg import ROIAlg
|
||
from AtmosphericDelayXmlInfo import CreateDict, CreateStadardXmlFile
|
||
from AtmosphericDelayDT import DataTransform
|
||
# from AtmosphericDealyTool import tropo_EAR5
|
||
from ConfigeHandle import Config as cf
|
||
import os
|
||
import shutil
|
||
import datetime
|
||
import sys
|
||
import tarfile
|
||
from tool.algorithm.algtools.logHandler import LogHandler
|
||
import logging
|
||
import numpy as np
|
||
import netCDF4 as Nc #解决打包错误
|
||
import cftime #解决打包错误
|
||
from tool.algorithm.image.ImageHandle import ImageHandler
|
||
from autorun import auto_run_main
|
||
import gc
|
||
|
||
|
||
EXE_NAME = cf.get('exe_name')
|
||
|
||
if cf.get('debug') == 'True':
|
||
DEBUG = True
|
||
else:
|
||
DEBUG = False
|
||
|
||
LogHandler.init_log_handler(r'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 AtmosphericMain:
|
||
"""
|
||
大气延迟校正主函数
|
||
"""
|
||
def __init__(self, alg_xml_path):
|
||
self.alg_xml_path = alg_xml_path
|
||
self.imageHandler = ImageHandler()
|
||
self.NcHandle = NcHandle()
|
||
self.ReadImage = ReadImage()
|
||
self.ROIAlg = ROIAlg()
|
||
self.__alg_xml_handler = ManageAlgXML(alg_xml_path) # 导入大气延迟的xml文件
|
||
self.__check_handler = CheckSource(self.__alg_xml_handler)
|
||
|
||
self.__workspace_path = None
|
||
self.__workspace_tem_dir_path = None
|
||
self.__input_paras = {}
|
||
|
||
self.__processing_paras = {}
|
||
self.__preprocessed_paras = {}
|
||
self.__preprocessed_paras2 = {}
|
||
self.__out_para = None
|
||
self.pro_img_path = None
|
||
self.tif_angle_mask_path = None
|
||
self.intersect_shp_path = None
|
||
# 坐标系
|
||
self.__proj = ''
|
||
# 影像投影变换矩阵
|
||
self.__geo = [0, 0, 0, 0, 0, 0]
|
||
|
||
def check_source(self):
|
||
"""
|
||
检查算法相关的配置文件，图像，辅助文件是否齐全
|
||
"""
|
||
self.env_str = os.getcwd()
|
||
logger.info("sysdir: %s", self.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 = ["MasterSarData", "AuxiliarySarData", "MasterNC", "AuxiliaryNC", "DEM", "box"]
|
||
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 = self.__workspace_path + EXE_NAME+"\\Temporary""\\"
|
||
self.__create_work_space()
|
||
self.__input_paras = self.__alg_xml_handler.get_input_paras() # 获取输入文件夹中的数据名、类型、地址
|
||
self.__processing_paras = self.__init_processing_paras(self.__input_paras) # 输出{文件名:地址}
|
||
aux_path = os.path.join(self.__workspace_origin_path, "AuxDir")
|
||
if os.path.exists(aux_path) is False:
|
||
os.mkdir(aux_path)
|
||
self.__processing_paras.update({'AuxDir': aux_path})
|
||
self.__out_para = os.path.join(self.__workspace_path, EXE_NAME, 'Output', r"AtmosphericDelayProduct.tar.gz")
|
||
self.__alg_xml_handler.write_out_para("AtmosphericDelayProduct", self.__out_para) # 写入输出参数
|
||
|
||
# 判断气象数据是tif格式还是file格式
|
||
name_list = ["MasterNC", "AuxiliaryNC"]
|
||
for name in name_list:
|
||
if self.__input_paras[name]["DataType"] in ["file", "File"]:
|
||
self.NcHandle.copy_tif(self.__processing_paras[name], self.__workspace_preprocessing_path, name)
|
||
elif self.__input_paras[name]["DataType"] == "nc":
|
||
pass
|
||
# self.NcHandle.check_nc(self.__processing_paras[name]) # 检查气象数据
|
||
# self.get_prodataset(name) # 读取气象数据，另存为tif
|
||
|
||
logger.info('check_source success!')
|
||
logger.info('progress bar :10％')
|
||
return True
|
||
|
||
def __init_processing_paras(self, names):
|
||
""""
|
||
param: names:字典列表，每个字典为一个输入产品的配置信息
|
||
"""
|
||
processing_paras = {}
|
||
|
||
for name in names:
|
||
para = self.__input_paras[name]
|
||
if para is None:
|
||
logger.error(name + "is None!")
|
||
return False
|
||
para_path = para['ParaValue']
|
||
if para['ParaType'] == 'File':
|
||
if para['DataType'] == "nc":
|
||
processing_paras.update({name: para_path})
|
||
elif para['DataType'] == 'file':
|
||
if name in ["MasterNC", "AuxiliaryNC"]:
|
||
processing_paras.update({name: para_path})
|
||
|
||
if name in ['MasterSarData', 'AuxiliarySarData']:
|
||
paths = para['ParaValue'].split(';')
|
||
|
||
# 获取文件名中的关键字:如"20190206"/"20190113"
|
||
if name == 'MasterSarData':
|
||
filename = os.path.basename(paths[0])
|
||
part = filename.split("_")
|
||
name = part[5]
|
||
self.mas_key_word = name[0:8]
|
||
elif name == 'AuxiliarySarData':
|
||
filename = os.path.basename(paths[0])
|
||
part = filename.split("_")
|
||
name = part[5]
|
||
self.aux_key_word = name[0:8]
|
||
|
||
for path in paths:
|
||
shutil.copy(path, self.__workspace_slc_path)
|
||
processing_paras.update({'slc': self.__workspace_slc_path})
|
||
if name == 'DEM':
|
||
if para['DataType'] == 'File':
|
||
processing_paras.update({'DEM': para['ParaValue']})
|
||
else:
|
||
para_path_list = para['ParaValue'].split(";")
|
||
if len(para_path_list) != 0:
|
||
dem_path = os.path.join(self.__workspace_origin_path, para['ParaName'])
|
||
if os.path.exists(dem_path) is False:
|
||
os.mkdir(dem_path)
|
||
for file_path in para_path_list:
|
||
tif_name = os.path.basename(file_path)
|
||
shutil.copy(file_path, os.path.join(dem_path, tif_name))
|
||
para_path = os.path.join(self.__workspace_origin_path, para['ParaName'])
|
||
processing_paras.update({'DEM': para_path})
|
||
# # 解压DEM到指定文件夹
|
||
# path = para['ParaValue']
|
||
# import zipfile
|
||
# zip_file = zipfile.ZipFile(path)
|
||
# zip_list = zip_file.namelist() # 得到压缩包里所有文件
|
||
# for f in zip_list:
|
||
# zip_file.extract(f, self.__workspace_dem_path) # 循环解压文件到指定目录
|
||
# if os.path.splitext(f)[1] == '.wgs84':
|
||
# dem_name = f
|
||
# processing_paras.update({'dem': os.path.join(self.__workspace_dem_path, f)})
|
||
# zip_file.close()
|
||
# self.verifyAndModifyWgsXml(self.__workspace_dem_path + '\\' + dem_name + '.xml',
|
||
# self.__workspace_dem_path + '\\' + dem_name)
|
||
if name == 'Orbits':
|
||
if para['DataType'] == 'File':
|
||
processing_paras.update({'orbits': para['ParaValue']})
|
||
else:
|
||
para_path_list = para['ParaValue'].split(";")
|
||
if len(para_path_list) != 0:
|
||
dem_path = os.path.join(self.__workspace_origin_path, para['ParaName'])
|
||
if os.path.exists(dem_path) is False:
|
||
os.mkdir(dem_path)
|
||
for file_path in para_path_list:
|
||
tif_name = os.path.basename(file_path)
|
||
shutil.copy(file_path, os.path.join(dem_path, tif_name))
|
||
para_path = os.path.join(self.__workspace_origin_path, para['ParaName'])
|
||
processing_paras.update({'orbits': para_path})
|
||
if name == 'box':
|
||
datas = para['ParaValue'].split(';')
|
||
if len(datas) != 4:
|
||
msg = 'para: box is error!box:' + para['ParaValue']
|
||
raise Exception(msg)
|
||
box = datas[0] + ' ' + datas[1] + ' ' + datas[2] + ' ' + datas[3]
|
||
processing_paras.update({'box': box})
|
||
return processing_paras
|
||
|
||
def verifyAndModifyWgsXml(self, xmlPath, demPath):
|
||
import xml.dom.minidom as xmldom
|
||
domobj = xmldom.parse(xmlPath)
|
||
rootNode = domobj.documentElement
|
||
pathInxml = ''
|
||
# 获得子标签
|
||
propertyElementObj = rootNode.getElementsByTagName("property")
|
||
for property in propertyElementObj:
|
||
if property.hasAttribute("name"):
|
||
if property.getAttribute("name") == "file_name":
|
||
pathNode = property.getElementsByTagName("value")[0]
|
||
pathInxml = pathNode.childNodes[0].data
|
||
print('pathInxml1:', pathInxml)
|
||
pathNode.childNodes[0].data = r"/".join(demPath.split("\\"))
|
||
pathInxml = pathNode.childNodes[0].data
|
||
print('pathInxml2:', pathInxml)
|
||
with open(xmlPath, 'w') as f:
|
||
# 缩进换行编码
|
||
domobj.writexml(f, addindent=' ', encoding='utf-8')
|
||
|
||
def __create_work_space(self):
|
||
"""
|
||
删除原有工作区文件夹,创建新工作区文件夹
|
||
"""
|
||
self.__workspace_preprocessing_path = self.__workspace_tem_dir_path+"preprocessing""\\"
|
||
self.__workspace_preprocessed_path = self.__workspace_tem_dir_path+"preprocessed""\\"
|
||
self.__workspace_cut_path = self.__workspace_tem_dir_path+"cut""\\"
|
||
|
||
self.__workspace_resame2_path = self.__workspace_tem_dir_path+"resame2""\\"
|
||
self.__workspace_roi_path = self.__workspace_tem_dir_path+"roi""\\"
|
||
self.__workspace_atmos_dely_path = self.__workspace_tem_dir_path + "atmos_delay""\\"
|
||
self.__workspace_resame_itself_path = self.__workspace_tem_dir_path+"resame_itself""\\"
|
||
|
||
self.__workspace_Output_path = self.__workspace_path + EXE_NAME + r"\Output""\\"
|
||
self.__workspace_iscePreprocessed_path = self.__workspace_path + EXE_NAME + r"\Temporary\isce\preprocessed""\\"
|
||
self.__workspace_isceProcessing_path = self.__workspace_path + EXE_NAME + r"\Temporary\isce\processing""\\"
|
||
self.__workspace_isce_path = os.path.join(self.__workspace_isceProcessing_path, 'isce_workspace')
|
||
self.__workspace_mintpy_path = os.path.join(self.__workspace_isceProcessing_path, 'mintpy_workspace')
|
||
self.__workspace_dem_path = os.path.join(self.__workspace_iscePreprocessed_path, 'dem')
|
||
self.__workspace_slc_path = self.__workspace_path + EXE_NAME + r"\Temporary\slc"
|
||
self.__workspace_origin_path = os.path.join(self.__workspace_path, EXE_NAME, "Temporary", "origin")
|
||
|
||
path_list = [self.__workspace_preprocessing_path, self.__workspace_preprocessed_path,
|
||
self.__workspace_cut_path, self.__workspace_resame2_path,
|
||
self.__workspace_roi_path, self.__workspace_atmos_dely_path,
|
||
self.__workspace_resame_itself_path, self.__workspace_Output_path,
|
||
self.__workspace_iscePreprocessed_path, self.__workspace_isceProcessing_path,
|
||
self.__workspace_dem_path, self.__workspace_isce_path,
|
||
self.__workspace_slc_path, self.__workspace_mintpy_path,
|
||
self.__workspace_origin_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!')
|
||
|
||
def get_prodataset(self, key_name):
|
||
"""
|
||
1、读取气象数据，另存为tif；2、输出每个气象数据文件
|
||
:param key_name :["MasterNC", "AuxiliaryNC"]
|
||
"""
|
||
nc_file = self.__processing_paras[key_name]
|
||
out_path = self.__workspace_preprocessing_path
|
||
self.NcHandle.tran_tif(nc_file, out_path, key_name)
|
||
|
||
def read_nc_into_array(self):
|
||
"""
|
||
1、读取气象数据（.nc）中气象参数geo、level、hum、temp;2、以数组的形式输出。
|
||
"""
|
||
nc_list = ["AuxiliaryNC_geo_h", "AuxiliaryNC_re_hum", "AuxiliaryNC_temp",
|
||
"MasterNC_geo_h", "MasterNC_re_hum", "MasterNC_temp"]
|
||
self.resampling_byitself(nc_list, self.__preprocessed_paras2)
|
||
|
||
arr_dict = {}
|
||
arr_name = ["a_geo", "a_re_hum", "a_temp", "m_geo", "m_re_hum", "m_temp"]
|
||
path = self.__workspace_resame_itself_path
|
||
for i in range(0, len(nc_list)):
|
||
in_path = path+nc_list[i]+".tif"
|
||
array = self.imageHandler.get_all_band_array(in_path)
|
||
arr_dict.update({arr_name[i]: array})
|
||
logger.info('read the weather data into arrays success!')
|
||
return arr_dict["a_geo"], arr_dict["a_re_hum"], arr_dict["a_temp"], arr_dict["m_geo"], arr_dict["m_re_hum"], arr_dict["m_temp"]
|
||
|
||
@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):
|
||
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_names1 =[]
|
||
para_names1 = ["m_filt_fine", "m_incLocal", "a_incLocal", "m_hgt"]
|
||
para_names22 = ["MasterNC_geo_h", "MasterNC_re_hum", "MasterNC_temp",
|
||
"AuxiliaryNC_geo_h", "AuxiliaryNC_re_hum", "AuxiliaryNC_temp"] # 气象数据名称
|
||
# 1/5、检查微波影像、气象数据坐标系
|
||
self.pro_img_path = self.check_img_projection(para_names1, para_names22)
|
||
|
||
# 2/5、气象数据相交区域矢量2;
|
||
self.intersect_shp_path = self.__workspace_preprocessed_path + "IntersectPolygon2.shp"
|
||
if self.cal_intersect_shp(self.intersect_shp_path, self.pro_img_path, para_names22) is False:
|
||
logger.info('ERA data not overlap!')
|
||
return False
|
||
logger.info('create intersect shp success!')
|
||
# 判断气象影像是否覆盖微波影像
|
||
if self.intersect_judge(self.pro_img_path, para_names22, para_names1) is False: # 检查气象数据是否覆盖干涉图
|
||
logger.info('ERA data not fully cover the interferogram!')
|
||
return False
|
||
|
||
# 3/5、矢量1裁剪微波影像; 矢量2裁剪气象影像
|
||
cutted_img_path2 = self.cut_img(para_names22, self.intersect_shp_path)
|
||
if cutted_img_path2 is None:
|
||
return False
|
||
logger.info('cut images success!')
|
||
|
||
self.force_del_file(self.__workspace_preprocessing_path) # 删除气象数据tif
|
||
|
||
# 4/5、气象影像重采样
|
||
self.__preprocessed_paras = {}
|
||
self.__preprocessed_paras.update({"m_filt_fine": self.pro_img_path["m_filt_fine"]})
|
||
self.__preprocessed_paras.update({"m_incLocal": self.pro_img_path["m_incLocal"]})
|
||
self.__preprocessed_paras.update({"a_incLocal": self.pro_img_path["a_incLocal"]})
|
||
self.__preprocessed_paras.update({"m_hgt": self.pro_img_path["m_hgt"]})
|
||
|
||
self.__preprocessed_paras2 = self.resampling_img(para_names22, cutted_img_path2, self.__workspace_resame2_path,
|
||
cutted_img_path2["MasterNC_geo_h"])
|
||
# 对主辅影像入射角进行重采样
|
||
incLocal_param = ["m_incLocal", "a_incLocal"]
|
||
m_hgt = ["m_hgt"]
|
||
m_incLocal_resample, a_incLocal_resample = self.reasmple_incLocal(incLocal_param, m_hgt)
|
||
self.__preprocessed_paras.update({"m_incLocal": m_incLocal_resample})
|
||
self.__preprocessed_paras.update({"a_incLocal": a_incLocal_resample})
|
||
logger.info("incLocal reasmple success!")
|
||
|
||
# 5/5、计算掩膜
|
||
self.check_LocalIncidenceAngle()
|
||
# self.tif_angle_mask_path = self.create_roi()
|
||
|
||
def reasmple_incLocal(self, para_name1, para_name2):
|
||
m_incLocal = self.__preprocessed_paras[para_name1[0]]
|
||
a_incLocal = self.__preprocessed_paras[para_name1[1]]
|
||
m_hgt = self.__preprocessed_paras[para_name2[0]]
|
||
m_incLocal_resample = os.path.join(self.__workspace_slc_path, "m_incLocal_resample.tif")
|
||
a_incLocal_resample = os.path.join(self.__workspace_slc_path, "a_incLocal_resample.tif")
|
||
pp.resampling_by_scale(m_incLocal, m_incLocal_resample, m_hgt)
|
||
pp.resampling_by_scale(a_incLocal, a_incLocal_resample, m_hgt)
|
||
return m_incLocal_resample, a_incLocal_resample
|
||
|
||
|
||
def check_img_projection(self, para_names1, para_names2):
|
||
"""
|
||
读取每一张图像,检查图像坐标系，如果是平面坐标，则转换为WGS84经纬度坐标
|
||
:param para_names1:需要检查的参数名称
|
||
:param para_names2:需要检查的参数名称
|
||
"""
|
||
if len(para_names1) == 0 or len(para_names2) == 0:
|
||
return False
|
||
pro_img_path = {}
|
||
for name1 in para_names1:
|
||
# dem、入射角
|
||
in_para = self.__processing_paras[name1]
|
||
proj = self.imageHandler.get_projection(in_para)
|
||
keyword = proj.split("[", 2)[0]
|
||
out_para = self.__workspace_preprocessed_path + "pro_"+os.path.basename(self.__processing_paras[name1])
|
||
if keyword == "GEOGCS":
|
||
shutil.copy(in_para, out_para)
|
||
elif keyword == "PROJCS":
|
||
pp.trans_epsg4326(out_para, in_para)
|
||
elif len(keyword) == 0 or keyword.strip() == "" or keyword.isspace() is True:
|
||
raise Exception('the coordinate of dem or incident is missing!')
|
||
pro_img_path.update({name1: out_para})
|
||
logger.info('image coordinate transformation finished!')
|
||
for name2 in para_names2:
|
||
in_para = self.__workspace_preprocessing_path+name2+".tif"
|
||
proj = self.imageHandler.get_projection(in_para)
|
||
keyword = proj.split("[", 2)[0]
|
||
out_para = self.__workspace_preprocessed_path+name2+".tif"
|
||
if keyword == "GEOGCS":
|
||
shutil.copy(in_para, out_para)
|
||
elif keyword == "PROJCS":
|
||
pp.trans_epsg4326(out_para, in_para)
|
||
elif len(keyword) == 0 or keyword.strip() == "" or keyword.isspace() is True:
|
||
raise Exception('data miss coordinate!')
|
||
pro_img_path.update({name2: out_para})
|
||
logger.info('ERA data coordinate transformation finished!')
|
||
return pro_img_path
|
||
|
||
def cal_intersect_shp(self, shp_path, para_dict, para_names):
|
||
"""
|
||
:param shp_path:相交区域矢量文件保存区域
|
||
:param para_dict: 字典 {影像名称:路径}
|
||
: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(para_dict[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, 4326) is False:
|
||
return False
|
||
logger.info('calculate intersect finished!')
|
||
return True
|
||
|
||
def intersect_judge(self, para_dict, para_names, para_names2):
|
||
"""
|
||
:param para_dict: 字典 {影像名称:路径}
|
||
:param para_names:判断相交影像的名称列表
|
||
:param para_names2:判断相交影像的名称列表
|
||
:return: True or False
|
||
"""
|
||
if len(para_names) == 0 or len(para_names2) == 0:
|
||
return False
|
||
scopes_01 = ()
|
||
for name in para_names:
|
||
scope_tuple_01 = (self.imageHandler.get_scope(para_dict[name]),)
|
||
scopes_01 += scope_tuple_01
|
||
intersect_polygon_01 = pp.intersect_polygon(scopes_01)
|
||
scopes_02 = ()
|
||
for name in para_names2:
|
||
scope_tuple_02 = (self.imageHandler.get_scope(para_dict[name]),)
|
||
scopes_02 += scope_tuple_02
|
||
intersect_polygon_02 = pp.intersect_polygon(scopes_02)
|
||
intersect_polygon_01 = np.array(intersect_polygon_01)
|
||
intersect_polygon_02 = np.array(intersect_polygon_02)
|
||
flag_left, flag_right, flag_down, flag_up = 0, 0, 0, 0
|
||
if intersect_polygon_02[0][0] >= intersect_polygon_01[0][0]:
|
||
flag_left = 1
|
||
if intersect_polygon_02[2][0] <= intersect_polygon_01[2][0]:
|
||
flag_right = 1
|
||
if intersect_polygon_02[0][1] >= intersect_polygon_01[0][1]:
|
||
flag_down = 1
|
||
if intersect_polygon_02[2][1] <= intersect_polygon_01[2][1]:
|
||
flag_up = 1
|
||
if flag_left == 0 or flag_right == 0 or flag_down == 0 or flag_up == 0:
|
||
return False
|
||
return True
|
||
|
||
def cut_img(self, para_names1, shp_path):
|
||
"""
|
||
使用矢量数据裁剪影像
|
||
:param para_names1:需要检查的参数名称
|
||
:param shp_path：裁剪的shp文件
|
||
"""
|
||
if len(para_names1) == 0:
|
||
return None
|
||
cutted_img_path = {}
|
||
# noinspection PyBroadException
|
||
try:
|
||
for name1 in para_names1:
|
||
input_path = self.pro_img_path[name1]
|
||
output_path = self.__workspace_cut_path + name1 + ".tif"
|
||
pp.cut_img(output_path, input_path, shp_path)
|
||
cutted_img_path.update({name1: output_path})
|
||
logger.info("cut %s success!", name1)
|
||
except BaseException:
|
||
logger.error('cut_img failed!')
|
||
return None
|
||
logger.info('cut image finished!')
|
||
return cutted_img_path
|
||
|
||
@staticmethod
|
||
def resampling_img(para_names1, img_paths, out_path, refer_img_path):
|
||
"""
|
||
以主影像为参考，对影像重采样
|
||
:param para_names1:需要检查的参数名称
|
||
:param img_paths：待重采样影像路径
|
||
:param out_path：保存数据的文件夹路径
|
||
:param refer_img_path：参考影像路径
|
||
"""
|
||
if len(para_names1) == 0 or len(img_paths) == 0:
|
||
return
|
||
prepro_imgs_path = {}
|
||
for name1 in para_names1:
|
||
img_path = img_paths[name1]
|
||
output_para = out_path+name1+".tif"
|
||
|
||
pp.resampling_by_scale(img_path, output_para, refer_img_path)
|
||
prepro_imgs_path.update({name1: output_para})
|
||
logger.info("resampling %s success!", name1)
|
||
logger.info('resampling finished !')
|
||
return prepro_imgs_path
|
||
|
||
def check_LocalIncidenceAngle(self):
|
||
"""
|
||
将角度的无效值设置为nan，把角度值转为弧度值
|
||
"""
|
||
# 主影像雷达入射角
|
||
M_Angle_path=self.__preprocessed_paras['m_incLocal']
|
||
out_M_Angle_path = os.path.splitext(M_Angle_path)[0] + r"_check.tif"
|
||
pp.check_LocalIncidenceAngle(out_M_Angle_path, M_Angle_path)
|
||
if os.path.isfile(M_Angle_path) is True:
|
||
os.remove(M_Angle_path)
|
||
self.__preprocessed_paras['m_incLocal'] = out_M_Angle_path
|
||
|
||
# 辅影像雷达入射角
|
||
A_Angle_path=self.__preprocessed_paras['a_incLocal']
|
||
out_A_Angle_path = os.path.splitext(A_Angle_path)[0] + r"_check.tif"
|
||
pp.check_LocalIncidenceAngle(out_A_Angle_path, A_Angle_path)
|
||
if os.path.isfile(A_Angle_path) is True:
|
||
os.remove(A_Angle_path)
|
||
self.__preprocessed_paras['a_incLocal'] = out_A_Angle_path
|
||
|
||
# def create_roi(self):
|
||
# """
|
||
# 计算 ROI
|
||
# """
|
||
# # 利用主、辅雷达入射角的有效范围生成MASK
|
||
# master_angle_mask_path = self.__workspace_roi_path + "mas_angle_mask.tif"
|
||
# ROIAlg.trans_tif2mask(master_angle_mask_path, self.__preprocessed_paras["m_incLocal"], np.nan)
|
||
# auxiliary_angle_mask_path = self.__workspace_roi_path + "aux_angle_mask.tif"
|
||
# ROIAlg.trans_tif2mask(auxiliary_angle_mask_path, self.__preprocessed_paras["a_incLocal"], np.nan)
|
||
#
|
||
# tif_angle_mask_path = self.__workspace_roi_path + "tif_angle_mask.tif"
|
||
# ROIAlg.combine_mask(tif_angle_mask_path, master_angle_mask_path, auxiliary_angle_mask_path)
|
||
#
|
||
# logger.info('combine_mask success!')
|
||
# logger.info('create ROI image success!')
|
||
#
|
||
# # 计算roi区域
|
||
# apath=self.__preprocessed_paras["m_incLocal"]
|
||
# ROIAlg.cal_roi(apath, apath, tif_angle_mask_path, background_value=1)
|
||
# apath2=self.__preprocessed_paras["a_incLocal"]
|
||
# ROIAlg.cal_roi(apath2, apath2, tif_angle_mask_path, background_value=1)
|
||
# logger.info('ROI cai image success!')
|
||
# return tif_angle_mask_path
|
||
|
||
def resampling_byitself(self, para_names, img_paths):
|
||
"""
|
||
将影像重采样为正方形
|
||
param: para_names 存放气象数据名称的列表
|
||
param: img_paths：字典{气象数据名称：路径}
|
||
"""
|
||
for name1 in para_names:
|
||
img_path = img_paths[name1]
|
||
dataset = gdal.Open(img_path)
|
||
ref_cols = dataset.RasterXSize # 列数
|
||
ref_rows = dataset.RasterYSize # 行数
|
||
band_count = dataset.RasterCount # 波段数
|
||
geotrans = list(dataset.GetGeoTransform())
|
||
output_path = self.__workspace_resame_itself_path+name1+".tif"
|
||
if ref_cols > ref_rows:
|
||
lagest_len = ref_cols
|
||
scale = ref_cols/ref_rows
|
||
geotrans[5] = geotrans[5] / scale # 像元高度变为原来的scale倍
|
||
else:
|
||
lagest_len = ref_rows
|
||
scale = ref_rows/ref_cols
|
||
geotrans[1] = geotrans[1] / scale # 像元宽度变为原来的scale倍
|
||
band1 = dataset.GetRasterBand(1)
|
||
data_type = band1.DataType
|
||
target = dataset.GetDriver().Create(output_path, xsize=lagest_len, ysize=lagest_len, bands=band_count,
|
||
eType=data_type)
|
||
target.SetProjection(dataset.GetProjection()) # 设置投影坐标
|
||
target.SetGeoTransform(geotrans) # 设置地理变换参数
|
||
|
||
total = band_count + 1
|
||
for index in range(1, total):
|
||
# 读取波段数据
|
||
data = dataset.GetRasterBand(index).ReadAsArray(buf_xsize=lagest_len, buf_ysize=lagest_len)
|
||
out_band = target.GetRasterBand(index)
|
||
|
||
no_data_value = dataset.GetRasterBand(index).GetNoDataValue() # 获取没有数据的点
|
||
if not (no_data_value is None):
|
||
out_band.SetNoDataValue(no_data_value)
|
||
|
||
out_band.WriteArray(data) # 写入数据到新影像中
|
||
out_band.FlushCache()
|
||
out_band.ComputeBandStats(False) # 计算统计信息
|
||
del dataset
|
||
del target
|
||
|
||
def resampe_image(self, aux_dry_wet_unresame_path, mas_dry_wet_unresame_path):
|
||
"""
|
||
将延迟图（ztd）恢复为非正方形的样子
|
||
para: aux_dry_wet_unresame_path：待裁剪的辅影像延迟图路径
|
||
para: mas_dry_wet_unresame_path：待裁剪的主影像延迟图路径
|
||
return :裁剪后的辅影像延迟图路径、主影像延迟图路径
|
||
"""
|
||
# 1.3 裁剪主辅影像延迟值
|
||
out_aux_dry_wet_resamed_path = self.__workspace_cut_path + "resamed_aux_dry_wet.tif"
|
||
out_mas_dry_wet_resamed_path = self.__workspace_cut_path + "resamed_mas_dry_wet.tif"
|
||
|
||
pp.resampling_by_scale(aux_dry_wet_unresame_path, out_aux_dry_wet_resamed_path,
|
||
self.__workspace_resame2_path + "MasterNC_geo_h.tif")
|
||
pp.resampling_by_scale(mas_dry_wet_unresame_path, out_mas_dry_wet_resamed_path,
|
||
self.__workspace_resame2_path + "MasterNC_geo_h.tif")
|
||
|
||
logger.info('resampling the aux correction and mas correction finished !')
|
||
return out_aux_dry_wet_resamed_path, out_mas_dry_wet_resamed_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 delay_conversion(self, out_aux_ztd_path, out_mas_ztd_path):
|
||
"""
|
||
延迟转换
|
||
param:out_aux_ztd_path 输出：辅影像延迟图路径
|
||
param:out_mas_ztd_path 输出：主影像延迟图路径
|
||
"""
|
||
self.__proj= ImageHandler().get_projection(out_mas_ztd_path)
|
||
self.__geo=ImageHandler().get_geotransform(out_mas_ztd_path)
|
||
|
||
out_aux_file = self.__preprocessed_paras["a_incLocal"]
|
||
out_mas_file = self.__preprocessed_paras["m_incLocal"]
|
||
|
||
# 计算延迟相位（数组）
|
||
xml_files = list(glob.glob(os.path.join(self.__workspace_isce_path, "reference", "*.xml")))
|
||
# mas_lambd = Ady().get_radarwavelength(os.path.join(self.__workspace_isce_path, "reference", "IW1.xml")) # 波长
|
||
mas_lambd = Ady().get_radarwavelength(xml_files[0]) # 波长
|
||
second_path = os.path.join(self.__workspace_isce_path, "secondarys")
|
||
floder_name = os.listdir(second_path)[0]
|
||
|
||
xmls = list(glob.glob(os.path.join(second_path, floder_name, "*.xml")))
|
||
aux_lambd = Ady().get_radarwavelength(xmls[0]) # 波长
|
||
|
||
mas_phase_bias_array = Ady().cal_phase_bias(out_mas_file, out_mas_ztd_path, mas_lambd) # 主影像延迟相位数组
|
||
aux_phase_bias_array = Ady().cal_phase_bias(out_aux_file, out_aux_ztd_path, aux_lambd) # 辅影像延迟相位数组
|
||
|
||
# 主辅影像延迟相位存为tif
|
||
mas_phase_bias_path = self.__workspace_atmos_dely_path + "MasterPhaseBias" + ".tif" # 输出的产品数据1/3
|
||
aux_phase_bias_path = self.__workspace_atmos_dely_path + "AuxiliaryPhaseBias" + ".tif" # 输出的产品数据2/3
|
||
|
||
ImageHandler().write_img(mas_phase_bias_path, self.__proj, self.__geo, mas_phase_bias_array)
|
||
ImageHandler().write_img(aux_phase_bias_path, self.__proj, self.__geo, aux_phase_bias_array)
|
||
|
||
self.imageHandler.write_quick_view(mas_phase_bias_path) # 输出的快视图1/3
|
||
self.imageHandler.write_quick_view(aux_phase_bias_path) # 输出的快视图2/3
|
||
return mas_phase_bias_path, aux_phase_bias_path,aux_lambd
|
||
|
||
def creat_xml(self, product_path, start):
|
||
"""
|
||
生成元文件案例
|
||
product_path: 大气延迟校正产品输出的影像文件路径
|
||
"""
|
||
xml_path = "./model_meta.xml"
|
||
tem_folder = self.__workspace_tem_dir_path
|
||
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, self.processinfo, out_path1, out_path2).calu_nature(start)
|
||
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()
|
||
return model_xml_path
|
||
|
||
def transformZtd(self, scopes, ori_sim, out_mas_ztd_path, geo_mas_ztd_path):
|
||
self._tr = TransImgL1A(ori_sim, scopes)
|
||
self._tr.l1a_2_geo(ori_sim, out_mas_ztd_path, geo_mas_ztd_path)
|
||
|
||
def get_scopes(self, ori_sim):
|
||
ori_sim_data = self.imageHandler.get_data(ori_sim)
|
||
lon = ori_sim_data[0, :, :]
|
||
lat = ori_sim_data[1, :, :]
|
||
|
||
min_lon = np.nanmin(lon)
|
||
max_lon = np.nanmax(lon)
|
||
min_lat = np.nanmin(lat)
|
||
max_lat = np.nanmax(lat)
|
||
|
||
scopes = [[min_lon, max_lat], [max_lon, max_lat], [min_lon, min_lat], [max_lon, min_lat]]
|
||
return scopes
|
||
|
||
def prod_handle(self):
|
||
"""
|
||
干涉测量相位大气延迟校正处理
|
||
"""
|
||
hh_flag, hv_flag, vh_flag, vv_flag = 0, 0, 0, 0
|
||
xml_files = list(glob.glob(os.path.join(self.__workspace_isce_path, "reference", "*.xml")))
|
||
if len(xml_files) == 0:
|
||
raise Exception("please check box!")
|
||
# polarization = Ady().get_polarization(os.path.join(self.__workspace_isce_path, "reference", "IW1.xml")) # 极化方式
|
||
polarization = Ady().get_polarization(xml_files[0]) # 极化方式
|
||
if 'hh' in polarization or 'HH' in polarization:
|
||
hh_flag = 1
|
||
elif 'hv' in polarization or 'HV' in polarization:
|
||
hv_flag = 1
|
||
elif 'vh' in polarization or 'VH' in polarization:
|
||
vh_flag = 1
|
||
elif 'vv' in polarization or 'VV' in polarization:
|
||
vv_flag = 1
|
||
self.processinfo = [hh_flag, hv_flag, vh_flag, vv_flag]
|
||
|
||
# 1.1 数据准备（三种气象数据的数组、dem路径、参考影像）
|
||
a_geo, a_re_hum, a_temp, m_geo, m_re_hum, m_temp = self.read_nc_into_array()
|
||
dem_file = self.__preprocessed_paras["m_hgt"]
|
||
base_file = self.__workspace_resame_itself_path + "MasterNC_geo_h.tif" # 用来获取裁剪后气象数据的经纬度、分辨率等信息
|
||
|
||
# 1.2 计算主、辅影像垂直分层湿干延迟值，并保存干湿延迟图
|
||
m_cdstack_dry_wet = Ady().aps_wrf_sar(base_file, m_temp, m_re_hum, m_geo, dem_file) # 主影像垂直分层干湿延迟值数组
|
||
m_dry_wet_path = self.__workspace_cut_path + "m_dry_wet.tif" # 主影像干湿延迟图保存路径
|
||
Ady().write_dry_wet_tif(dem_file, base_file, m_cdstack_dry_wet, m_dry_wet_path)
|
||
|
||
a_cdstack_dry_wet = Ady().aps_wrf_sar(base_file, a_temp, a_re_hum, a_geo, dem_file) # 辅影像垂直分层干湿延迟值数组
|
||
a_dry_wet_path = self.__workspace_cut_path + "a_dry_wet.tif" # 辅影像干湿延迟图保存路径
|
||
# Ady().write_tif(dem_file, a_cdstack_dry_wet, a_dry_wet_path)
|
||
Ady().write_dry_wet_tif(dem_file, base_file, a_cdstack_dry_wet, a_dry_wet_path)
|
||
|
||
# 1.3 主、辅影像干湿延迟图重采样为原来（非正方形）的形状
|
||
aux_dry_wet_resamed_path, mas_dry_wet_resamed_path = self.resampe_image(a_dry_wet_path,
|
||
m_dry_wet_path) # 参照气象数据分辨率进行重采样
|
||
# del a_cdstack_dry_wet, m_cdstack_dry_wet, a_geo, a_re_hum, a_temp, m_geo, m_re_hum, m_temp
|
||
|
||
# self.del_folder(self.__workspace_resame_itself_path)
|
||
# self.del_folder(self.__workspace_resame2_path)
|
||
|
||
logger.info('start calc ztd ...')
|
||
gc.collect() # 回收内存
|
||
# 1.4 输出主辅影像ztd数组
|
||
a_ztd = Ady().calc_ztd(a_dry_wet_path, dem_file) # 辅影像的ztd数组
|
||
m_ztd = Ady().calc_ztd(m_dry_wet_path, dem_file) # 主影像的ztd数组
|
||
out_aux_ztd_path = self.__workspace_cut_path + "aux_ztd.tif"
|
||
out_mas_ztd_path = self.__workspace_cut_path + "mas_ztd.tif"
|
||
Ady().write_ztd_tif(dem_file, base_file, a_ztd, out_aux_ztd_path) # ztd数组->ztd影像
|
||
Ady().write_ztd_tif(dem_file, base_file, m_ztd, out_mas_ztd_path) # ztd数组->ztd影像
|
||
|
||
geo_mas_ztd_path = self.__workspace_atmos_dely_path + "MasterZtd.tif"
|
||
geo_aux_ztd_path = self.__workspace_atmos_dely_path + "AuxiliaryZtd.tif"
|
||
|
||
# todo 将ztd结果进行地理编码并放到结果包中
|
||
# 生成坐标映射文件
|
||
limit_lat = os.path.join(self.__workspace_isce_path, "merged", "geom_reference", "lat.rdr")
|
||
limit_lon = os.path.join(self.__workspace_isce_path, "merged", "geom_reference", "lon.rdr")
|
||
ori_sim = self.__workspace_cut_path + 'ori_sim.tif'
|
||
self.imageHandler.band_merge(limit_lon, limit_lat, ori_sim)
|
||
if not os.path.exists(ori_sim):
|
||
logger.error("create ori_sim failed")
|
||
scopes = self.imageHandler.get_scopes(ori_sim)
|
||
self.transformZtd(scopes, ori_sim, out_mas_ztd_path, geo_mas_ztd_path) # 主影像ztd进行地理编码
|
||
self.transformZtd(scopes, ori_sim, out_aux_ztd_path, geo_aux_ztd_path) # 辅影像ztd进行地理编码
|
||
|
||
# 2.1 计算延迟转换, 输出两个影像相位误差图
|
||
M_PhaseError_path, A_PhaseError_path,aux_lambd = self.delay_conversion(out_aux_ztd_path, out_mas_ztd_path)
|
||
|
||
# 干涉图
|
||
Filt_path = self.__preprocessed_paras["m_filt_fine"]
|
||
|
||
# 2.2 延迟相位去除
|
||
self.atmospheric_product_path = self.__workspace_atmos_dely_path + "AtmosphericDelay"+".tif"
|
||
|
||
bias_result = Ady().complex_phase_removal(M_PhaseError_path, A_PhaseError_path, Filt_path, aux_lambd)
|
||
ImageHandler.write_img(self.atmospheric_product_path, self.__proj, self.__geo, bias_result) # 剔除相位延迟后的数组->剔除相位延迟干涉图
|
||
# self.force_del_file(self.__workspace_cut_path)
|
||
logger.info('progress bar :77％')
|
||
|
||
def isce_stackSentinel(self, slc_dir, dem_path, aux_dir, orbits_dir, isce_work_space, key_word, isce_exe_dir,box):
|
||
"""执行 stackSentinel.exe"""
|
||
os.chdir(isce_exe_dir)
|
||
# cmd = "stackSentinel.exe -s {} -d {} -a {} -o {} -w {} -n {} -m {} -c {} --exeabsolute_dir {} -b {}".format(slc_dir, dem_path, aux_dir, orbits_dir, isce_work_space,"'1 2 3'", key_word, '1', isce_exe_dir,box)
|
||
# cmd = "stackSentinel.exe -s {} -d {} -a {} -o {} -w {} -n {} -m {} -W {} -c {} --exeabsolute_dir {} -b {}".format(slc_dir, dem_path, aux_dir, orbits_dir, isce_work_space,"'1 2 3'", key_word,'interferogram','6', isce_exe_dir,box)
|
||
cmd = "stackSentinel.exe -s {} -d {} -a {} -o {} -w {} -n '1 2 3' -m {} -W 'interferogram' -c '1' --exeabsolute_dir {} -b {}".format(slc_dir, dem_path, aux_dir, orbits_dir, isce_work_space, key_word, isce_exe_dir, box)
|
||
# cmd = "stackSentinel.exe -s {} -d {} -a {} -o {} -w {} -n '1 2 3' -m {} -W 'interferogram' -c '1' --exeabsolute_dir {} -b '37.3 39.4 108.10 111.82'".format(slc_dir, dem_path, aux_dir, orbits_dir, isce_work_space, key_word, isce_exe_dir)
|
||
|
||
#stackSentinel.exe -s I:/MicroProduct/daqiceshishuju_20221109/slc -d I:/MicroProduct/daqiceshishuju_20221109/DEM/demLat_N37_N40_Lon_E107_E112.dem.wgs84 -a I:/MicroProduct/daqiceshishuju_20221109/AuxDir -o I:/MicroProduct/daqiceshishuju_20221109/orbits -w I:/MicroProduct/daqiceshishuju_20221109/daqi_workspace/isce_workspace -n '1 2 3' -m 20190206 -W interferogram -c 6 --exeabsolute_dir I:/MicroProduct/microproduct/atmosphericDelay/ISCEApp -b '37.3 39.4 108.10 111.82'
|
||
logger.info('stackSentinel_cmd:{}'.format(cmd))
|
||
result = os.system(cmd)
|
||
logger.info('cmd_result:{}'.format(result))
|
||
logger.info('stackSentinel finish!')
|
||
|
||
def isce_autorun(self, isce_exe_dir, isce_work_space):
|
||
"""执行 autorun.py"""
|
||
#cmd = 'python {}/autorun.py -e {} -o {}'.format(isce_exe_dir, isce_exe_dir, isce_work_space)
|
||
#注释by-shj-解决打包错误
|
||
# cmd = 'python autorun.py -e {} -o {} -cygwinFlag False'.format(isce_exe_dir, isce_work_space)
|
||
# logger.info('autorun_cmd:{}'.format(cmd))
|
||
# result = os.system(cmd)
|
||
|
||
|
||
isce_exe_dir = r"/".join(os.path.join(self.env_str, "ISCEApp").split("\\"))
|
||
os.chdir(isce_exe_dir)
|
||
cmd = ['-e', isce_exe_dir, '-o', isce_work_space]
|
||
logger.info('autorun_cmd:{}'.format(cmd))
|
||
auto_run_main(cmd)
|
||
os.chdir(self.env_str)
|
||
# logger.info('cmd_result:{}'.format(result))
|
||
# if result != 0:
|
||
# raise Exception('autorun.py run failed!')
|
||
logger.info('autorun_cmd success!')
|
||
|
||
def isce_write(self, in_dict, out_dict):
|
||
"""isce(.rdr、.vrt)-->.tif"""
|
||
for i in range(0, len(in_dict)):
|
||
output_file = os.path.join(self.__workspace_slc_path, out_dict[i])
|
||
limit_lat = os.path.join(self.__workspace_isce_path, "merged", "geom_reference", "lat.rdr")
|
||
limit_lon = os.path.join(self.__workspace_isce_path, "merged", "geom_reference", "lon.rdr")
|
||
if in_dict[i]!="filt_fine.int.vrt" and in_dict[i] != "filt_fine.unw.vrt":
|
||
input_file = os.path.join(self.__workspace_isce_path, "merged", "geom_reference", in_dict[i])
|
||
if os.path.isfile(input_file):
|
||
DataTransform().write_tif(output_file, limit_lat, limit_lon, input_file)
|
||
else:
|
||
raise Exception('can not find file: ', input_file)
|
||
elif in_dict[i] == "filt_fine.int.vrt":
|
||
filename = os.listdir(os.path.join(self.__workspace_isce_path, "merged", "interferograms"))[0]
|
||
input_file = os.path.join(self.__workspace_isce_path, "merged", "interferograms", filename, in_dict[i])
|
||
if os.path.isfile(input_file):
|
||
DataTransform().write_tif(output_file, limit_lat, limit_lon, input_file)
|
||
else:
|
||
raise Exception('can not find file: ', input_file)
|
||
elif in_dict[i] == "filt_fine.unw.vrt":
|
||
filename = os.listdir(os.path.join(self.__workspace_isce_path, "merged", "interferograms"))[0]
|
||
input_file = os.path.join(self.__workspace_isce_path, "merged", "interferograms", filename, in_dict[i])
|
||
if os.path.isfile(input_file):
|
||
DataTransform().unw_write_tif(output_file, limit_lat, limit_lon, input_file)
|
||
else:
|
||
raise Exception('can not find file: ', input_file)
|
||
name = out_dict[i].split(".")[0]
|
||
self.__processing_paras.update({name: output_file})
|
||
|
||
def isce_run_steps(self, run_steps, target):
|
||
for i in range(0, len(run_steps)):
|
||
uwm_file = os.path.join(self.__workspace_isce_path, "run_files", run_steps[i])
|
||
shutil.move(uwm_file, target)
|
||
|
||
def isce_del_dir(self, path_list):
|
||
for path in path_list:
|
||
if os.path.exists(path):
|
||
if os.path.isdir(path):
|
||
shutil.rmtree(path)
|
||
os.makedirs(path)
|
||
logger.info('create new floder success!')
|
||
|
||
def get_interv_img(self):
|
||
"""根据干涉形变算法生成干涉图"""
|
||
# 执行isce2.5生成干涉图
|
||
# 生成工作流配置文件
|
||
slc_dir = r"/".join(self.__processing_paras['slc'].split("\\"))
|
||
|
||
dem_dir = r"/".join(self.__processing_paras['DEM'].split("\\"))
|
||
aux_dir = r"/".join(os.path.join(self.env_str, r'AuxDir').split("\\"))
|
||
orbits_dir = r"/".join(self.__processing_paras['orbits'].split("\\"))
|
||
isce_work_space = r"/".join(self.__workspace_isce_path.split("\\"))
|
||
box = "'" + self.__processing_paras['box'] + "'"
|
||
isce_exe_dir = r"/".join(os.path.join(self.env_str, "ISCEApp").split("\\"))
|
||
target = isce_work_space
|
||
|
||
os.chdir(isce_exe_dir)
|
||
#前处理 转换tif影像为 wgs84格式
|
||
dem_dir = '/cygdrive/' + dem_dir.replace(":/", "/")
|
||
out_dem_dir = self.__workspace_dem_path
|
||
out_dem_dir = '/cygdrive/' + out_dem_dir.replace(":\\", "/")
|
||
cmd = "demhgt2wgs.exe --tif_path {} --hgt_path {} --ASTGTM2".format(dem_dir, out_dem_dir)
|
||
logger.info('demhgt2wgs_cmd:{}'.format(cmd))
|
||
result = os.system(cmd)
|
||
logger.info('cmd_result:{}'.format(result))
|
||
import glob
|
||
in_tif_paths = list(glob.glob(os.path.join(self.__workspace_dem_path, '*.wgs84')))
|
||
if in_tif_paths == []:
|
||
raise Exception('demhgt2wgs.exe run failed!')
|
||
dem_path = r"/".join(in_tif_paths[0].split("\\"))
|
||
logger.info('demhgt2wgs finish!')
|
||
logger.info('progress bar: 5%')
|
||
|
||
|
||
# 一、 辅影像
|
||
# 1.1 执行stackSentinel.exe 生成运行配置文件
|
||
# key_word = "20190206"
|
||
key_word = self.aux_key_word
|
||
self.isce_stackSentinel(slc_dir, dem_path, aux_dir, orbits_dir, isce_work_space, key_word, isce_exe_dir,box)
|
||
|
||
# 1.2 执行autorun.py 生成干涉图
|
||
run_steps = ["run_13_generate_burst_igram", "run_14_merge_burst_igram",
|
||
"run_15_filter_coherence", "run_16_unwrap"]
|
||
self.isce_run_steps(run_steps, target)
|
||
self.isce_autorun(isce_exe_dir, isce_work_space)
|
||
|
||
# 1.3 转为tif
|
||
in_dict = ["incLocal.rdr"]
|
||
out_dict = ["a_incLocal.tif"]
|
||
self.isce_write(in_dict, out_dict)
|
||
|
||
# 1.4 删除文件夹
|
||
path_list = [self.__workspace_isce_path, self.__workspace_mintpy_path]
|
||
self.isce_del_dir(path_list)
|
||
logger.info('progress bar :35％')
|
||
|
||
# 二、 主影像
|
||
# 2.1 执行stackSentinel.exe 生成运行配置文件 #
|
||
# key_word = "20190113"
|
||
key_word = self.mas_key_word
|
||
self.isce_stackSentinel(slc_dir, dem_path, aux_dir, orbits_dir, isce_work_space, key_word, isce_exe_dir,box)
|
||
|
||
# 2.2 执行autorun.py生成干涉图
|
||
run_steps = ["run_16_unwrap"]
|
||
self.isce_run_steps(run_steps, target)
|
||
self.isce_autorun(isce_exe_dir, isce_work_space)
|
||
|
||
# 2.3 转为tif
|
||
in_dict = ["hgt.rdr", "incLocal.rdr", "filt_fine.int.vrt"] # 高程、局部入射角、干涉图
|
||
out_dict = ["m_hgt.tif", "m_incLocal.tif", "m_filt_fine.tif"]
|
||
self.isce_write(in_dict, out_dict)
|
||
|
||
incLocal_resame_path = os.path.join(self.__workspace_slc_path,
|
||
"resame_" + os.path.split(self.__processing_paras["a_incLocal"])[1])
|
||
a_proj, a_geotrans, a_arr = ImageHandler().read_img(self.__processing_paras["a_incLocal"])
|
||
m_proj, m_geotrans, m_arr = ImageHandler().read_img(self.__processing_paras["m_incLocal"])
|
||
if a_arr.shape != m_arr.shape:
|
||
logger.error('master inclocal.rdr is different from Auxiliary inclocal.rdr!')
|
||
ImageHandler().write_img(incLocal_resame_path, a_proj, a_geotrans, a_arr, no_data='null') # todo
|
||
self.__processing_paras["a_incLocal"] = incLocal_resame_path
|
||
logger.info('progress bar :75％')
|
||
return True
|
||
|
||
def tif_int(self):
|
||
"""GeoTif转成isce格式"""
|
||
img_path = self.atmospheric_product_path
|
||
AtmosDelay_path = os.path.join(self.__workspace_slc_path, "resamed_"+os.path.basename(img_path))
|
||
refer_img_path = os.path.join(self.__workspace_slc_path, "m_hgt.tif")
|
||
pp.resampling_by_scale(img_path, AtmosDelay_path, refer_img_path)
|
||
|
||
filename = os.listdir(os.path.join(self.__workspace_isce_path, "merged", "interferograms"))[0]
|
||
vrt_path = os.path.join(self.__workspace_isce_path, "merged", "interferograms", filename, "filt_fine.int.vrt")
|
||
|
||
# 数组存到.vrt中
|
||
Atmo_arr=DataTransform().Cover(AtmosDelay_path)
|
||
Atmo_arr = Atmo_arr.astype(np.complex64)
|
||
DataTransform().tran(vrt_path, Atmo_arr)
|
||
DataTransform().tran_af(vrt_path)
|
||
return True
|
||
|
||
def phase_unwrapping(self):
|
||
"""执行相位解缠"""
|
||
# 2.1移除run_files 中的steps
|
||
isce_work_space = r"/".join(self.__workspace_isce_path.split("\\"))
|
||
folder_run_files = os.path.join(self.__workspace_isce_path, "run_files")
|
||
run_steps = os.listdir(folder_run_files)
|
||
steps = run_steps[0:15]
|
||
self.isce_run_steps(steps, isce_work_space)
|
||
run_16_unwrap_file = os.path.join(self.__workspace_isce_path, "run_16_unwrap")
|
||
shutil.move(run_16_unwrap_file, folder_run_files)
|
||
|
||
# 2.2 执行autorun.py生成解缠图
|
||
isce_work_space = r"/".join(self.__workspace_isce_path.split("\\"))
|
||
isce_exe_dir = r"/".join(os.path.join(self.env_str, "ISCEApp").split("\\"))
|
||
self.isce_autorun(isce_exe_dir, isce_work_space)
|
||
|
||
# 2.3 解缠图转tif
|
||
in_dict = ["filt_fine.unw.vrt"]
|
||
out_dict = ["filt_fine_unw.tif"]
|
||
self.isce_write(in_dict, out_dict)
|
||
unw_file = os.path.join(self.__workspace_slc_path, "filt_fine_unw.tif")
|
||
unwrapped_file= os.path.join(self.__workspace_atmos_dely_path, "filt_fine_unw.tif")
|
||
shutil.copy(unw_file, unwrapped_file)
|
||
self.imageHandler.write_quick_view(unwrapped_file) # 输出的快视图1/3
|
||
logger.info('progress bar :95％')
|
||
return True
|
||
|
||
def creat_xml_format(self, start):
|
||
# 3.1生成元文件案例
|
||
model_xml_path = self.creat_xml(self.atmospheric_product_path, start)
|
||
SrcImageName_M = os.path.split(self.__input_paras["MasterSarData"]['ParaValue'])[1].split('.tar.gz')[0]
|
||
SrcImageName_A = os.path.split(self.__input_paras["AuxiliarySarData"]['ParaValue'])[1].split('.tar.gz')[0]
|
||
SrcImageName = SrcImageName_M+";"+SrcImageName_A
|
||
meta_xml_path = self.__workspace_atmos_dely_path+EXE_NAME+"Product.meta.xml"
|
||
self.image_meta_xml = None
|
||
CreateMetafile(self.image_meta_xml, self.alg_xml_path, model_xml_path, meta_xml_path).process(SrcImageName)
|
||
# 3.2生成压缩包
|
||
self.make_targz(self.__out_para, self.__workspace_atmos_dely_path)
|
||
logger.info('make targz finish')
|
||
# self.del_folder(self.__workspace_preprocessing_path)
|
||
# self.del_folder(self.__workspace_preprocessed_path)
|
||
# self.del_folder(self.__workspace_cut_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 = 'AtmosphericDelay.xml'
|
||
else:
|
||
xml_path = sys.argv[1]
|
||
AtmosphericMain = AtmosphericMain(xml_path)
|
||
if AtmosphericMain.check_source() is False:
|
||
raise Exception('check_source() failed!')
|
||
if AtmosphericMain.get_interv_img() is False:
|
||
raise Exception('get_interv_img() failed!')
|
||
if AtmosphericMain.preprocess_handle() is False:
|
||
raise Exception('preprocess_handle() failed!')
|
||
if AtmosphericMain.prod_handle() is False:
|
||
raise Exception('prod_handle() failed!')
|
||
if AtmosphericMain.tif_int() is False:
|
||
raise Exception('prod_handle() failed!')
|
||
if AtmosphericMain.phase_unwrapping() is False:
|
||
raise Exception('phase_unwrapping() failed!')
|
||
if AtmosphericMain.creat_xml_format(start) is not True:
|
||
raise Exception('creat_xml_format() failed!')
|
||
logger.info('successful production of ' + EXE_NAME + ' Product'+'.mate.xml')
|
||
except Exception:
|
||
logger.exception("run-time error!")
|
||
finally:
|
||
AtmosphericMain.del_temp_workspace()
|
||
# pass
|
||
|
||
end = datetime.datetime.now()
|
||
msg = 'running use time: %s ' % (end - start)
|
||
logger.info(msg)
|