From 9c7a96ace9540ed6d379a34152409542cf10a166 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr> Date: Thu, 30 Jan 2020 15:25:06 +0000 Subject: [PATCH] ENH : Refactoring for python_src directory --- python_src/diapOTB.py | 478 +++------ python_src/diapOTB_OLD.py | 532 ++++++++++ python_src/diapOTB_S1IW.py | 847 +++------------- python_src/diapOTB_S1IW_OLD.py | 905 ++++++++++++++++++ python_src/processings/DInSar.py | 553 +++++++++++ python_src/processings/Ground.py | 186 ++++ python_src/processings/Metadata_Correction.py | 93 ++ python_src/processings/Pre_Processing.py | 204 ++++ python_src/processings/__init__.py | 0 python_src/utils/DiapOTB_applications.py | 584 +++++++++++ python_src/utils/__init__.py | 0 python_src/utils/func_utils.py | 235 +++++ 12 files changed, 3551 insertions(+), 1066 deletions(-) create mode 100644 python_src/diapOTB_OLD.py create mode 100644 python_src/diapOTB_S1IW_OLD.py create mode 100644 python_src/processings/DInSar.py create mode 100644 python_src/processings/Ground.py create mode 100644 python_src/processings/Metadata_Correction.py create mode 100644 python_src/processings/Pre_Processing.py create mode 100644 python_src/processings/__init__.py create mode 100644 python_src/utils/DiapOTB_applications.py create mode 100644 python_src/utils/__init__.py create mode 100644 python_src/utils/func_utils.py diff --git a/python_src/diapOTB.py b/python_src/diapOTB.py index a39c35d..1f57666 100644 --- a/python_src/diapOTB.py +++ b/python_src/diapOTB.py @@ -20,8 +20,6 @@ # limitations under the License. # -from __future__ import print_function - __author__ = "POC-THALES" __version__ = "0.1" __status__ = "Developpement" @@ -31,48 +29,16 @@ __last_modified__ = "27/10/2017" # Imports import sys import logging -import json -from jsonschema import validate import os import argparse import h5py -import otbApplication as otb - - -# Streamer to our log file -class StreamToLogger(object): - """ - Fake file-like stream object that redirects writes to a logger instance. - """ - def __init__(self, logger, log_level=logging.INFO): - self.logger = logger - self.log_level = log_level - - def write(self, buf): - for line in buf.rstrip().splitlines(): - self.logger.log(self.log_level, line.rstrip()) - - def flush(self): - for handler in self.logger.handlers: - handler.flush() - - -def validate_json(json, schema): - try: - validate(json, schema) - except Exception as valid_err: - print("Invalid JSON: {err}".format(err=valid_err)) - return False - else: - # Realise votre travail - print("Valid JSON") - return True - -# string to bool -def str2bool(v): - return v.lower() in ("yes", "true", "t", "1") +from processings import Pre_Processing +from processings import Ground +from processings import DInSar +from processings import Metadata_Correction +from utils import func_utils # Main if __name__ == "__main__": @@ -83,50 +49,11 @@ if __name__ == "__main__": args = parser.parse_args() print(args.configfile) - # Logger initialization - logger = logging.getLogger(__name__) - logger.setLevel(logging.INFO) - - LogFormatter = logging.Formatter('%(filename)s :: %(levelname)s :: %(message)s') - - # Create console handler with a high log level (warning level) - stream_handler = logging.StreamHandler() - stream_handler.setLevel(logging.WARNING) + func_utils.init_logger() + + dataConfig = func_utils.load_configfile(args.configfile) - # Add Handlers - logger.addHandler(stream_handler) - - # Read and Load the configuration file - try: - with open(args.configfile, 'r') as f: - dataConfig = json.load(f) - - except Exception as err: - logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=args.configfile)) - quit() - - # Load schema (into DiapOTB install) - diapOTB_install = os.getenv('DIAPOTB_INSTALL') - if diapOTB_install is not None and os.path.exists(diapOTB_install): - schemas_path = os.path.join(diapOTB_install, "json_schemas") - if os.path.exists(schemas_path): - schema_S1SM = os.path.join(schemas_path, "schema_S1SM.json") - - try: - with open(schema_S1SM, "r") as sch: - dataSchema = json.load(sch) - except Exception as err: - logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=schema_S1SM)) - quit() - - # Check Json file - jsonIsValid = validate_json(dataConfig, dataSchema) - - if not jsonIsValid : - logger.critical("Error, provided config file does not match requirements") - quit() - # Get dictionaries dict_Global = dataConfig['Global'] dict_PreProcessing = dataConfig['Pre_Processing'] @@ -174,17 +101,17 @@ if __name__ == "__main__": gain_interfero = dict_DInSAR['parameter']['Interferogram_gain'] if (geoGrid_threshold < 0) or (geoGrid_threshold > 1) : - logger.critical("Error, Wrong Threshold for fine deformation grid") + func_utils.log(logging.CRITICAL, "Error, Wrong Threshold for fine deformation grid") # Check if images exist if not os.path.exists(master_Image) : - logger.critical("{img} does not exist. Check its path.".format(img=master_Image)) + func_utils.log(logging.CRITICAL, "Error, {img} does not exist. Check its path.".format(img=master_Image)) quit() if not os.path.exists(slave_Image) : - logger.critical("{img} does not exist. Check its path.".format(img=slave_Image)) + func_utils.log(logging.CRITICAL, "Error, {img} does not exist. Check its path.".format(img=slave_Image)) quit() if not os.path.exists(dem) : - logger.critical("{img} does not exist. Check its path.".format(img=dem)) + func_utils.log(logging.CRITICAL, "Error, {img} does not exist. Check its path.".format(img=dem)) quit() if not os.path.exists(output_dir): print("The output directory does not exist and will be created") @@ -193,58 +120,43 @@ if __name__ == "__main__": print("The output directory exists. Some files can be overwritten") - # File handler for the logger - # Create file handler which logs even info messages (used as stdout redirection) - file_handler = logging.FileHandler(os.path.join(output_dir, 'info.log'), 'a') - file_handler.setLevel(logging.INFO) - file_handler.setFormatter(LogFormatter) - - # Add Handlers - logger.addHandler(file_handler) - - # Redirect stdout and stderr to logger - s1 = StreamToLogger(logger, logging.INFO) - stdout_saveWrite = sys.stdout.write # Save stdout.write to print some info into the console - stdout_saveFlush = sys.stdout.flush # Save stdout.flush to print some info into the console - sys.stdout.write = s1.write # Replace stdout.write by our StreamToLogger - sys.stdout.flush = s1.flush # Replace stdout.flush by our StreamToLogger - stdout_save = s1 # Different object - stdout_save.write = stdout_saveWrite # Restore stdout.write into stdout_save - stdout_save.flush = stdout_saveFlush # Restore stdout.write into stdout_save + # Init file handler (all normaly print on std is redirected into info.log) + # To use previous print on std, use printOnStd + func_utils.init_fileLog(output_dir) # Recap of input parameter into info.log - logger.info("########### Input Parameters for the current execution ############## ") - logger.info(" Pre_Processing : ") - logger.info("ml_range : {param}".format(param=ml_range)) - logger.info("ml_azimut : {param}".format(param=ml_azimut)) - logger.info("ml_gain : {param}".format(param=ml_gain)) - logger.info("dop_file : {param}".format(param=dop_file)) + func_utils.log(logging.INFO, "########### Input Parameters for the current execution ############## ") + func_utils.log(logging.INFO, " Pre_Processing : ") + func_utils.log(logging.INFO, "ml_range : {param}".format(param=ml_range)) + func_utils.log(logging.INFO, "ml_azimut : {param}".format(param=ml_azimut)) + func_utils.log(logging.INFO, "ml_gain : {param}".format(param=ml_gain)) + func_utils.log(logging.INFO, "dop_file : {param}".format(param=dop_file)) # Metadata_Correction - logger.info(" Metadata_Correction : ") - logger.info("activateMetadataCorrection : {param}".format(param=activateMetadataCorrection)) + func_utils.log(logging.INFO, " Metadata_Correction : ") + func_utils.log(logging.INFO, "activateMetadataCorrection : {param}".format(param=activateMetadataCorrection)) if activateMetadataCorrection : - logger.info("ml_simu_range : {param}".format(param=ml_simu_range)) - logger.info("ml_simu_azimut : {param}".format(param=ml_simu_azimut)) - logger.info("ml_simu_gain : {param}".format(param=ml_simu_gain)) - logger.info("ml_correlSimu_range : {param}".format(param=ml_correlSimu_range)) - logger.info("ml_correlSimu_azimut : {param}".format(param=ml_correlSimu_azimut)) - logger.info("correlSimu_gridstep_range : {param}".format(param=correlSimu_gridstep_range)) - logger.info("correlSimu_gridstep_azimut : {param}".format(param=correlSimu_gridstep_azimut)) - logger.info("fine_metadata_file : {param}".format(param=fine_metadata_file)) + func_utils.log(logging.INFO, "ml_simu_range : {param}".format(param=ml_simu_range)) + func_utils.log(logging.INFO, "ml_simu_azimut : {param}".format(param=ml_simu_azimut)) + func_utils.log(logging.INFO, "ml_simu_gain : {param}".format(param=ml_simu_gain)) + func_utils.log(logging.INFO, "ml_correlSimu_range : {param}".format(param=ml_correlSimu_range)) + func_utils.log(logging.INFO, "ml_correlSimu_azimut : {param}".format(param=ml_correlSimu_azimut)) + func_utils.log(logging.INFO, "correlSimu_gridstep_range : {param}".format(param=correlSimu_gridstep_range)) + func_utils.log(logging.INFO, "correlSimu_gridstep_azimut : {param}".format(param=correlSimu_gridstep_azimut)) + func_utils.log(logging.INFO, "fine_metadata_file : {param}".format(param=fine_metadata_file)) # DIn_SAR - logger.info(" DIn_SAR : ") - logger.info("geoGrid_gridstep_range : {param}".format(param=geoGrid_gridstep_range)) - logger.info("geoGrid_gridstep_azimut : {param}".format(param=geoGrid_gridstep_azimut)) - logger.info("geoGrid_threshold : {param}".format(param=geoGrid_threshold)) - logger.info("geoGrid_gap : {param}".format(param=geoGrid_gap)) - logger.info("ml_geoGrid_range : {param}".format(param=ml_geoGrid_range)) - logger.info("ml_geoGrid_azimut : {param}".format(param=ml_geoGrid_azimut)) - logger.info("gain_interfero : {param}".format(param=gain_interfero)) + func_utils.log(logging.INFO, " DIn_SAR : ") + func_utils.log(logging.INFO, "geoGrid_gridstep_range : {param}".format(param=geoGrid_gridstep_range)) + func_utils.log(logging.INFO, "geoGrid_gridstep_azimut : {param}".format(param=geoGrid_gridstep_azimut)) + func_utils.log(logging.INFO, "geoGrid_threshold : {param}".format(param=geoGrid_threshold)) + func_utils.log(logging.INFO, "geoGrid_gap : {param}".format(param=geoGrid_gap)) + func_utils.log(logging.INFO, "ml_geoGrid_range : {param}".format(param=ml_geoGrid_range)) + func_utils.log(logging.INFO, "ml_geoGrid_azimut : {param}".format(param=ml_geoGrid_azimut)) + func_utils.log(logging.INFO, "gain_interfero : {param}".format(param=gain_interfero)) - logger.info("########### Input Images for the current execution ############## ") + func_utils.log(logging.INFO, "########### Input Images for the current execution ############## ") master_Image_base = os.path.basename(master_Image) slave_Image_base = os.path.basename(slave_Image) @@ -253,8 +165,8 @@ if __name__ == "__main__": master_ext = master_Image.split(".")[-1:] slave_ext = slave_Image.split(".")[-1:] - logger.info("master_ext = {ext}".format(ext=master_ext[0])) - logger.info("slave_ext = {ext}".format(ext=slave_ext[0])) + func_utils.log(logging.INFO, "master_ext = {ext}".format(ext=master_ext[0])) + func_utils.log(logging.INFO, "slave_ext = {ext}".format(ext=slave_ext[0])) if master_ext[0] == "h5" : @@ -263,17 +175,17 @@ if __name__ == "__main__": if len(lDataSet_master) != 1 : - logger.critical("Error, H5 input files does not contain the expected dataset") + func_utils.log(logging.CRITICAL, "Error, H5 input files does not contain the expected dataset") quit() if lDataSet_master[0] != "S01" : - logger.critical("Error, H5 input files does not contain the expected dataset") + func_utils.log(logging.CRITICAL, "Error, H5 input files does not contain the expected dataset") quit() master_S01 = dict(master_H5['S01']) if not 'SBI' in master_S01: - logger.critical("H5 input files does not contain the expected dataset") + func_utils.log(logging.CRITICAL, "Error, H5 input files does not contain the expected dataset") quit() # Change the name of master and slave image to read directly the //S01/SBI @@ -287,246 +199,124 @@ if __name__ == "__main__": lDataSet_slave = list(slave_H5.keys()) if len(lDataSet_slave) != 1 : - logger.critical("H5 input files does not contain the expected dataset") + func_utils.log(logging.CRITICAL, "Error, H5 input files does not contain the expected dataset") quit() if lDataSet_slave[0] != "S01" : - logger.critical("H5 input files does not contain the expected dataset") + func_utils.log(logging.CRITICAL, "Error, H5 input files does not contain the expected dataset") quit() slave_S01 = dict(slave_H5['S01']) if not 'SBI' in slave_S01 : - logger.critical("H5 input files does not contain the expected dataset") + func_utils.log(logging.CRITICAL, "Error, H5 input files does not contain the expected dataset") quit() slave_Image = "HDF5:" + slave_Image + "://S01/SBI" - logger.info("master_Image = {img}".format(img=master_Image)) - logger.info("slave_Image = {img}".format(img=slave_Image)) - logger.info("dem : {param}".format(param=dem)) + func_utils.log(logging.INFO, "master_Image = {img}".format(img=master_Image)) + func_utils.log(logging.INFO, "slave_Image = {img}".format(img=slave_Image)) + func_utils.log(logging.INFO, "dem : {param}".format(param=dem)) - print("\n Beginning of DiapOTB processing \n", file=stdout_save) - logger.info("############ Beginning of DiapOTB processing ##############") + + func_utils.printOnStd("\n Beginning of DiapOTB processing (S1 SM or Cosmo mode) \n") + func_utils.log(logging.INFO, "############ Beginning of DiapOTB processing (S1 SM or Cosmo mode) ##############") ####################### Pre Processing Chain ########################## - ######## SARDoppler Application ####### - print("\n Doppler Application \n", file=stdout_save) - logger.info("Doppler Application") # Master - dopFile = open(os.path.join(output_dir, dop_file), "w") - dopFile.write("Doppler for master image : " + os.path.basename(master_Image_base)+ "\n") - dopFile.close() - appDoppler0Master = otb.Registry.CreateApplication("SARDoppler0") - appDoppler0Master.SetParameterString("insar", master_Image) - appDoppler0Master.SetParameterString("outfile", os.path.join(output_dir, dop_file)) - appDoppler0Master.SetParameterString("ram", "4000") - appDoppler0Master.ExecuteAndWriteOutput() - + func_utils.printOnStd("\n Master Pre_Processing chain \n") + func_utils.log(logging.INFO, "Master Pre_Processing Application") + + paramPreMaster = {} + paramPreMaster['ml_range'] = ml_range + paramPreMaster['ml_azimut'] = ml_azimut + paramPreMaster['ml_gain'] = ml_gain + paramPreMaster['dop_file'] = dop_file + + dop0Master = Pre_Processing.extractToMultilook(master_Image, master_Image_base, paramPreMaster, "Others", + output_dir) # Slave - dopFile = open(os.path.join(output_dir, dop_file), "a") - dopFile.write("Doppler for slave image : " + os.path.basename(slave_Image_base) + "\n") - dopFile.close() - appDoppler0Slave = otb.Registry.CreateApplication("SARDoppler0") - appDoppler0Slave.SetParameterString("insar", slave_Image) - appDoppler0Slave.SetParameterString("outfile", os.path.join(output_dir, dop_file)) - appDoppler0Slave.SetParameterString("ram", "4000") - appDoppler0Slave.ExecuteAndWriteOutput() - - - ####### SARMultiLook Application ####### - print("\n MultiLook Application \n", file=stdout_save) - logger.info("MultiLook Application") - # Master - master_Image_ML = os.path.splitext(master_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" - appMultiLookMaster = otb.Registry.CreateApplication("SARMultiLook") - appMultiLookMaster.SetParameterString("incomplex", master_Image) - appMultiLookMaster.SetParameterString("out", os.path.join(output_dir, master_Image_ML)) - appMultiLookMaster.SetParameterInt("mlran",ml_range) - appMultiLookMaster.SetParameterInt("mlazi",ml_azimut) - appMultiLookMaster.SetParameterFloat("mlgain", ml_gain) - appMultiLookMaster.SetParameterString("ram", "4000") - appMultiLookMaster.ExecuteAndWriteOutput() + func_utils.printOnStd("\n Slave Pre_Processing chain \n") + func_utils.log(logging.INFO, "Slave Pre_Processing Application") - # Slave - slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" - appMultiLookSlave = otb.Registry.CreateApplication("SARMultiLook") - appMultiLookSlave.SetParameterString("incomplex", slave_Image) - appMultiLookSlave.SetParameterString("out", os.path.join(output_dir, slave_Image_ML)) - appMultiLookSlave.SetParameterInt("mlran",ml_range) - appMultiLookSlave.SetParameterInt("mlazi",ml_azimut) - appMultiLookSlave.SetParameterFloat("mlgain", ml_gain) - appMultiLookSlave.SetParameterString("ram", "4000") - appMultiLookSlave.ExecuteAndWriteOutput() + paramPreSlave = {} + paramPreSlave['ml_range'] = ml_range + paramPreSlave['ml_azimut'] = ml_azimut + paramPreSlave['ml_gain'] = ml_gain + paramPreSlave['dop_file'] = dop_file + dop0Slave = Pre_Processing.extractToMultilook(slave_Image, slave_Image_base, paramPreSlave, "Others", + output_dir) # ######################## Metadata Correction Chain ############################# - if activateMetadataCorrection : - ######## SARDEMToAmplitude Application (Simu_SAR step) ####### - print("\n SARDEMToAmplitude Application \n", file=stdout_save) - logger.info("SARDEMToAmplitude Application") - amplitude_simu_image = os.path.splitext(master_Image_base)[0] + "_simuSAR" + "_ml" + str(ml_simu_azimut) + str(ml_simu_range) + ".tif" - appDEMToAmplitude = otb.Registry.CreateApplication("SARDEMToAmplitude") - appDEMToAmplitude.SetParameterString("insar", master_Image) - appDEMToAmplitude.SetParameterString("indem", dem) - appDEMToAmplitude.SetParameterString("out", os.path.join(output_dir, amplitude_simu_image)) - appDEMToAmplitude.SetParameterInt("mlran",ml_simu_range) - appDEMToAmplitude.SetParameterInt("mlazi",ml_simu_azimut) - appDEMToAmplitude.SetParameterFloat("mlgain", ml_simu_gain) - appDEMToAmplitude.SetParameterInt("nodata", -32768) - appDEMToAmplitude.SetParameterString("ram", "4000") - #appDEMToAmplitude.Execute() - appDEMToAmplitude.ExecuteAndWriteOutput() - - - ######## SARCorrelationGrid Application (Correl step) ####### - print("\n SARCorrelationGrid Application \n", file=stdout_save) - logger.info("SARCorrelationGrid Application") - correl_grid = "correl_simu" + "_gridstep" + str(correlSimu_gridstep_azimut) + str(correlSimu_gridstep_range) + ".tif" - appCorGrid = otb.Registry.CreateApplication("SARCorrelationGrid") - appCorGrid.SetParameterString("inmaster", os.path.join(output_dir, master_Image_ML)) - appCorGrid.SetParameterString("inslave", os.path.join(output_dir, amplitude_simu_image)) - #appCorGrid.SetParameterInputImage("inslave", appDEMToAmplitude.GetParameterOutputImage("out")) # Input image - appCorGrid.SetParameterString("out", os.path.join(output_dir, correl_grid)) - appCorGrid.SetParameterInt("mlran",ml_correlSimu_range) - appCorGrid.SetParameterInt("mlazi",ml_correlSimu_azimut) - appCorGrid.SetParameterInt("gridsteprange", correlSimu_gridstep_range) - appCorGrid.SetParameterInt("gridstepazimut", correlSimu_gridstep_azimut) - appCorGrid.SetParameterString("ram", "4000") - #appCorGrid.Execute() - appCorGrid.ExecuteAndWriteOutput() - - - ######## SARFineMetadata Application (Correct_snrt step) ####### - print("\n SARFineMetadata Application \n", file=stdout_save) - logger.info("SARFineMetadata Application") - appFineMetadata = otb.Registry.CreateApplication("SARFineMetadata") - appFineMetadata.SetParameterString("insar", master_Image) - appFineMetadata.SetParameterString("ingrid", os.path.join(output_dir, correl_grid)) - #appFineMetadata.SetParameterInputImage("ingrid", appCorGrid.GetParameterOutputImage("out")) # Input image - appFineMetadata.SetParameterFloat("stepmax", 0.1) - appFineMetadata.SetParameterFloat("threshold", 0.3) - appFineMetadata.SetParameterString("outfile", os.path.join(output_dir, fine_metadata_file)) - appFineMetadata.ExecuteAndWriteOutput() - + if activateMetadataCorrection : + # TO DO + func_utils.printOnStd("Metadata Correction Chain not available for now") + # paramMetadata = {} + # paramMetadata['ml_range'] = ml_simu_range + # paramMetadata['ml_azimut'] = ml_simu_azimut + # paramMetadata['ml_gain'] = ml_simu_gain + # paramMetadata['geoGrid_gridstep_range'] = correlSimu_gridstep_range + # paramMetadata['geoGrid_gridstep_azimut'] = correlSimu_gridstep_azimut + # paramMetadata['nodata'] = -32768 + # paramMetadata['fine_metadata_file'] = fine_metadata_file + # Metadata_Correction.fineMetadata(master_Image, master_Image_base, dem, paramMetadata, output_dir) - - ######################## DIn_SAR Chain ############################# - ######## SARDEMProjection Application ####### - print("\n SARDEMProjection Application \n", file=stdout_save) - logger.info("SARDEMProjection Application") + + ######################### Ground Chain ############################# # Master - demProj_Master = "demProj_Master.tif" - appDEMProjectionMaster = otb.Registry.CreateApplication("SARDEMProjection") - appDEMProjectionMaster.SetParameterString("insar", master_Image) - appDEMProjectionMaster.SetParameterString("indem", dem) - if activateMetadataCorrection : - appDEMProjectionMaster.SetParameterString("infilemetadata", os.path.join(output_dir, fine_metadata_file)) - appDEMProjectionMaster.SetParameterString("withxyz", "true") - appDEMProjectionMaster.SetParameterInt("nodata", -32768) - appDEMProjectionMaster.SetParameterString("out", os.path.join(output_dir, demProj_Master)) - appDEMProjectionMaster.SetParameterString("ram", "4000") - appDEMProjectionMaster.ExecuteAndWriteOutput() + func_utils.printOnStd("\n Master Ground chain \n") + func_utils.log(logging.INFO, "Master Ground Application") + + paramGroundMaster = {} + paramGroundMaster['nodata'] = -32768 + paramGroundMaster['withxyz'] = "true" + paramGroundMaster['for_slave_Image'] = False + + Ground.demProjectionAndCartesianEstimation(master_Image, master_Image_base, dem, paramGroundMaster, + "Others", output_dir) + # Slave - demProj_Slave = "demProj_Slave.tif" - appDEMProjectionSlave = otb.Registry.CreateApplication("SARDEMProjection") - appDEMProjectionSlave.SetParameterString("insar", slave_Image) - appDEMProjectionSlave.SetParameterString("indem", dem) - appDEMProjectionSlave.SetParameterString("withxyz", "true"); - appDEMProjectionSlave.SetParameterInt("nodata", -32768) - appDEMProjectionSlave.SetParameterString("out", os.path.join(output_dir, demProj_Slave)) - appDEMProjectionSlave.SetParameterString("ram", "4000") - appDEMProjectionSlave.ExecuteAndWriteOutput() - - - ######## SARFineDeformationGrid Application (geo_grid step) ####### - print("\n SARFineDeformationGrid Application \n", file=stdout_save) - logger.info("SARFineDeformationGrid Application") - fine_grid = "fineDeformationGrid.tif" - appFineDeformationGrid = otb.Registry.CreateApplication("SARFineDeformationGrid") - appFineDeformationGrid.SetParameterString("indem", dem) - appFineDeformationGrid.SetParameterString("insarmaster", master_Image) - appFineDeformationGrid.SetParameterString("insarslave", slave_Image) - appFineDeformationGrid.SetParameterString("inmlmaster", os.path.join(output_dir, master_Image_ML)) - appFineDeformationGrid.SetParameterString("inmlslave", os.path.join(output_dir, slave_Image_ML)) - appFineDeformationGrid.SetParameterString("indemprojmaster", os.path.join(output_dir, demProj_Master)) - appFineDeformationGrid.SetParameterString("indemprojslave", os.path.join(output_dir, demProj_Slave)) - appFineDeformationGrid.SetParameterInt("mlran", ml_geoGrid_range) - appFineDeformationGrid.SetParameterInt("mlazi", ml_geoGrid_azimut) - appFineDeformationGrid.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appFineDeformationGrid.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold) - appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap) - appFineDeformationGrid.SetParameterString("out", os.path.join(output_dir, fine_grid)) - # Projection is not reliable for cosmo sensor => correction of all values with correlation grid - if satellite == "cosmo" or satellite == "CSK" : - appFineDeformationGrid.SetParameterString("advantage", "correlation") - appFineDeformationGrid.SetParameterString("ram", "4000") - appFineDeformationGrid.ExecuteAndWriteOutput() + func_utils.printOnStd("\n Slave Ground chain \n") + func_utils.log(logging.INFO, "Slave Ground Application") + paramGroundSlave = {} + paramGroundSlave['nodata'] = -32768 + paramGroundSlave['withxyz'] = "true" + paramGroundSlave['for_slave_Image'] = True - ######## SARCoRegistration Application (changeo step) ####### - print("\n SARCoRegistration Application \n", file=stdout_save) - logger.info("SARCoRegistration Application") - slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_coregistrated.tif" - appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") - appCoRegistration.SetParameterString("insarmaster", master_Image) - appCoRegistration.SetParameterString("insarslave", slave_Image) - appCoRegistration.SetParameterString("ingrid", os.path.join(output_dir, fine_grid)) - appCoRegistration.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appCoRegistration.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appCoRegistration.SetParameterFloat("doppler0", appDoppler0Slave.GetParameterFloat("doppler0")) - appCoRegistration.SetParameterInt("sizetiles", 50) - appCoRegistration.SetParameterInt("margin", 7) - appCoRegistration.SetParameterInt("nbramps", 257) - appCoRegistration.SetParameterString("ram", "4000") - appCoRegistration.SetParameterString("out", os.path.join(output_dir, slave_Image_CoRe)) - appCoRegistration.ExecuteAndWriteOutput() + Ground.demProjectionAndCartesianEstimation(slave_Image, slave_Image_base, dem, paramGroundSlave, + "Others", output_dir) + - - ######## SARCartesianMeanEstimation Application ####### - print("\n SARCartesianMeanEstimation Application \n") - logger.info("SARCartesianMeanEstimation Application") - # Master - master_cartesian_mean = "CartMeanMaster.tif" - master_cartesianperline_mean = "CartMeanPerLineMaster.tif" - appCartMeanMaster = otb.Registry.CreateApplication("SARCartesianMeanEstimation") - appCartMeanMaster.SetParameterString("insar", master_Image) - appCartMeanMaster.SetParameterString("indem", dem) - appCartMeanMaster.SetParameterString("indemproj", os.path.join(output_dir, demProj_Master)) - appCartMeanMaster.SetParameterInt("indirectiondemc", appDEMProjectionMaster.GetParameterInt("directiontoscandemc")) - appCartMeanMaster.SetParameterInt("indirectiondeml", appDEMProjectionMaster.GetParameterInt("directiontoscandeml")) - appCartMeanMaster.SetParameterInt("mlran", 1) - appCartMeanMaster.SetParameterInt("mlazi", 1) - appCartMeanMaster.SetParameterString("ram", "4000") - appCartMeanMaster.SetParameterString("out", os.path.join(output_dir, master_cartesian_mean)) - appCartMeanMaster.ExecuteAndWriteOutput() - - - ######## SARRobustInterferogram Application (interf step) ####### - print("\n SARRobustInterferogram Application \n", file=stdout_save) - logger.info("SARRobustInterferogram Application") - interferogram = "interferogram.tif" - appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") - appInterferogram.SetParameterString("insarmaster", master_Image) - appInterferogram.SetParameterString("incoregistratedslave", os.path.join(output_dir, slave_Image_CoRe)) - appInterferogram.SetParameterString("insarslave", slave_Image) - appInterferogram.SetParameterString("ingrid", os.path.join(output_dir, fine_grid)) - appInterferogram.SetParameterString("incartmeanmaster", os.path.join(output_dir, master_cartesian_mean)) - appInterferogram.SetParameterInt("mlran", ml_range) - appInterferogram.SetParameterInt("mlazi", ml_azimut) - appInterferogram.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appInterferogram.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appInterferogram.SetParameterFloat("gain", gain_interfero) - appInterferogram.SetParameterString("ram", "4000") - appInterferogram.SetParameterString("out", os.path.join(output_dir, interferogram)) - appInterferogram.ExecuteAndWriteOutput() - - - print("\n End of DiapOTB processing \n", file=stdout_save) - logger.info("############# End of DiapOTB processing ##############") + ######################## DIn_SAR Chain ############################# + func_utils.printOnStd("\n DIn_SAR chain \n") + func_utils.log(logging.INFO, "DIn_SAR chain") + + advantage = "projection" + if satellite == "cosmo" or satellite == "CSK" : + advantage = "correlation" + + # Create param + param = {} + param['ml_azimut'] = ml_azimut + param['ml_range'] = ml_range + param['ml_geoGrid_azimut'] = ml_geoGrid_azimut + param['ml_geoGrid_range'] = ml_geoGrid_range + param['geoGrid_gridstep_range'] = geoGrid_gridstep_range + param['geoGrid_gridstep_azimut'] = geoGrid_gridstep_azimut + param['geoGrid_threshold'] = geoGrid_threshold + param['geoGrid_gap'] = geoGrid_gap + param['doppler0'] = dop0Slave + param['gain_interfero'] = gain_interfero + param['advantage'] = advantage + + DInSar.gridToInterferogram(dem, master_Image, master_Image_base, slave_Image, slave_Image_base, output_dir, output_dir, param, "Others", output_dir) + + func_utils.printOnStd("\n End of DiapOTB processing (S1 SM or Cosmo mode) \n") + func_utils.log(logging.INFO, "############ End of DiapOTB processing (S1 SM or Cosmo mode) ##############") diff --git a/python_src/diapOTB_OLD.py b/python_src/diapOTB_OLD.py new file mode 100644 index 0000000..a39c35d --- /dev/null +++ b/python_src/diapOTB_OLD.py @@ -0,0 +1,532 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +from __future__ import print_function + +__author__ = "POC-THALES" +__version__ = "0.1" +__status__ = "Developpement" +__date__ = "27/10/2017" +__last_modified__ = "27/10/2017" + +# Imports +import sys +import logging +import json +from jsonschema import validate +import os +import argparse +import h5py + +import otbApplication as otb + + +# Streamer to our log file +class StreamToLogger(object): + """ + Fake file-like stream object that redirects writes to a logger instance. + """ + def __init__(self, logger, log_level=logging.INFO): + self.logger = logger + self.log_level = log_level + + def write(self, buf): + for line in buf.rstrip().splitlines(): + self.logger.log(self.log_level, line.rstrip()) + + def flush(self): + for handler in self.logger.handlers: + handler.flush() + + +def validate_json(json, schema): + try: + validate(json, schema) + except Exception as valid_err: + print("Invalid JSON: {err}".format(err=valid_err)) + return False + else: + # Realise votre travail + print("Valid JSON") + return True + +# string to bool +def str2bool(v): + return v.lower() in ("yes", "true", "t", "1") + + +# Main +if __name__ == "__main__": + + # Check arguments + parser = argparse.ArgumentParser() + parser.add_argument("configfile", help="input conguration file for the application DiapOTB") + args = parser.parse_args() + print(args.configfile) + + # Logger initialization + logger = logging.getLogger(__name__) + logger.setLevel(logging.INFO) + + LogFormatter = logging.Formatter('%(filename)s :: %(levelname)s :: %(message)s') + + # Create console handler with a high log level (warning level) + stream_handler = logging.StreamHandler() + stream_handler.setLevel(logging.WARNING) + + # Add Handlers + logger.addHandler(stream_handler) + + + # Read and Load the configuration file + try: + with open(args.configfile, 'r') as f: + dataConfig = json.load(f) + + except Exception as err: + logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=args.configfile)) + quit() + + # Load schema (into DiapOTB install) + diapOTB_install = os.getenv('DIAPOTB_INSTALL') + if diapOTB_install is not None and os.path.exists(diapOTB_install): + schemas_path = os.path.join(diapOTB_install, "json_schemas") + if os.path.exists(schemas_path): + schema_S1SM = os.path.join(schemas_path, "schema_S1SM.json") + + try: + with open(schema_S1SM, "r") as sch: + dataSchema = json.load(sch) + except Exception as err: + logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=schema_S1SM)) + quit() + + # Check Json file + jsonIsValid = validate_json(dataConfig, dataSchema) + + if not jsonIsValid : + logger.critical("Error, provided config file does not match requirements") + quit() + + # Get dictionaries + dict_Global = dataConfig['Global'] + dict_PreProcessing = dataConfig['Pre_Processing'] + dict_Metadata_Correction = dataConfig['Metadata_Correction'] + dict_DInSAR = dataConfig['DIn_SAR'] + + # Get elements from configuration file + # Global + master_Image = dict_Global['in']['Master_Image_Path'] + slave_Image = dict_Global['in']['Slave_Image_Path'] + dem = dict_Global['in']['DEM_Path'] + output_dir = dict_Global['out']['output_dir'] + + satellite = "default" + mode = "default" + + if 'sensor' in dict_Global: + satellite = dict_Global['sensor']['satellite'] + mode = dict_Global['sensor']['mode'] + + # Pre_Processing + ml_range = dict_PreProcessing['parameter']['ML_range'] + ml_azimut = dict_PreProcessing['parameter']['ML_azimut'] + ml_gain = dict_PreProcessing['parameter']['ML_gain'] + dop_file = dict_PreProcessing['out']['doppler_file'] + + # Metadata_Correction + activateMetadataCorrection = dict_Metadata_Correction['parameter']['activate'] + ml_simu_range = ml_range + ml_simu_azimut = ml_azimut + ml_simu_gain = 1. + ml_correlSimu_range = ml_range + ml_correlSimu_azimut = ml_azimut + correlSimu_gridstep_range = dict_Metadata_Correction['parameter']['GridStep_range'] + correlSimu_gridstep_azimut = dict_Metadata_Correction['parameter']['GridStep_azimut'] + fine_metadata_file = dict_Metadata_Correction['out']['fine_metadata_file'] + + # DIn_SAR + geoGrid_gridstep_range = dict_DInSAR['parameter']['GridStep_range'] + geoGrid_gridstep_azimut = dict_DInSAR['parameter']['GridStep_azimut'] + geoGrid_threshold = dict_DInSAR['parameter']['Grid_Threshold'] + geoGrid_gap = dict_DInSAR['parameter']['Grid_Gap'] + ml_geoGrid_range = ml_range + ml_geoGrid_azimut = ml_azimut + gain_interfero = dict_DInSAR['parameter']['Interferogram_gain'] + + if (geoGrid_threshold < 0) or (geoGrid_threshold > 1) : + logger.critical("Error, Wrong Threshold for fine deformation grid") + + # Check if images exist + if not os.path.exists(master_Image) : + logger.critical("{img} does not exist. Check its path.".format(img=master_Image)) + quit() + if not os.path.exists(slave_Image) : + logger.critical("{img} does not exist. Check its path.".format(img=slave_Image)) + quit() + if not os.path.exists(dem) : + logger.critical("{img} does not exist. Check its path.".format(img=dem)) + quit() + if not os.path.exists(output_dir): + print("The output directory does not exist and will be created") + os.makedirs(output_dir) + else : + print("The output directory exists. Some files can be overwritten") + + + # File handler for the logger + # Create file handler which logs even info messages (used as stdout redirection) + file_handler = logging.FileHandler(os.path.join(output_dir, 'info.log'), 'a') + file_handler.setLevel(logging.INFO) + file_handler.setFormatter(LogFormatter) + + # Add Handlers + logger.addHandler(file_handler) + + # Redirect stdout and stderr to logger + s1 = StreamToLogger(logger, logging.INFO) + stdout_saveWrite = sys.stdout.write # Save stdout.write to print some info into the console + stdout_saveFlush = sys.stdout.flush # Save stdout.flush to print some info into the console + sys.stdout.write = s1.write # Replace stdout.write by our StreamToLogger + sys.stdout.flush = s1.flush # Replace stdout.flush by our StreamToLogger + stdout_save = s1 # Different object + stdout_save.write = stdout_saveWrite # Restore stdout.write into stdout_save + stdout_save.flush = stdout_saveFlush # Restore stdout.write into stdout_save + + # Recap of input parameter into info.log + logger.info("########### Input Parameters for the current execution ############## ") + logger.info(" Pre_Processing : ") + logger.info("ml_range : {param}".format(param=ml_range)) + logger.info("ml_azimut : {param}".format(param=ml_azimut)) + logger.info("ml_gain : {param}".format(param=ml_gain)) + logger.info("dop_file : {param}".format(param=dop_file)) + + # Metadata_Correction + logger.info(" Metadata_Correction : ") + logger.info("activateMetadataCorrection : {param}".format(param=activateMetadataCorrection)) + if activateMetadataCorrection : + logger.info("ml_simu_range : {param}".format(param=ml_simu_range)) + logger.info("ml_simu_azimut : {param}".format(param=ml_simu_azimut)) + logger.info("ml_simu_gain : {param}".format(param=ml_simu_gain)) + logger.info("ml_correlSimu_range : {param}".format(param=ml_correlSimu_range)) + logger.info("ml_correlSimu_azimut : {param}".format(param=ml_correlSimu_azimut)) + logger.info("correlSimu_gridstep_range : {param}".format(param=correlSimu_gridstep_range)) + logger.info("correlSimu_gridstep_azimut : {param}".format(param=correlSimu_gridstep_azimut)) + logger.info("fine_metadata_file : {param}".format(param=fine_metadata_file)) + + # DIn_SAR + logger.info(" DIn_SAR : ") + logger.info("geoGrid_gridstep_range : {param}".format(param=geoGrid_gridstep_range)) + logger.info("geoGrid_gridstep_azimut : {param}".format(param=geoGrid_gridstep_azimut)) + logger.info("geoGrid_threshold : {param}".format(param=geoGrid_threshold)) + logger.info("geoGrid_gap : {param}".format(param=geoGrid_gap)) + logger.info("ml_geoGrid_range : {param}".format(param=ml_geoGrid_range)) + logger.info("ml_geoGrid_azimut : {param}".format(param=ml_geoGrid_azimut)) + logger.info("gain_interfero : {param}".format(param=gain_interfero)) + + + logger.info("########### Input Images for the current execution ############## ") + + master_Image_base = os.path.basename(master_Image) + slave_Image_base = os.path.basename(slave_Image) + + # Check extension (if .h5 => HDF5 file => Cosmo Sensor) + master_ext = master_Image.split(".")[-1:] + slave_ext = slave_Image.split(".")[-1:] + + logger.info("master_ext = {ext}".format(ext=master_ext[0])) + logger.info("slave_ext = {ext}".format(ext=slave_ext[0])) + + + if master_ext[0] == "h5" : + master_H5 = h5py.File(master_Image, 'r') + lDataSet_master = list(master_H5.keys()) + + + if len(lDataSet_master) != 1 : + logger.critical("Error, H5 input files does not contain the expected dataset") + quit() + + if lDataSet_master[0] != "S01" : + logger.critical("Error, H5 input files does not contain the expected dataset") + quit() + + master_S01 = dict(master_H5['S01']) + + if not 'SBI' in master_S01: + logger.critical("H5 input files does not contain the expected dataset") + quit() + + # Change the name of master and slave image to read directly the //S01/SBI + master_Image = "HDF5:" + master_Image + "://S01/SBI" + # Adapt sattelite + satellite = "cosmo" + + + if slave_ext[0] == "h5" : + slave_H5 = h5py.File(slave_Image, 'r') + lDataSet_slave = list(slave_H5.keys()) + + if len(lDataSet_slave) != 1 : + logger.critical("H5 input files does not contain the expected dataset") + quit() + + if lDataSet_slave[0] != "S01" : + logger.critical("H5 input files does not contain the expected dataset") + quit() + + slave_S01 = dict(slave_H5['S01']) + + if not 'SBI' in slave_S01 : + logger.critical("H5 input files does not contain the expected dataset") + quit() + + slave_Image = "HDF5:" + slave_Image + "://S01/SBI" + + logger.info("master_Image = {img}".format(img=master_Image)) + logger.info("slave_Image = {img}".format(img=slave_Image)) + logger.info("dem : {param}".format(param=dem)) + + print("\n Beginning of DiapOTB processing \n", file=stdout_save) + logger.info("############ Beginning of DiapOTB processing ##############") + + ####################### Pre Processing Chain ########################## + ######## SARDoppler Application ####### + print("\n Doppler Application \n", file=stdout_save) + logger.info("Doppler Application") + # Master + dopFile = open(os.path.join(output_dir, dop_file), "w") + dopFile.write("Doppler for master image : " + os.path.basename(master_Image_base)+ "\n") + dopFile.close() + appDoppler0Master = otb.Registry.CreateApplication("SARDoppler0") + appDoppler0Master.SetParameterString("insar", master_Image) + appDoppler0Master.SetParameterString("outfile", os.path.join(output_dir, dop_file)) + appDoppler0Master.SetParameterString("ram", "4000") + appDoppler0Master.ExecuteAndWriteOutput() + + # Slave + dopFile = open(os.path.join(output_dir, dop_file), "a") + dopFile.write("Doppler for slave image : " + os.path.basename(slave_Image_base) + "\n") + dopFile.close() + appDoppler0Slave = otb.Registry.CreateApplication("SARDoppler0") + appDoppler0Slave.SetParameterString("insar", slave_Image) + appDoppler0Slave.SetParameterString("outfile", os.path.join(output_dir, dop_file)) + appDoppler0Slave.SetParameterString("ram", "4000") + appDoppler0Slave.ExecuteAndWriteOutput() + + + ####### SARMultiLook Application ####### + print("\n MultiLook Application \n", file=stdout_save) + logger.info("MultiLook Application") + # Master + master_Image_ML = os.path.splitext(master_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + appMultiLookMaster = otb.Registry.CreateApplication("SARMultiLook") + appMultiLookMaster.SetParameterString("incomplex", master_Image) + appMultiLookMaster.SetParameterString("out", os.path.join(output_dir, master_Image_ML)) + appMultiLookMaster.SetParameterInt("mlran",ml_range) + appMultiLookMaster.SetParameterInt("mlazi",ml_azimut) + appMultiLookMaster.SetParameterFloat("mlgain", ml_gain) + appMultiLookMaster.SetParameterString("ram", "4000") + appMultiLookMaster.ExecuteAndWriteOutput() + + # Slave + slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + appMultiLookSlave = otb.Registry.CreateApplication("SARMultiLook") + appMultiLookSlave.SetParameterString("incomplex", slave_Image) + appMultiLookSlave.SetParameterString("out", os.path.join(output_dir, slave_Image_ML)) + appMultiLookSlave.SetParameterInt("mlran",ml_range) + appMultiLookSlave.SetParameterInt("mlazi",ml_azimut) + appMultiLookSlave.SetParameterFloat("mlgain", ml_gain) + appMultiLookSlave.SetParameterString("ram", "4000") + appMultiLookSlave.ExecuteAndWriteOutput() + + + + + # ######################## Metadata Correction Chain ############################# + if activateMetadataCorrection : + ######## SARDEMToAmplitude Application (Simu_SAR step) ####### + print("\n SARDEMToAmplitude Application \n", file=stdout_save) + logger.info("SARDEMToAmplitude Application") + amplitude_simu_image = os.path.splitext(master_Image_base)[0] + "_simuSAR" + "_ml" + str(ml_simu_azimut) + str(ml_simu_range) + ".tif" + appDEMToAmplitude = otb.Registry.CreateApplication("SARDEMToAmplitude") + appDEMToAmplitude.SetParameterString("insar", master_Image) + appDEMToAmplitude.SetParameterString("indem", dem) + appDEMToAmplitude.SetParameterString("out", os.path.join(output_dir, amplitude_simu_image)) + appDEMToAmplitude.SetParameterInt("mlran",ml_simu_range) + appDEMToAmplitude.SetParameterInt("mlazi",ml_simu_azimut) + appDEMToAmplitude.SetParameterFloat("mlgain", ml_simu_gain) + appDEMToAmplitude.SetParameterInt("nodata", -32768) + appDEMToAmplitude.SetParameterString("ram", "4000") + #appDEMToAmplitude.Execute() + appDEMToAmplitude.ExecuteAndWriteOutput() + + + ######## SARCorrelationGrid Application (Correl step) ####### + print("\n SARCorrelationGrid Application \n", file=stdout_save) + logger.info("SARCorrelationGrid Application") + correl_grid = "correl_simu" + "_gridstep" + str(correlSimu_gridstep_azimut) + str(correlSimu_gridstep_range) + ".tif" + appCorGrid = otb.Registry.CreateApplication("SARCorrelationGrid") + appCorGrid.SetParameterString("inmaster", os.path.join(output_dir, master_Image_ML)) + appCorGrid.SetParameterString("inslave", os.path.join(output_dir, amplitude_simu_image)) + #appCorGrid.SetParameterInputImage("inslave", appDEMToAmplitude.GetParameterOutputImage("out")) # Input image + appCorGrid.SetParameterString("out", os.path.join(output_dir, correl_grid)) + appCorGrid.SetParameterInt("mlran",ml_correlSimu_range) + appCorGrid.SetParameterInt("mlazi",ml_correlSimu_azimut) + appCorGrid.SetParameterInt("gridsteprange", correlSimu_gridstep_range) + appCorGrid.SetParameterInt("gridstepazimut", correlSimu_gridstep_azimut) + appCorGrid.SetParameterString("ram", "4000") + #appCorGrid.Execute() + appCorGrid.ExecuteAndWriteOutput() + + + ######## SARFineMetadata Application (Correct_snrt step) ####### + print("\n SARFineMetadata Application \n", file=stdout_save) + logger.info("SARFineMetadata Application") + appFineMetadata = otb.Registry.CreateApplication("SARFineMetadata") + appFineMetadata.SetParameterString("insar", master_Image) + appFineMetadata.SetParameterString("ingrid", os.path.join(output_dir, correl_grid)) + #appFineMetadata.SetParameterInputImage("ingrid", appCorGrid.GetParameterOutputImage("out")) # Input image + appFineMetadata.SetParameterFloat("stepmax", 0.1) + appFineMetadata.SetParameterFloat("threshold", 0.3) + appFineMetadata.SetParameterString("outfile", os.path.join(output_dir, fine_metadata_file)) + appFineMetadata.ExecuteAndWriteOutput() + + + + + + ######################## DIn_SAR Chain ############################# + ######## SARDEMProjection Application ####### + print("\n SARDEMProjection Application \n", file=stdout_save) + logger.info("SARDEMProjection Application") + # Master + demProj_Master = "demProj_Master.tif" + appDEMProjectionMaster = otb.Registry.CreateApplication("SARDEMProjection") + appDEMProjectionMaster.SetParameterString("insar", master_Image) + appDEMProjectionMaster.SetParameterString("indem", dem) + if activateMetadataCorrection : + appDEMProjectionMaster.SetParameterString("infilemetadata", os.path.join(output_dir, fine_metadata_file)) + appDEMProjectionMaster.SetParameterString("withxyz", "true") + appDEMProjectionMaster.SetParameterInt("nodata", -32768) + appDEMProjectionMaster.SetParameterString("out", os.path.join(output_dir, demProj_Master)) + appDEMProjectionMaster.SetParameterString("ram", "4000") + appDEMProjectionMaster.ExecuteAndWriteOutput() + + # Slave + demProj_Slave = "demProj_Slave.tif" + appDEMProjectionSlave = otb.Registry.CreateApplication("SARDEMProjection") + appDEMProjectionSlave.SetParameterString("insar", slave_Image) + appDEMProjectionSlave.SetParameterString("indem", dem) + appDEMProjectionSlave.SetParameterString("withxyz", "true"); + appDEMProjectionSlave.SetParameterInt("nodata", -32768) + appDEMProjectionSlave.SetParameterString("out", os.path.join(output_dir, demProj_Slave)) + appDEMProjectionSlave.SetParameterString("ram", "4000") + appDEMProjectionSlave.ExecuteAndWriteOutput() + + + ######## SARFineDeformationGrid Application (geo_grid step) ####### + print("\n SARFineDeformationGrid Application \n", file=stdout_save) + logger.info("SARFineDeformationGrid Application") + fine_grid = "fineDeformationGrid.tif" + appFineDeformationGrid = otb.Registry.CreateApplication("SARFineDeformationGrid") + appFineDeformationGrid.SetParameterString("indem", dem) + appFineDeformationGrid.SetParameterString("insarmaster", master_Image) + appFineDeformationGrid.SetParameterString("insarslave", slave_Image) + appFineDeformationGrid.SetParameterString("inmlmaster", os.path.join(output_dir, master_Image_ML)) + appFineDeformationGrid.SetParameterString("inmlslave", os.path.join(output_dir, slave_Image_ML)) + appFineDeformationGrid.SetParameterString("indemprojmaster", os.path.join(output_dir, demProj_Master)) + appFineDeformationGrid.SetParameterString("indemprojslave", os.path.join(output_dir, demProj_Slave)) + appFineDeformationGrid.SetParameterInt("mlran", ml_geoGrid_range) + appFineDeformationGrid.SetParameterInt("mlazi", ml_geoGrid_azimut) + appFineDeformationGrid.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appFineDeformationGrid.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold) + appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap) + appFineDeformationGrid.SetParameterString("out", os.path.join(output_dir, fine_grid)) + # Projection is not reliable for cosmo sensor => correction of all values with correlation grid + if satellite == "cosmo" or satellite == "CSK" : + appFineDeformationGrid.SetParameterString("advantage", "correlation") + appFineDeformationGrid.SetParameterString("ram", "4000") + appFineDeformationGrid.ExecuteAndWriteOutput() + + + ######## SARCoRegistration Application (changeo step) ####### + print("\n SARCoRegistration Application \n", file=stdout_save) + logger.info("SARCoRegistration Application") + slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_coregistrated.tif" + appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") + appCoRegistration.SetParameterString("insarmaster", master_Image) + appCoRegistration.SetParameterString("insarslave", slave_Image) + appCoRegistration.SetParameterString("ingrid", os.path.join(output_dir, fine_grid)) + appCoRegistration.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appCoRegistration.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appCoRegistration.SetParameterFloat("doppler0", appDoppler0Slave.GetParameterFloat("doppler0")) + appCoRegistration.SetParameterInt("sizetiles", 50) + appCoRegistration.SetParameterInt("margin", 7) + appCoRegistration.SetParameterInt("nbramps", 257) + appCoRegistration.SetParameterString("ram", "4000") + appCoRegistration.SetParameterString("out", os.path.join(output_dir, slave_Image_CoRe)) + appCoRegistration.ExecuteAndWriteOutput() + + + ######## SARCartesianMeanEstimation Application ####### + print("\n SARCartesianMeanEstimation Application \n") + logger.info("SARCartesianMeanEstimation Application") + # Master + master_cartesian_mean = "CartMeanMaster.tif" + master_cartesianperline_mean = "CartMeanPerLineMaster.tif" + appCartMeanMaster = otb.Registry.CreateApplication("SARCartesianMeanEstimation") + appCartMeanMaster.SetParameterString("insar", master_Image) + appCartMeanMaster.SetParameterString("indem", dem) + appCartMeanMaster.SetParameterString("indemproj", os.path.join(output_dir, demProj_Master)) + appCartMeanMaster.SetParameterInt("indirectiondemc", appDEMProjectionMaster.GetParameterInt("directiontoscandemc")) + appCartMeanMaster.SetParameterInt("indirectiondeml", appDEMProjectionMaster.GetParameterInt("directiontoscandeml")) + appCartMeanMaster.SetParameterInt("mlran", 1) + appCartMeanMaster.SetParameterInt("mlazi", 1) + appCartMeanMaster.SetParameterString("ram", "4000") + appCartMeanMaster.SetParameterString("out", os.path.join(output_dir, master_cartesian_mean)) + appCartMeanMaster.ExecuteAndWriteOutput() + + + ######## SARRobustInterferogram Application (interf step) ####### + print("\n SARRobustInterferogram Application \n", file=stdout_save) + logger.info("SARRobustInterferogram Application") + interferogram = "interferogram.tif" + appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") + appInterferogram.SetParameterString("insarmaster", master_Image) + appInterferogram.SetParameterString("incoregistratedslave", os.path.join(output_dir, slave_Image_CoRe)) + appInterferogram.SetParameterString("insarslave", slave_Image) + appInterferogram.SetParameterString("ingrid", os.path.join(output_dir, fine_grid)) + appInterferogram.SetParameterString("incartmeanmaster", os.path.join(output_dir, master_cartesian_mean)) + appInterferogram.SetParameterInt("mlran", ml_range) + appInterferogram.SetParameterInt("mlazi", ml_azimut) + appInterferogram.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appInterferogram.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appInterferogram.SetParameterFloat("gain", gain_interfero) + appInterferogram.SetParameterString("ram", "4000") + appInterferogram.SetParameterString("out", os.path.join(output_dir, interferogram)) + appInterferogram.ExecuteAndWriteOutput() + + + print("\n End of DiapOTB processing \n", file=stdout_save) + logger.info("############# End of DiapOTB processing ##############") diff --git a/python_src/diapOTB_S1IW.py b/python_src/diapOTB_S1IW.py index d5a8364..65c3b3d 100644 --- a/python_src/diapOTB_S1IW.py +++ b/python_src/diapOTB_S1IW.py @@ -20,8 +20,6 @@ # limitations under the License. # -from __future__ import print_function - __author__ = "POC-THALES" __version__ = "0.1" __status__ = "Developpement" @@ -30,87 +28,16 @@ __last_modified__ = "11/12/2018" # Imports import logging -import json -from jsonschema import validate import os import sys import argparse import re -import otbApplication as otb - -# Streamer to our log file -class StreamToLogger(object): - """ - Fake file-like stream object that redirects writes to a logger instance. - """ - def __init__(self, logger, log_level=logging.INFO): - self.logger = logger - self.log_level = log_level - - def write(self, buf): - for line in buf.rstrip().splitlines(): - self.logger.log(self.log_level, line.rstrip()) - - def flush(self): - for handler in self.logger.handlers: - handler.flush() - - -def validate_json(json, schema): - try: - validate(json, schema) - except Exception as valid_err: - print("Invalid JSON: {}".format(valid_err)) - return False - else: - # Realise votre travail - print("Valid JSON") - return True - -# string to bool -def str2bool(v): - return v.lower() in ("yes", "true", "t", "1") - -# Bursts selection (according to anx time values) -def selectBurst(dictMaster, dictSlave, firstBurst, lastBurst, nbBurstSlave, validBurstMaster, validBurstSlave): - - key1Burst = "support_data.geom.bursts.burst[" - - # Initialize the output lists (empty lists) - validBurstMaster.clear() - validBurstSlave.clear() - - # Loop on Master bursts - for id_B in range(firstBurst, lastBurst+1): - keyBurstMaster = key1Burst + str(id_B) + "].azimuth_anx_time" - - # Get the anx time for the current burst (into Master Image) - anxMaster = float(dictMaster[keyBurstMaster]) - - # Loop on slave bursts to find the closest anx time - minDiff = 200 - id_B_save = id_B - for id_B_slave in range(0, nbBurstSlave): - keyBurstSlave = key1Burst + str(id_B_slave) + "].azimuth_anx_time" - - # Get anx time for slave burst - anxSlave = float(dictSlave[keyBurstSlave]) - - # Comparaison between master and slave - diff = abs(anxMaster - anxSlave) - - if minDiff > diff : - minDiff = diff - id_B_save = id_B_slave - - # Check if difference between the anx time is valid (must be inferior to 1) - if minDiff < 1. : - # Fill lists with master Burst_id and the selected slave burst_id - validBurstMaster.append(id_B) - validBurstSlave.append(id_B_save) - +from processings import Pre_Processing +from processings import Ground +from processings import DInSar +from utils import func_utils # Main if __name__ == "__main__": @@ -121,50 +48,11 @@ if __name__ == "__main__": args = parser.parse_args() print(args.configfile) - # Logger initialization - logger = logging.getLogger(__name__) - logger.setLevel(logging.INFO) - - LogFormatter = logging.Formatter('%(filename)s :: %(levelname)s :: %(message)s') - - # Create console handler with a high log level (warning level) - stream_handler = logging.StreamHandler() - stream_handler.setLevel(logging.WARNING) - - # Add Handlers - logger.addHandler(stream_handler) - - # Read and Load the configuration file - try: - with open(args.configfile, 'r') as f: - dataConfig = json.load(f) - - except Exception as err: - logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=args.configfile)) - quit() + func_utils.init_logger() + + dataConfig = func_utils.load_configfile(args.configfile, "S1_IW") - # Load schema (into DiapOTB install) - diapOTB_install = os.getenv('DIAPOTB_INSTALL') - if diapOTB_install is not None and os.path.exists(diapOTB_install): - schemas_path = os.path.join(diapOTB_install, "json_schemas") - if os.path.exists(schemas_path): - schema_S1IW = os.path.join(schemas_path, "schema_S1IW.json") - - try: - with open(schema_S1IW, "r") as sch: - dataSchema = json.load(sch) - except Exception as err: - logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=schema_S1SM)) - quit() - - # Check Json file - jsonIsValid = validate_json(dataConfig, dataSchema) - - if not jsonIsValid : - logger.critical("Error, provided config file does not match requirements") - quit() - - # Get dictionaries + # Get dictionaries dict_Global = dataConfig['Global'] dict_PreProcessing = dataConfig['Pre_Processing'] dict_Ground = dataConfig['Ground'] @@ -210,18 +98,18 @@ if __name__ == "__main__": # Check Threshold if (geoGrid_threshold < 0) or (geoGrid_threshold > 1) : - logger.critical("Error, Wrong Threshold for fine deformation grid") + func_utils.log(logging.CRITICAL, "Error, Wrong Threshold for fine deformation grid") geoGrid_threshold = 0.3 # Check if images/dem exist if not os.path.exists(master_Image) : - logger.critical("{img} does not exist. Check its path.".format(img=master_Image)) + func_utils.log(logging.CRITICAL, "Error, {img} does not exist. Check its path.".format(img=master_Image)) quit() if not os.path.exists(slave_Image) : - logger.critical("{img} does not exist. Check its path.".format(img=slave_Image)) + func_utils.log(logging.CRITICAL, "Error, {img} does not exist. Check its path.".format(img=slave_Image)) quit() if not os.path.exists(dem) : - logger.critical("{img} does not exist. Check its path.".format(img=dem)) + func_utils.log(logging.CRITICAL, "Error, {img} does not exist. Check its path.".format(img=dem)) quit() if not os.path.exists(output_dir): print("The output directory does not exist and will be created") @@ -229,101 +117,59 @@ if __name__ == "__main__": else : print("The output directory exists. Some files can be overwritten") - - # File handler for the logger - # Create file handler which logs even info messages (used as stdout redirection) - file_handler = logging.FileHandler(os.path.join(output_dir, 'info.log'), 'a') - file_handler.setLevel(logging.INFO) - file_handler.setFormatter(LogFormatter) - - # Add Handlers - logger.addHandler(file_handler) - - # Redirect stdout and stderr to logger - s1 = StreamToLogger(logger, logging.INFO) - stdout_saveWrite = sys.stdout.write # Save stdout.write to print some info into the console - stdout_saveFlush = sys.stdout.flush # Save stdout.flush to print some info into the console - sys.stdout.write = s1.write # Replace stdout.write by our StreamToLogger - sys.stdout.flush = s1.flush # Replace stdout.flush by our StreamToLogger - stdout_save = s1 # Different object - stdout_save.write = stdout_saveWrite # Restore stdout.write into stdout_save - stdout_save.flush = stdout_saveFlush # Restore stdout.write into stdout_save - - # Recap of input parameter into info.log - logger.info("########### Input Parameters for the current execution ############## ") - logger.info(" Pre_Processing : ") - logger.info("ml_range : {param}".format(param=ml_range)) - logger.info("ml_azimut : {param}".format(param=ml_azimut)) - logger.info("ml_gain : {param}".format(param=ml_gain)) - logger.info("dop_file : {param}".format(param=dop_file)) + + # Init file handler (all normaly print on std is redirected into info.log) + # To use previous print on std, use printOnStd + func_utils.init_fileLog(output_dir) + + # Recap of input parameter into info.log + func_utils.log(logging.INFO, "########### Input Parameters for the current execution ############## ") + func_utils.log(logging.INFO, " Pre_Processing : ") + func_utils.log(logging.INFO, "ml_range : {param}".format(param=ml_range)) + func_utils.log(logging.INFO, "ml_azimut : {param}".format(param=ml_azimut)) + func_utils.log(logging.INFO, "ml_gain : {param}".format(param=ml_gain)) + func_utils.log(logging.INFO, "dop_file : {param}".format(param=dop_file)) # DIn_SAR - logger.info(" DIn_SAR : ") - logger.info("geoGrid_gridstep_range : {param}".format(param=geoGrid_gridstep_range)) - logger.info("geoGrid_gridstep_azimut : {param}".format(param=geoGrid_gridstep_azimut)) - logger.info("geoGrid_threshold : {param}".format(param=geoGrid_threshold)) - logger.info("geoGrid_gap : {param}".format(param=geoGrid_gap)) - logger.info("ml_geoGrid_range : {param}".format(param=ml_geoGrid_range)) - logger.info("ml_geoGrid_azimut : {param}".format(param=ml_geoGrid_azimut)) - logger.info("gain_interfero : {param}".format(param=gain_interfero)) - logger.info("esd_AutoMode : {param}".format(param=esd_AutoMode)) - logger.info("esd_NbIter : {param}".format(param=esd_NbIter)) - - logger.info("########### Input Images for the current execution ############## ") - logger.info("master_Image : {param}".format(param=master_Image)) - logger.info("slave_Image : {param}".format(param=slave_Image)) - logger.info("dem : {param}".format(param=dem)) + func_utils.log(logging.INFO, " DIn_SAR : ") + func_utils.log(logging.INFO, "geoGrid_gridstep_range : {param}".format(param=geoGrid_gridstep_range)) + func_utils.log(logging.INFO, "geoGrid_gridstep_azimut : {param}".format(param=geoGrid_gridstep_azimut)) + func_utils.log(logging.INFO, "geoGrid_threshold : {param}".format(param=geoGrid_threshold)) + func_utils.log(logging.INFO, "geoGrid_gap : {param}".format(param=geoGrid_gap)) + func_utils.log(logging.INFO, "ml_geoGrid_range : {param}".format(param=ml_geoGrid_range)) + func_utils.log(logging.INFO, "ml_geoGrid_azimut : {param}".format(param=ml_geoGrid_azimut)) + func_utils.log(logging.INFO, "gain_interfero : {param}".format(param=gain_interfero)) + func_utils.log(logging.INFO, "esd_AutoMode : {param}".format(param=esd_AutoMode)) + func_utils.log(logging.INFO, "esd_NbIter : {param}".format(param=esd_NbIter)) + + func_utils.log(logging.INFO, "########### Input Images for the current execution ############## ") + func_utils.log(logging.INFO, "master_Image : {param}".format(param=master_Image)) + func_utils.log(logging.INFO, "slave_Image : {param}".format(param=slave_Image)) + func_utils.log(logging.INFO, "dem : {param}".format(param=dem)) # check Burst index master_Image_base = os.path.basename(master_Image) slave_Image_base = os.path.basename(slave_Image) # Retrieve some information about our input images - ReadImageInfo = otb.Registry.CreateApplication("ReadImageInfo") - ReadImageInfo.SetParameterString("in", master_Image) - ReadImageInfo.SetParameterString("keywordlist", "true") - ReadImageInfo.Execute() - - keywordMaster = ReadImageInfo.GetParameterString("keyword") - keywordlistMaster = ReadImageInfo.GetParameterString("keyword").split("\n") - keywordlistMaster = filter(None, keywordlistMaster) - dictKWLMaster = { i.split(':')[0] : re.sub(r"[\n\t\s]*", "", i.split(':')[1]) for i in keywordlistMaster } - - ReadImageInfo.SetParameterString("in", slave_Image) - ReadImageInfo.SetParameterString("keywordlist", "true") - ReadImageInfo.Execute() - - keywordSlave = ReadImageInfo.GetParameterString("keyword") - keywordlistSlave = ReadImageInfo.GetParameterString("keyword").split("\n") - keywordlistSlave = filter(None, keywordlistSlave) - dictKWLSlave = { i.split(':')[0] : re.sub(r"[\n\t\s]*", "", i.split(':')[1]) for i in keywordlistSlave } - + dictKWLMaster = func_utils.getImageKWL(master_Image) + dictKWLSlave = func_utils.getImageKWL(slave_Image) + # Check header version if int(dictKWLMaster['header.version']) < 3 or int(dictKWLSlave['header.version']) < 3 : - logger.critical("Upgrade your geom file") + func_utils.log(logging.CRITICAL, "Error, Upgrade your geom file") quit() # Get information about DEM (spacing, size ..) - ReadDEMInfo = otb.Registry.CreateApplication("ReadImageInfo") - ReadDEMInfo.SetParameterString("in", dem) - ReadDEMInfo.SetParameterString("keywordlist", "true") - ReadDEMInfo.Execute() - - spacingXDEM = ReadDEMInfo.GetParameterFloat("spacingx") - estimatedGroundSpacingXDEM = ReadDEMInfo.GetParameterFloat("estimatedgroundspacingx") - spacingYDEM = ReadDEMInfo.GetParameterFloat("spacingy") - estimatedGroundSpacingYDEM = ReadDEMInfo.GetParameterFloat("estimatedgroundspacingy") + dictDEMInformation = func_utils.getDEMInformation(dem) # Choose advantage for correlation or projection according to DEM resolution advantage = "projection" # By default projection - if estimatedGroundSpacingXDEM > 40. or estimatedGroundSpacingYDEM > 40. : + if dictDEMInformation['estimatedGroundSpacingXDEM'] > 40. or dictDEMInformation['estimatedGroundSpacingYDEM'] > 40. : advantage = "correlation" # Correlation if resolution > 40 m - logger.warning("Resolution of the input DEM is inferior to 40 meters : A correlation will be used to correct all deformation grids") - + func_utils.log(logging.WARNING, "Resolution of the input DEM is inferior to 40 meters : A correlation will be used to correct all deformation grids") - # Selection of bursts - keyBurst = "support_data.geom.bursts.burst[" + str(0) + "].azimuth_anx_time" # Check the index of bursts minNbBurst = min([int(dictKWLMaster['support_data.geom.bursts.number']), int(dictKWLSlave['support_data.geom.bursts.number'])]) @@ -342,564 +188,121 @@ if __name__ == "__main__": firstBurst = min(burstList) lastBurst = max(burstList) except Exception as err: - logger.critical("Wrong burst index") + func_utils.log(logging.CRITICAL, "Error, Wrong burst index") quit() if minNbBurst < firstBurst or minNbBurst < lastBurst or lastBurst < 0 or firstBurst < 0 : - logger.critical("Wrong burst index") + func_utils.log(logging.CRITICAL,"Error, Wrong burst index") quit() - validBurstMaster = [] - validBurstSlave = [] + nbBurstSlave = int(dictKWLSlave['support_data.geom.bursts.number']) - selectBurst(dictKWLMaster, dictKWLSlave, firstBurst, lastBurst, nbBurstSlave, validBurstMaster, validBurstSlave) + validBurstMaster, validBurstSlave = func_utils.selectBurst(dictKWLMaster, dictKWLSlave, firstBurst, lastBurst, nbBurstSlave) + if len(validBurstMaster) == 0 or len(validBurstSlave) == 0 : - logger.critical("Wrong burst index (slave index does not match with master index)") + func_utils.log(logging.CRITICAL, "Error, Wrong burst index (slave index does not match with master index)") quit() - # Update firstBurst and lastBurst with selected Burst for master image - firstBurst = validBurstMaster[0] - lastBurst = validBurstMaster[len(validBurstMaster)-1] - # Create directory for each burst for burstId in range(validBurstMaster[0], validBurstMaster[len(validBurstMaster)-1]+1): if not os.path.exists(os.path.join(output_dir, "burst" + str(burstId))): os.makedirs(os.path.join(output_dir, "burst" + str(burstId))) - print("\n Beginning of DiapOTB processing (S1 IW mode) \n", file=stdout_save) - logger.info("############ Beginning of DiapOTB processing (S1 IW mode) ##############") + func_utils.printOnStd("\n Beginning of DiapOTB processing (S1 IW mode) \n") + func_utils.log(logging.INFO, "############ Beginning of DiapOTB processing (S1 IW mode) ##############") ####################### Pre Processing Chain ########################## - # Initialisation of doppler file - dopFile = open(os.path.join(output_dir, dop_file), "w") - dopFile.close() - # Master - print("\n Master Pre-Processing \n", file=stdout_save) - logger.info("Master Pre-Processing") - - dop0Master = [] + func_utils.printOnStd("\n Master Pre_Processing chain \n") + func_utils.log(logging.INFO, "Master Pre_Processing Application") + + paramPreMaster = {} + paramPreMaster['ml_range'] = ml_range + paramPreMaster['ml_azimut'] = ml_azimut + paramPreMaster['ml_gain'] = ml_gain + paramPreMaster['dop_file'] = dop_file + paramPreMaster['validBurstMaster'] = validBurstMaster + + dop0Master = Pre_Processing.extractToMultilook(master_Image, master_Image_base, paramPreMaster, "S1_IW", + output_dir) - for id_loop in range(0, len(validBurstMaster)): - #burstId = id_loop + firstBurst - burstId = validBurstMaster[id_loop] - - print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) - logger.info("BurstId = {ind}".format(ind=burstId)) - - burst_dir = os.path.join(output_dir, "burst" + str(burstId)) - - ######## SARBurstExtraction Application ####### - print("\n Burst Extraction Application \n", file=stdout_save) - logger.info("Burst Extraction Application") - burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" - appBurstExtractionMaster = otb.Registry.CreateApplication("SARBurstExtraction") - appBurstExtractionMaster.SetParameterString("in", master_Image) - appBurstExtractionMaster.SetParameterString("out", os.path.join(burst_dir, burstM)) - appBurstExtractionMaster.SetParameterInt("burstindex", burstId) - appBurstExtractionMaster.SetParameterString("allpixels", "true") - appBurstExtractionMaster.SetParameterString("ram", "4000") - appBurstExtractionMaster.ExecuteAndWriteOutput() - - ######## SARDeramp Application ####### - print("\n Deramping Application \n", file=stdout_save) - logger.info("Deramping Application") - burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" - appDerampMaster = otb.Registry.CreateApplication("SARDeramp") - appDerampMaster.SetParameterString("in", os.path.join(burst_dir, burstM)) - appDerampMaster.SetParameterString("out", os.path.join(burst_dir, burstDerampM)) - appDerampMaster.SetParameterString("ram", "4000") - appDerampMaster.ExecuteAndWriteOutput() - - ######## SARDoppler Application ####### - print("\n Doppler Application \n", file=stdout_save) - logger.info("Doppler Application") - dopFile = open(os.path.join(output_dir, dop_file), "a") - dopFile.write("Doppler for master image : " + os.path.basename(master_Image_base)+ " for burst_index = " + str(burstId) + "\n") - dopFile.close() - appDoppler0Master = otb.Registry.CreateApplication("SARDoppler0") - appDoppler0Master.SetParameterString("insar", os.path.join(burst_dir, burstDerampM)) - appDoppler0Master.SetParameterString("outfile", os.path.join(output_dir, dop_file)) - appDoppler0Master.SetParameterString("ram", "4000") - appDoppler0Master.ExecuteAndWriteOutput() - - dop0Master.append(appDoppler0Master.GetParameterFloat('doppler0')) - - ####### SARMultiLook Application ####### - print("\n MultiLook Application \n", file=stdout_save) - logger.info("MultiLook Application") - master_Image_ML = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" - appMultiLookMaster = otb.Registry.CreateApplication("SARMultiLook") - appMultiLookMaster.SetParameterString("incomplex", os.path.join(burst_dir, burstDerampM)) - appMultiLookMaster.SetParameterString("out", os.path.join(burst_dir, master_Image_ML)) - appMultiLookMaster.SetParameterInt("mlran",ml_range) - appMultiLookMaster.SetParameterInt("mlazi",ml_azimut) - appMultiLookMaster.SetParameterFloat("mlgain", ml_gain) - appMultiLookMaster.SetParameterString("ram", "4000") - appMultiLookMaster.ExecuteAndWriteOutput() - # Slave - print("\n Slave Pre-Processing \n", file=stdout_save) - logger.info("Slave Pre-Processing") - - dop0Slave = [] - - for id_loop in range(0, len(validBurstMaster)): - #burstId = id_loop + firstBurst - burstId = validBurstMaster[id_loop] - burstId_slave = validBurstSlave[id_loop] - - print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) - logger.info("BurstId = {ind}".format(ind=burstId)) - - - burst_dir = os.path.join(output_dir, "burst" + str(burstId)) - - ######## SARBurstExtraction Application ####### - print("\n Burst Extraction Application \n", file=stdout_save) - logger.info("Burst Extraction Application") - burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" - appBurstExtractionSlave = otb.Registry.CreateApplication("SARBurstExtraction") - appBurstExtractionSlave.SetParameterString("in", slave_Image) - appBurstExtractionSlave.SetParameterString("out", os.path.join(burst_dir, burstS)) - appBurstExtractionSlave.SetParameterInt("burstindex", burstId_slave) - appBurstExtractionSlave.SetParameterString("allpixels", "true") - appBurstExtractionSlave.SetParameterString("ram", "4000") - appBurstExtractionSlave.ExecuteAndWriteOutput() + func_utils.printOnStd("\n Slave Pre_Processing chain \n") + func_utils.log(logging.INFO, "Slave Pre_Processing Application") + + paramPreSlave = {} + paramPreSlave['ml_range'] = ml_range + paramPreSlave['ml_azimut'] = ml_azimut + paramPreSlave['ml_gain'] = ml_gain + paramPreSlave['dop_file'] = dop_file + paramPreSlave['validBurstMaster'] = validBurstMaster + paramPreSlave['validBurstSlave'] = validBurstSlave + + dop0Slave = Pre_Processing.extractToMultilook(slave_Image, slave_Image_base, paramPreSlave, "S1_IW", + output_dir) - ######## SARDeramp Application ####### - print("\n Deramping Application \n", file=stdout_save) - logger.info("Deramping Application") - burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" - appDerampSlave = otb.Registry.CreateApplication("SARDeramp") - appDerampSlave.SetParameterString("in", os.path.join(burst_dir, burstS)) - appDerampSlave.SetParameterString("out", os.path.join(burst_dir, burstDerampS)) - appDerampSlave.SetParameterString("ram", "4000") - appDerampSlave.ExecuteAndWriteOutput() - - ######## SARDoppler Application ####### - print("\n Doppler Application \n", file=stdout_save) - logger.info("Doppler Application") - dopFile = open(os.path.join(output_dir, dop_file), "a") - dopFile.write("Doppler for slave image : " + os.path.basename(slave_Image_base) + " for burst_index = " + str(burstId_slave) + "\n") - dopFile.close() - appDoppler0Slave = otb.Registry.CreateApplication("SARDoppler0") - appDoppler0Slave.SetParameterString("insar", os.path.join(burst_dir, burstDerampS)) - appDoppler0Slave.SetParameterString("outfile", os.path.join(output_dir, dop_file)) - appDoppler0Slave.SetParameterString("ram", "4000") - appDoppler0Slave.ExecuteAndWriteOutput() - - dop0Slave.append(appDoppler0Slave.GetParameterFloat('doppler0')) - - ####### SARMultiLook Application ####### - print("\n MultiLook Application \n", file=stdout_save) - logger.info("MultiLook Application") - slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" - appMultiLookSlave = otb.Registry.CreateApplication("SARMultiLook") - appMultiLookSlave.SetParameterString("incomplex", os.path.join(burst_dir, burstDerampS)) - appMultiLookSlave.SetParameterString("out", os.path.join(burst_dir, slave_Image_ML)) - appMultiLookSlave.SetParameterInt("mlran",ml_range) - appMultiLookSlave.SetParameterInt("mlazi",ml_azimut) - appMultiLookSlave.SetParameterFloat("mlgain", ml_gain) - appMultiLookSlave.SetParameterString("ram", "4000") - appMultiLookSlave.ExecuteAndWriteOutput() - ######################### Ground Chain ############################# # Master - print("\n Master Ground chain \n", file=stdout_save) - logger.info("Master Ground Application") - - gainMaster = [] - directionDEMMaster = [] + func_utils.printOnStd("\n Master Ground chain \n") + func_utils.log(logging.INFO, "Master Ground Application") + + paramGroundMaster = {} + paramGroundMaster['nodata'] = -32768 + paramGroundMaster['withxyz'] = "true" + paramGroundMaster['validBurstMaster'] = validBurstMaster - for id_loop in range(0, len(validBurstMaster)): - #burstId = id_loop + firstBurst - burstId = validBurstMaster[id_loop] - - print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) - logger.info("BurstId = {ind}".format(ind=burstId)) - - burst_dir = os.path.join(output_dir, "burst" + str(burstId)) - - burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" - + Ground.demProjectionAndCartesianEstimation(master_Image, master_Image_base, dem, paramGroundMaster, + "S1_IW", output_dir) - ######## SARDEMProjection Application ####### - print("\n SARDEMProjection Application \n", file=stdout_save) - logger.info("SARDEMProjection Application") - demProj_Master = "demProj" + "_burst" + str(burstId) +"_Master.tif" - appDEMProjectionMaster = otb.Registry.CreateApplication("SARDEMProjection") - appDEMProjectionMaster.SetParameterString("insar", os.path.join(burst_dir, burstDerampM)) - appDEMProjectionMaster.SetParameterString("indem", dem) - appDEMProjectionMaster.SetParameterString("withxyz", "true") - appDEMProjectionMaster.SetParameterInt("nodata", -32768) - appDEMProjectionMaster.SetParameterString("out", os.path.join(burst_dir, demProj_Master)) - appDEMProjectionMaster.SetParameterString("ram", "4000") - appDEMProjectionMaster.ExecuteAndWriteOutput() - - gainMaster.append(appDEMProjectionMaster.GetParameterFloat('gain')) - directionDEMMaster.append([appDEMProjectionMaster.GetParameterInt('directiontoscandemc'), appDEMProjectionMaster.GetParameterInt('directiontoscandeml')]) - - ######### SARCartesianMeanEstimation Application ####### - print("\n SARCartesianMeanEstimation Application \n", file=stdout_save) - logger.info("SARCartesianMeanEstimation Application") - master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" - appCartMeanMaster = otb.Registry.CreateApplication("SARCartesianMeanEstimation") - appCartMeanMaster.SetParameterString("insar", os.path.join(burst_dir, burstDerampM)) - appCartMeanMaster.SetParameterString("indem", dem) - appCartMeanMaster.SetParameterString("indemproj", os.path.join(burst_dir, demProj_Master)) - appCartMeanMaster.SetParameterInt("indirectiondemc", appDEMProjectionMaster.GetParameterInt("directiontoscandemc")) - appCartMeanMaster.SetParameterInt("indirectiondeml", appDEMProjectionMaster.GetParameterInt("directiontoscandeml")) - appCartMeanMaster.SetParameterInt("mlran", 1) - appCartMeanMaster.SetParameterInt("mlazi", 1) - appCartMeanMaster.SetParameterString("ram", "4000") - appCartMeanMaster.SetParameterString("out", os.path.join(burst_dir, master_cartesian_mean)) - appCartMeanMaster.ExecuteAndWriteOutput() - - - print (directionDEMMaster) # Slave - print("\n Slave Ground chain \n", file=stdout_save) - logger.info("Slave Ground Application") - - gainSlave = [] - directionDEMSlave = [] - - for id_loop in range(0, len(validBurstMaster)): - #burstId = id_loop + firstBurst - burstId = validBurstMaster[id_loop] - burstId_slave = validBurstSlave[id_loop] - - print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) - logger.info("BurstId = {ind}".format(ind=burstId)) - - burst_dir = os.path.join(output_dir, "burst" + str(burstId)) - - burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" - - print("\n SARDEMProjection Application \n", file=stdout_save) - logger.info("SARDEMProjection Application") - demProj_Slave = "demProj" + "_burst" + str(burstId_slave) + "_Slave.tif" - appDEMProjectionSlave = otb.Registry.CreateApplication("SARDEMProjection") - appDEMProjectionSlave.SetParameterString("insar", os.path.join(burst_dir, burstDerampS)) - appDEMProjectionSlave.SetParameterString("indem", dem) - appDEMProjectionSlave.SetParameterString("withxyz", "true"); - appDEMProjectionSlave.SetParameterInt("nodata", -32768) - appDEMProjectionSlave.SetParameterString("out", os.path.join(burst_dir, demProj_Slave)) - appDEMProjectionSlave.SetParameterString("ram", "4000") - appDEMProjectionSlave.ExecuteAndWriteOutput() - + func_utils.printOnStd("\n Slave Ground chain \n") + func_utils.log(logging.INFO, "Slave Ground Application") + + paramGroundSlave = {} + paramGroundSlave['nodata'] = -32768 + paramGroundSlave['withxyz'] = "true" + paramGroundSlave['validBurstMaster'] = validBurstMaster + paramGroundSlave['validBurstSlave'] = validBurstSlave + Ground.demProjectionAndCartesianEstimation(slave_Image, slave_Image_base, dem, paramGroundSlave, + "S1_IW", output_dir) + + ######################## DIn_SAR Chain ############################# - print("\n DIn_SAR chain \n", file=stdout_save) - logger.info("DIn_SAR chain") + func_utils.printOnStd("\n DIn_SAR chain \n") + func_utils.log(logging.INFO, "DIn_SAR chain") counter = 0 list_of_Interferograms = [] - list_of_Grids = [] - - for id_loop in range(0, len(validBurstMaster)): - #burstId = id_loop + firstBurst - burstId = validBurstMaster[id_loop] - burstId_slave = validBurstSlave[id_loop] - - print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) - logger.info("BurstId = {ind}".format(ind=burstId)) - - burst_dir = os.path.join(output_dir, "burst" + str(burstId)) - - burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" - burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" - - burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" - burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" - - master_Image_ML = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" - slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" - - - demProj_Master = "demProj" + "_burst" + str(burstId) +"_Master.tif" - demProj_Slave = "demProj" + "_burst" + str(burstId_slave) + "_Slave.tif" - master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" - - ######## Step 1 : SARFineDeformationGrid ####### - ######## SARFineDeformationGrid Application (geo_grid step) ####### - print("\n SARFineDeformationGrid Application \n", file=stdout_save) - logger.info("SARFineDeformationGrid Application") - fine_grid = "fineDeformationGrid"+ "_burst" + str(burstId) + ".tif" - appFineDeformationGrid = otb.Registry.CreateApplication("SARFineDeformationGrid") - appFineDeformationGrid.SetParameterString("indem", dem) - appFineDeformationGrid.SetParameterString("insarmaster", os.path.join(burst_dir, burstDerampM)) - appFineDeformationGrid.SetParameterString("insarslave", os.path.join(burst_dir, burstDerampS)) - appFineDeformationGrid.SetParameterString("inmlmaster", os.path.join(burst_dir, master_Image_ML)) - appFineDeformationGrid.SetParameterString("inmlslave", os.path.join(burst_dir, slave_Image_ML)) - appFineDeformationGrid.SetParameterString("indemprojmaster", os.path.join(burst_dir, demProj_Master)) - appFineDeformationGrid.SetParameterString("indemprojslave", os.path.join(burst_dir, demProj_Slave)) - appFineDeformationGrid.SetParameterInt("mlran", ml_geoGrid_range) - appFineDeformationGrid.SetParameterInt("mlazi", ml_geoGrid_azimut) - appFineDeformationGrid.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appFineDeformationGrid.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold) - appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap) - appFineDeformationGrid.SetParameterString("advantage", advantage) - appFineDeformationGrid.SetParameterString("out", os.path.join(burst_dir, fine_grid)) - appFineDeformationGrid.SetParameterString("ram", "4000") - appFineDeformationGrid.ExecuteAndWriteOutput() - - list_of_Grids.append(os.path.join(burst_dir, fine_grid)) - - ####### Step 3 : SARCoRegistration + SARDeramp + SARRobustInterferogram ####### - ####### SARCoRegistration Application (changeo step) ####### - print("\n SARCoRegistration Application \n", file=stdout_save) - logger.info("SARCoRegistration Application") - slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated.tif" - appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") - appCoRegistration.SetParameterString("insarmaster", os.path.join(burst_dir, burstDerampM)) - appCoRegistration.SetParameterString("insarslave", os.path.join(burst_dir, burstDerampS)) - appCoRegistration.SetParameterString("ingrid", os.path.join(burst_dir, fine_grid)) - appCoRegistration.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appCoRegistration.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appCoRegistration.SetParameterFloat("doppler0", dop0Slave[burstId-firstBurst]) - appCoRegistration.SetParameterInt("sizetiles", 50) - appCoRegistration.SetParameterInt("margin", 7) - appCoRegistration.SetParameterInt("nbramps", 5121) #256*2*10+1 - appCoRegistration.SetParameterString("ram", "4000") - appCoRegistration.SetParameterString("out", os.path.join(burst_dir, slave_Image_CoRe)) - appCoRegistration.ExecuteAndWriteOutput() - - ######## SARDeramp Application (reramp mode) ####### - print("\n Deramping Application \n", file=stdout_save) - logger.info("Deramping Application") - # Slave (CoRegistrated) - burstCoReRerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated_reramp" + ".tif" - appDerampSlave = otb.Registry.CreateApplication("SARDeramp") - appDerampSlave.SetParameterString("in", os.path.join(burst_dir, slave_Image_CoRe)) - appDerampSlave.SetParameterString("out", os.path.join(burst_dir, burstCoReRerampS)) - appDerampSlave.SetParameterString("reramp", "true") - appDerampSlave.SetParameterString("shift", "true") - appDerampSlave.SetParameterString("ingrid", os.path.join(burst_dir, fine_grid)) - appDerampSlave.SetParameterString("inslave", os.path.join(burst_dir, burstDerampS)) - appDerampSlave.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appDerampSlave.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appDerampSlave.SetParameterString("ram", "4000") - appDerampSlave.ExecuteAndWriteOutput() - - ######## SARRobustInterferogram Application (interf step) ####### - print("\n SARRobustInterferogram Application \n", file=stdout_save) - logger.info("SARRobustInterferogram Application") - interferogram_path = os.path.join(burst_dir, "interferogram" + "_burst" + str(burstId) + ".tif") - - if esd_NbIter > 0 : - esd_dir = os.path.join(burst_dir, "esd") - - if not os.path.exists(esd_dir): - os.makedirs(esd_dir) - - interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") - - appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") - appInterferogram.SetParameterString("insarmaster", os.path.join(burst_dir, burstM)) - appInterferogram.SetParameterString("incoregistratedslave", os.path.join(burst_dir, burstCoReRerampS)) - appInterferogram.SetParameterString("insarslave", os.path.join(burst_dir, burstS)) - appInterferogram.SetParameterString("ingrid", os.path.join(burst_dir, fine_grid)) - appInterferogram.SetParameterString("incartmeanmaster", os.path.join(burst_dir, master_cartesian_mean)) - appInterferogram.SetParameterInt("mlran", 1) - appInterferogram.SetParameterInt("mlazi", 1) - appInterferogram.SetParameterInt("marginran", 1) - appInterferogram.SetParameterInt("marginazi", 1) - appInterferogram.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appInterferogram.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appInterferogram.SetParameterFloat("gain", gain_interfero) - appInterferogram.SetParameterString("ram", "4000") - appInterferogram.SetParameterString("out", interferogram_path) - appInterferogram.ExecuteAndWriteOutput() - - list_of_Interferograms.append(interferogram_path) - counter = counter + 1 - - - ######## Step 4 : ESD Loop (SARESD + SARGridOffset + Step 3) ####### - print ("Beginning of ESD Loop", file=stdout_save) - logger.info("Beginning of ESD Loop") + # Create param + param = {} + param['ml_azimut'] = ml_azimut + param['ml_range'] = ml_range + param['validBurstMaster'] = validBurstMaster + param['validBurstSlave'] = validBurstSlave + param['ml_geoGrid_azimut'] = ml_geoGrid_azimut + param['ml_geoGrid_range'] = ml_geoGrid_range + param['geoGrid_gridstep_range'] = geoGrid_gridstep_range + param['geoGrid_gridstep_azimut'] = geoGrid_gridstep_azimut + param['geoGrid_threshold'] = geoGrid_threshold + param['geoGrid_gap'] = geoGrid_gap + param['doppler0'] = dop0Slave + param['gain_interfero'] = gain_interfero + param['advantage'] = advantage + param['esd_NbIter'] = esd_NbIter + param['esd_AutoMode'] = esd_AutoMode + + list_of_Grids, list_of_Interferograms = DInSar.gridToInterferogram(dem, master_Image, master_Image_base, slave_Image, slave_Image_base, output_dir, output_dir, param, 'S1_IW', output_dir) - if esd_NbIter > 0 : - - azimut_shift_esd = [] - azimut_shift_esd_global = [0.] * (lastBurst-firstBurst) - - for iter_esd in range(1, esd_NbIter+1) : - - print ("Iteration number = {number}".format(number=iter_esd), file=stdout_save) - logger.info("Iteration number = {number}".format(number=iter_esd)) - - # Clear all azimut shifts - azimut_shift_esd[:] = [] + print(list_of_Grids) + print(list_of_Interferograms) - for id_loop in range(0, len(validBurstMaster)): - #burstId = id_loop + firstBurst - burstId = validBurstMaster[id_loop] - - print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) - logger.info("BurstId = {ind}".format(ind=burstId)) - - # Paths - burst_dir = os.path.join(output_dir, "burst" + str(burstId)) - esd_dir = os.path.join(burst_dir, "esd") - esd_path = os.path.join(esd_dir, "esdOut" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") - - # Check number of burst index (not SARESD for lastBurst) - if burstId < lastBurst : - ######## SARESD Application ####### - print("\n ESD Application \n", file=stdout_save) - logger.info("Beginning of ESD Loop") - appESD = otb.Registry.CreateApplication("SARESD") - appESD.SetParameterString("ininterfup", list_of_Interferograms[burstId-firstBurst]) - appESD.SetParameterString("ininterflow", list_of_Interferograms[burstId-firstBurst+1]) - appESD.SetParameterString("insar", master_Image) - appESD.SetParameterInt("burstindex", burstId) - appESD.SetParameterFloat("threshold", 0.3) - appESD.SetParameterInt("mlazi", 1) - appESD.SetParameterString("ram", "4000") - appESD.SetParameterString("out", esd_path) - appESD.Execute() - - azimut_shift_esd.append(appESD.GetParameterFloat('azishift')) - - # Clear our list of interferograms - list_of_Interferograms[:] = [] - - for id_loop in range(0, len(validBurstMaster)): - #burstId = id_loop + firstBurst - burstId = validBurstMaster[id_loop] - burstId_slave = validBurstSlave[id_loop] - - # Paths - burst_dir = os.path.join(output_dir, "burst" + str(burstId)) - esd_dir = os.path.join(burst_dir, "esd") - gridOffset_path = os.path.join(esd_dir, "gridOffOut" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") - - burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" - burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" - burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" - burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" - master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" - - # Adjust azimut shift according to the burstId - aziShift = 0. - if burstId == (lastBurst - 1) : - # Only accumulation between iterations - azimut_shift_esd_global[burstId-firstBurst] += azimut_shift_esd[burstId-firstBurst] - aziShift = azimut_shift_esd_global[burstId-firstBurst] - elif burstId == lastBurst : - # Same as the lastBurst -1 - aziShift = azimut_shift_esd_global[burstId - 1 - firstBurst] - else : - # Accumulation of means between the current burstId and the next index - azimut_shift_esd_global[burstId-firstBurst] += ((azimut_shift_esd[burstId-firstBurst] + - azimut_shift_esd[burstId-firstBurst+1])/2) - aziShift = azimut_shift_esd_global[burstId-firstBurst] - - # Offset of deformation grid with the azimut shift - ######## SARGridOffset Application ####### - print("\n GridOffset Application \n", file=stdout_save) - logger.info("GridOffset Application") - appGridOff = otb.Registry.CreateApplication("SARGridOffset") - appGridOff.SetParameterString("ingrid", list_of_Grids[burstId-firstBurst]) - appGridOff.SetParameterFloat("offsetran", 0.) - appGridOff.SetParameterFloat("offsetazi", aziShift) - appGridOff.SetParameterString("ram", "4000") - appGridOff.SetParameterString("out", gridOffset_path) - appGridOff.Execute() - - ######## Step 3 : SARCoRegistration + SARDeramp + SARRobustInterferogram ####### - ######## SARCoRegistration Application (changeo step) ####### - print("\n SARCoRegistration Application \n", file=stdout_save) - logger.info("SARCoRegistration Application") - slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated.tif" - appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") - appCoRegistration.SetParameterString("insarmaster", os.path.join(burst_dir, burstDerampM)) - appCoRegistration.SetParameterString("insarslave", os.path.join(burst_dir, burstDerampS)) - appCoRegistration.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) - appCoRegistration.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appCoRegistration.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appCoRegistration.SetParameterFloat("doppler0", dop0Slave[burstId-firstBurst]) - appCoRegistration.SetParameterInt("sizetiles", 50) - appCoRegistration.SetParameterInt("margin", 7) - appCoRegistration.SetParameterInt("nbramps", 5121) #256*2*10+1 - appCoRegistration.SetParameterString("ram", "4000") - appCoRegistration.SetParameterString("out", os.path.join(esd_dir, slave_Image_CoRe)) - appCoRegistration.Execute() - - ######## SARDeramp Application (reramp mode) ####### - print("\n Deramping Application \n", file=stdout_save) - logger.info("Deramping Application") - # Slave (CoRegistrated) - burstCoReRerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated_reramp" + ".tif" - appDerampSlave = otb.Registry.CreateApplication("SARDeramp") - appDerampSlave.SetParameterInputImage("in", - appCoRegistration.GetParameterOutputImage("out")) - appDerampSlave.SetParameterString("out", os.path.join(esd_dir, burstCoReRerampS)) - appDerampSlave.SetParameterString("reramp", "true") - appDerampSlave.SetParameterString("shift", "true") - appDerampSlave.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) - appDerampSlave.SetParameterString("inslave", os.path.join(burst_dir, burstDerampS)) - appDerampSlave.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appDerampSlave.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appDerampSlave.SetParameterString("ram", "4000") - appDerampSlave.Execute() - - ######## SARRobustInterferogram Application (interf step) ####### - print("\n SARRobustInterferogram Application \n", file=stdout_save) - logger.info("SARRobustInterferogram Application") - interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(iter_esd) + ".tif") - - appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") - appInterferogram.SetParameterString("insarmaster", os.path.join(burst_dir, burstM)) - appInterferogram.SetParameterInputImage("incoregistratedslave", - appDerampSlave.GetParameterOutputImage("out")) - appInterferogram.SetParameterString("insarslave", os.path.join(burst_dir, burstS)) - appInterferogram.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) - appInterferogram.SetParameterString("incartmeanmaster", - os.path.join(burst_dir, master_cartesian_mean)) - appInterferogram.SetParameterInt("mlran", 1) - appInterferogram.SetParameterInt("mlazi", 1) - appInterferogram.SetParameterInt("marginran", 1) - appInterferogram.SetParameterInt("marginazi", 1) - appInterferogram.SetParameterInt("gridsteprange", geoGrid_gridstep_range) - appInterferogram.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) - appInterferogram.SetParameterFloat("gain", gain_interfero) - appInterferogram.SetParameterString("ram", "4000") - appInterferogram.SetParameterString("out", interferogram_path) - appInterferogram.ExecuteAndWriteOutput() - - # Store the new interferogram paths - list_of_Interferograms.append(interferogram_path) - - # Check azimuth shift to end esd loop if automatic mode enabled - if esd_AutoMode : - absMax_aziShift = abs(max(azimut_shift_esd_global, key=abs)) - if absMax_aziShift < 0.01 : - break - - ######## Step 5 : SARConcatenateBursts ####### - ######### SARConcatenateBursts (for interferograms) ######### - print("\n SARConcatenateBursts Application (for interferograms) \n", file=stdout_save) - logger.info("SARConcatenateBursts Application") - interferogram_swath = "interferogram_swath.tif" - appCon = otb.Registry.CreateApplication("SARConcatenateBursts") - appCon.SetParameterStringList("il", list_of_Interferograms) - appCon.SetParameterString("insar", master_Image) - appCon.SetParameterInt("burstindex", firstBurst) - appCon.SetParameterString("out", os.path.join(output_dir, interferogram_swath)) - appCon.SetParameterString("ram", "4000") - appCon.ExecuteAndWriteOutput() - - print("\n End of DiapOTB processing (S1 IW mode) \n", file=stdout_save) - logger.info("############ End of DiapOTB processing (S1 IW mode) ##############") + func_utils.printOnStd("\n End of DiapOTB processing (S1 IW mode) \n") + func_utils.log(logging.INFO, "############ End of DiapOTB processing (S1 IW mode) ##############") diff --git a/python_src/diapOTB_S1IW_OLD.py b/python_src/diapOTB_S1IW_OLD.py new file mode 100644 index 0000000..d5a8364 --- /dev/null +++ b/python_src/diapOTB_S1IW_OLD.py @@ -0,0 +1,905 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +from __future__ import print_function + +__author__ = "POC-THALES" +__version__ = "0.1" +__status__ = "Developpement" +__date__ = "11/12/2018" +__last_modified__ = "11/12/2018" + +# Imports +import logging +import json +from jsonschema import validate +import os +import sys +import argparse +import re + +import otbApplication as otb + +# Streamer to our log file +class StreamToLogger(object): + """ + Fake file-like stream object that redirects writes to a logger instance. + """ + def __init__(self, logger, log_level=logging.INFO): + self.logger = logger + self.log_level = log_level + + def write(self, buf): + for line in buf.rstrip().splitlines(): + self.logger.log(self.log_level, line.rstrip()) + + def flush(self): + for handler in self.logger.handlers: + handler.flush() + + +def validate_json(json, schema): + try: + validate(json, schema) + except Exception as valid_err: + print("Invalid JSON: {}".format(valid_err)) + return False + else: + # Realise votre travail + print("Valid JSON") + return True + +# string to bool +def str2bool(v): + return v.lower() in ("yes", "true", "t", "1") + +# Bursts selection (according to anx time values) +def selectBurst(dictMaster, dictSlave, firstBurst, lastBurst, nbBurstSlave, validBurstMaster, validBurstSlave): + + key1Burst = "support_data.geom.bursts.burst[" + + # Initialize the output lists (empty lists) + validBurstMaster.clear() + validBurstSlave.clear() + + # Loop on Master bursts + for id_B in range(firstBurst, lastBurst+1): + keyBurstMaster = key1Burst + str(id_B) + "].azimuth_anx_time" + + # Get the anx time for the current burst (into Master Image) + anxMaster = float(dictMaster[keyBurstMaster]) + + # Loop on slave bursts to find the closest anx time + minDiff = 200 + id_B_save = id_B + for id_B_slave in range(0, nbBurstSlave): + keyBurstSlave = key1Burst + str(id_B_slave) + "].azimuth_anx_time" + + # Get anx time for slave burst + anxSlave = float(dictSlave[keyBurstSlave]) + + # Comparaison between master and slave + diff = abs(anxMaster - anxSlave) + + if minDiff > diff : + minDiff = diff + id_B_save = id_B_slave + + # Check if difference between the anx time is valid (must be inferior to 1) + if minDiff < 1. : + # Fill lists with master Burst_id and the selected slave burst_id + validBurstMaster.append(id_B) + validBurstSlave.append(id_B_save) + + + +# Main +if __name__ == "__main__": + + # Check arguments + parser = argparse.ArgumentParser() + parser.add_argument("configfile", help="input conguration file for the application DiapOTB") + args = parser.parse_args() + print(args.configfile) + + # Logger initialization + logger = logging.getLogger(__name__) + logger.setLevel(logging.INFO) + + LogFormatter = logging.Formatter('%(filename)s :: %(levelname)s :: %(message)s') + + # Create console handler with a high log level (warning level) + stream_handler = logging.StreamHandler() + stream_handler.setLevel(logging.WARNING) + + # Add Handlers + logger.addHandler(stream_handler) + + # Read and Load the configuration file + try: + with open(args.configfile, 'r') as f: + dataConfig = json.load(f) + + except Exception as err: + logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=args.configfile)) + quit() + + # Load schema (into DiapOTB install) + diapOTB_install = os.getenv('DIAPOTB_INSTALL') + if diapOTB_install is not None and os.path.exists(diapOTB_install): + schemas_path = os.path.join(diapOTB_install, "json_schemas") + if os.path.exists(schemas_path): + schema_S1IW = os.path.join(schemas_path, "schema_S1IW.json") + + try: + with open(schema_S1IW, "r") as sch: + dataSchema = json.load(sch) + except Exception as err: + logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=schema_S1SM)) + quit() + + # Check Json file + jsonIsValid = validate_json(dataConfig, dataSchema) + + if not jsonIsValid : + logger.critical("Error, provided config file does not match requirements") + quit() + + # Get dictionaries + dict_Global = dataConfig['Global'] + dict_PreProcessing = dataConfig['Pre_Processing'] + dict_Ground = dataConfig['Ground'] + dict_DInSAR = dataConfig['DIn_SAR'] + + + # Get elements from dictionaries + # Global + master_Image = dict_Global['in']['Master_Image_Path'] + slave_Image = dict_Global['in']['Slave_Image_Path'] + dem = dict_Global['in']['DEM_Path'] + output_dir = dict_Global['out']['output_dir'] + + # Pre_Processing + ml_range = dict_PreProcessing['parameter']['ML_range'] + ml_azimut = dict_PreProcessing['parameter']['ML_azimut'] + ml_gain = dict_PreProcessing['parameter']['ML_gain'] + dop_file = dict_PreProcessing['out']['doppler_file'] + + # Ground + + # DIn_SAR + geoGrid_gridstep_range = dict_DInSAR['parameter']['GridStep_range'] + geoGrid_gridstep_azimut = dict_DInSAR['parameter']['GridStep_azimut'] + geoGrid_threshold = dict_DInSAR['parameter']['Grid_Threshold'] + geoGrid_gap = dict_DInSAR['parameter']['Grid_Gap'] + ml_geoGrid_range = ml_range + ml_geoGrid_azimut = ml_azimut + gain_interfero = dict_DInSAR['parameter']['Interferogram_gain'] + # esd loop + esd_AutoMode = False # automatic mode to apply a threshold inside the esd loop + esd_NbIter = 0 + + if 'ESD_iter' in dict_DInSAR['parameter']: + esd_NbIter = dict_DInSAR['parameter']['ESD_iter'] + if not isinstance(esd_NbIter, int) : + esd_AutoMode = True + esd_NbIter = 10 # 10 iterations maximum for automatic mode + else : + esd_AutoMode = True + esd_NbIter = 10 # 10 iterations maximum for automatic mode + + + # Check Threshold + if (geoGrid_threshold < 0) or (geoGrid_threshold > 1) : + logger.critical("Error, Wrong Threshold for fine deformation grid") + geoGrid_threshold = 0.3 + + # Check if images/dem exist + if not os.path.exists(master_Image) : + logger.critical("{img} does not exist. Check its path.".format(img=master_Image)) + quit() + if not os.path.exists(slave_Image) : + logger.critical("{img} does not exist. Check its path.".format(img=slave_Image)) + quit() + if not os.path.exists(dem) : + logger.critical("{img} does not exist. Check its path.".format(img=dem)) + quit() + if not os.path.exists(output_dir): + print("The output directory does not exist and will be created") + os.makedirs(output_dir) + else : + print("The output directory exists. Some files can be overwritten") + + + # File handler for the logger + # Create file handler which logs even info messages (used as stdout redirection) + file_handler = logging.FileHandler(os.path.join(output_dir, 'info.log'), 'a') + file_handler.setLevel(logging.INFO) + file_handler.setFormatter(LogFormatter) + + # Add Handlers + logger.addHandler(file_handler) + + # Redirect stdout and stderr to logger + s1 = StreamToLogger(logger, logging.INFO) + stdout_saveWrite = sys.stdout.write # Save stdout.write to print some info into the console + stdout_saveFlush = sys.stdout.flush # Save stdout.flush to print some info into the console + sys.stdout.write = s1.write # Replace stdout.write by our StreamToLogger + sys.stdout.flush = s1.flush # Replace stdout.flush by our StreamToLogger + stdout_save = s1 # Different object + stdout_save.write = stdout_saveWrite # Restore stdout.write into stdout_save + stdout_save.flush = stdout_saveFlush # Restore stdout.write into stdout_save + + # Recap of input parameter into info.log + logger.info("########### Input Parameters for the current execution ############## ") + logger.info(" Pre_Processing : ") + logger.info("ml_range : {param}".format(param=ml_range)) + logger.info("ml_azimut : {param}".format(param=ml_azimut)) + logger.info("ml_gain : {param}".format(param=ml_gain)) + logger.info("dop_file : {param}".format(param=dop_file)) + + # DIn_SAR + logger.info(" DIn_SAR : ") + logger.info("geoGrid_gridstep_range : {param}".format(param=geoGrid_gridstep_range)) + logger.info("geoGrid_gridstep_azimut : {param}".format(param=geoGrid_gridstep_azimut)) + logger.info("geoGrid_threshold : {param}".format(param=geoGrid_threshold)) + logger.info("geoGrid_gap : {param}".format(param=geoGrid_gap)) + logger.info("ml_geoGrid_range : {param}".format(param=ml_geoGrid_range)) + logger.info("ml_geoGrid_azimut : {param}".format(param=ml_geoGrid_azimut)) + logger.info("gain_interfero : {param}".format(param=gain_interfero)) + logger.info("esd_AutoMode : {param}".format(param=esd_AutoMode)) + logger.info("esd_NbIter : {param}".format(param=esd_NbIter)) + + logger.info("########### Input Images for the current execution ############## ") + logger.info("master_Image : {param}".format(param=master_Image)) + logger.info("slave_Image : {param}".format(param=slave_Image)) + logger.info("dem : {param}".format(param=dem)) + + # check Burst index + master_Image_base = os.path.basename(master_Image) + slave_Image_base = os.path.basename(slave_Image) + + # Retrieve some information about our input images + ReadImageInfo = otb.Registry.CreateApplication("ReadImageInfo") + ReadImageInfo.SetParameterString("in", master_Image) + ReadImageInfo.SetParameterString("keywordlist", "true") + ReadImageInfo.Execute() + + keywordMaster = ReadImageInfo.GetParameterString("keyword") + keywordlistMaster = ReadImageInfo.GetParameterString("keyword").split("\n") + keywordlistMaster = filter(None, keywordlistMaster) + dictKWLMaster = { i.split(':')[0] : re.sub(r"[\n\t\s]*", "", i.split(':')[1]) for i in keywordlistMaster } + + ReadImageInfo.SetParameterString("in", slave_Image) + ReadImageInfo.SetParameterString("keywordlist", "true") + ReadImageInfo.Execute() + + keywordSlave = ReadImageInfo.GetParameterString("keyword") + keywordlistSlave = ReadImageInfo.GetParameterString("keyword").split("\n") + keywordlistSlave = filter(None, keywordlistSlave) + dictKWLSlave = { i.split(':')[0] : re.sub(r"[\n\t\s]*", "", i.split(':')[1]) for i in keywordlistSlave } + + # Check header version + if int(dictKWLMaster['header.version']) < 3 or int(dictKWLSlave['header.version']) < 3 : + logger.critical("Upgrade your geom file") + quit() + + # Get information about DEM (spacing, size ..) + ReadDEMInfo = otb.Registry.CreateApplication("ReadImageInfo") + ReadDEMInfo.SetParameterString("in", dem) + ReadDEMInfo.SetParameterString("keywordlist", "true") + ReadDEMInfo.Execute() + + spacingXDEM = ReadDEMInfo.GetParameterFloat("spacingx") + estimatedGroundSpacingXDEM = ReadDEMInfo.GetParameterFloat("estimatedgroundspacingx") + spacingYDEM = ReadDEMInfo.GetParameterFloat("spacingy") + estimatedGroundSpacingYDEM = ReadDEMInfo.GetParameterFloat("estimatedgroundspacingy") + + # Choose advantage for correlation or projection according to DEM resolution + advantage = "projection" # By default projection + + if estimatedGroundSpacingXDEM > 40. or estimatedGroundSpacingYDEM > 40. : + advantage = "correlation" # Correlation if resolution > 40 m + logger.warning("Resolution of the input DEM is inferior to 40 meters : A correlation will be used to correct all deformation grids") + + + # Selection of bursts + keyBurst = "support_data.geom.bursts.burst[" + str(0) + "].azimuth_anx_time" + + # Check the index of bursts + minNbBurst = min([int(dictKWLMaster['support_data.geom.bursts.number']), int(dictKWLSlave['support_data.geom.bursts.number'])]) + + firstBurst = 0 + lastBurst = minNbBurst + burstIndexOk = True + + try: + if 'parameter' in dict_Global: + if 'burst_index' in dict_Global['parameter']: + burstList = dict_Global['parameter']['burst_index'].split('-'); + burstList = [int(i) for i in burstList] + + if len(burstList) == 2 : + firstBurst = min(burstList) + lastBurst = max(burstList) + except Exception as err: + logger.critical("Wrong burst index") + quit() + + if minNbBurst < firstBurst or minNbBurst < lastBurst or lastBurst < 0 or firstBurst < 0 : + logger.critical("Wrong burst index") + quit() + + validBurstMaster = [] + validBurstSlave = [] + nbBurstSlave = int(dictKWLSlave['support_data.geom.bursts.number']) + selectBurst(dictKWLMaster, dictKWLSlave, firstBurst, lastBurst, nbBurstSlave, validBurstMaster, validBurstSlave) + + if len(validBurstMaster) == 0 or len(validBurstSlave) == 0 : + logger.critical("Wrong burst index (slave index does not match with master index)") + quit() + + # Update firstBurst and lastBurst with selected Burst for master image + firstBurst = validBurstMaster[0] + lastBurst = validBurstMaster[len(validBurstMaster)-1] + + + # Create directory for each burst + for burstId in range(validBurstMaster[0], validBurstMaster[len(validBurstMaster)-1]+1): + if not os.path.exists(os.path.join(output_dir, "burst" + str(burstId))): + os.makedirs(os.path.join(output_dir, "burst" + str(burstId))) + + + print("\n Beginning of DiapOTB processing (S1 IW mode) \n", file=stdout_save) + logger.info("############ Beginning of DiapOTB processing (S1 IW mode) ##############") + + ####################### Pre Processing Chain ########################## + # Initialisation of doppler file + dopFile = open(os.path.join(output_dir, dop_file), "w") + dopFile.close() + + # Master + print("\n Master Pre-Processing \n", file=stdout_save) + logger.info("Master Pre-Processing") + + dop0Master = [] + + for id_loop in range(0, len(validBurstMaster)): + #burstId = id_loop + firstBurst + burstId = validBurstMaster[id_loop] + + print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) + logger.info("BurstId = {ind}".format(ind=burstId)) + + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + + ######## SARBurstExtraction Application ####### + print("\n Burst Extraction Application \n", file=stdout_save) + logger.info("Burst Extraction Application") + burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" + appBurstExtractionMaster = otb.Registry.CreateApplication("SARBurstExtraction") + appBurstExtractionMaster.SetParameterString("in", master_Image) + appBurstExtractionMaster.SetParameterString("out", os.path.join(burst_dir, burstM)) + appBurstExtractionMaster.SetParameterInt("burstindex", burstId) + appBurstExtractionMaster.SetParameterString("allpixels", "true") + appBurstExtractionMaster.SetParameterString("ram", "4000") + appBurstExtractionMaster.ExecuteAndWriteOutput() + + ######## SARDeramp Application ####### + print("\n Deramping Application \n", file=stdout_save) + logger.info("Deramping Application") + burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + appDerampMaster = otb.Registry.CreateApplication("SARDeramp") + appDerampMaster.SetParameterString("in", os.path.join(burst_dir, burstM)) + appDerampMaster.SetParameterString("out", os.path.join(burst_dir, burstDerampM)) + appDerampMaster.SetParameterString("ram", "4000") + appDerampMaster.ExecuteAndWriteOutput() + + ######## SARDoppler Application ####### + print("\n Doppler Application \n", file=stdout_save) + logger.info("Doppler Application") + dopFile = open(os.path.join(output_dir, dop_file), "a") + dopFile.write("Doppler for master image : " + os.path.basename(master_Image_base)+ " for burst_index = " + str(burstId) + "\n") + dopFile.close() + appDoppler0Master = otb.Registry.CreateApplication("SARDoppler0") + appDoppler0Master.SetParameterString("insar", os.path.join(burst_dir, burstDerampM)) + appDoppler0Master.SetParameterString("outfile", os.path.join(output_dir, dop_file)) + appDoppler0Master.SetParameterString("ram", "4000") + appDoppler0Master.ExecuteAndWriteOutput() + + dop0Master.append(appDoppler0Master.GetParameterFloat('doppler0')) + + ####### SARMultiLook Application ####### + print("\n MultiLook Application \n", file=stdout_save) + logger.info("MultiLook Application") + master_Image_ML = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + appMultiLookMaster = otb.Registry.CreateApplication("SARMultiLook") + appMultiLookMaster.SetParameterString("incomplex", os.path.join(burst_dir, burstDerampM)) + appMultiLookMaster.SetParameterString("out", os.path.join(burst_dir, master_Image_ML)) + appMultiLookMaster.SetParameterInt("mlran",ml_range) + appMultiLookMaster.SetParameterInt("mlazi",ml_azimut) + appMultiLookMaster.SetParameterFloat("mlgain", ml_gain) + appMultiLookMaster.SetParameterString("ram", "4000") + appMultiLookMaster.ExecuteAndWriteOutput() + + # Slave + print("\n Slave Pre-Processing \n", file=stdout_save) + logger.info("Slave Pre-Processing") + + dop0Slave = [] + + for id_loop in range(0, len(validBurstMaster)): + #burstId = id_loop + firstBurst + burstId = validBurstMaster[id_loop] + burstId_slave = validBurstSlave[id_loop] + + print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) + logger.info("BurstId = {ind}".format(ind=burstId)) + + + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + + ######## SARBurstExtraction Application ####### + print("\n Burst Extraction Application \n", file=stdout_save) + logger.info("Burst Extraction Application") + burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" + appBurstExtractionSlave = otb.Registry.CreateApplication("SARBurstExtraction") + appBurstExtractionSlave.SetParameterString("in", slave_Image) + appBurstExtractionSlave.SetParameterString("out", os.path.join(burst_dir, burstS)) + appBurstExtractionSlave.SetParameterInt("burstindex", burstId_slave) + appBurstExtractionSlave.SetParameterString("allpixels", "true") + appBurstExtractionSlave.SetParameterString("ram", "4000") + appBurstExtractionSlave.ExecuteAndWriteOutput() + + ######## SARDeramp Application ####### + print("\n Deramping Application \n", file=stdout_save) + logger.info("Deramping Application") + burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" + appDerampSlave = otb.Registry.CreateApplication("SARDeramp") + appDerampSlave.SetParameterString("in", os.path.join(burst_dir, burstS)) + appDerampSlave.SetParameterString("out", os.path.join(burst_dir, burstDerampS)) + appDerampSlave.SetParameterString("ram", "4000") + appDerampSlave.ExecuteAndWriteOutput() + + ######## SARDoppler Application ####### + print("\n Doppler Application \n", file=stdout_save) + logger.info("Doppler Application") + dopFile = open(os.path.join(output_dir, dop_file), "a") + dopFile.write("Doppler for slave image : " + os.path.basename(slave_Image_base) + " for burst_index = " + str(burstId_slave) + "\n") + dopFile.close() + appDoppler0Slave = otb.Registry.CreateApplication("SARDoppler0") + appDoppler0Slave.SetParameterString("insar", os.path.join(burst_dir, burstDerampS)) + appDoppler0Slave.SetParameterString("outfile", os.path.join(output_dir, dop_file)) + appDoppler0Slave.SetParameterString("ram", "4000") + appDoppler0Slave.ExecuteAndWriteOutput() + + dop0Slave.append(appDoppler0Slave.GetParameterFloat('doppler0')) + + ####### SARMultiLook Application ####### + print("\n MultiLook Application \n", file=stdout_save) + logger.info("MultiLook Application") + slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + appMultiLookSlave = otb.Registry.CreateApplication("SARMultiLook") + appMultiLookSlave.SetParameterString("incomplex", os.path.join(burst_dir, burstDerampS)) + appMultiLookSlave.SetParameterString("out", os.path.join(burst_dir, slave_Image_ML)) + appMultiLookSlave.SetParameterInt("mlran",ml_range) + appMultiLookSlave.SetParameterInt("mlazi",ml_azimut) + appMultiLookSlave.SetParameterFloat("mlgain", ml_gain) + appMultiLookSlave.SetParameterString("ram", "4000") + appMultiLookSlave.ExecuteAndWriteOutput() + + + ######################### Ground Chain ############################# + # Master + print("\n Master Ground chain \n", file=stdout_save) + logger.info("Master Ground Application") + + gainMaster = [] + directionDEMMaster = [] + + for id_loop in range(0, len(validBurstMaster)): + #burstId = id_loop + firstBurst + burstId = validBurstMaster[id_loop] + + print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) + logger.info("BurstId = {ind}".format(ind=burstId)) + + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + + burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + + + ######## SARDEMProjection Application ####### + print("\n SARDEMProjection Application \n", file=stdout_save) + logger.info("SARDEMProjection Application") + demProj_Master = "demProj" + "_burst" + str(burstId) +"_Master.tif" + appDEMProjectionMaster = otb.Registry.CreateApplication("SARDEMProjection") + appDEMProjectionMaster.SetParameterString("insar", os.path.join(burst_dir, burstDerampM)) + appDEMProjectionMaster.SetParameterString("indem", dem) + appDEMProjectionMaster.SetParameterString("withxyz", "true") + appDEMProjectionMaster.SetParameterInt("nodata", -32768) + appDEMProjectionMaster.SetParameterString("out", os.path.join(burst_dir, demProj_Master)) + appDEMProjectionMaster.SetParameterString("ram", "4000") + appDEMProjectionMaster.ExecuteAndWriteOutput() + + gainMaster.append(appDEMProjectionMaster.GetParameterFloat('gain')) + directionDEMMaster.append([appDEMProjectionMaster.GetParameterInt('directiontoscandemc'), appDEMProjectionMaster.GetParameterInt('directiontoscandeml')]) + + ######### SARCartesianMeanEstimation Application ####### + print("\n SARCartesianMeanEstimation Application \n", file=stdout_save) + logger.info("SARCartesianMeanEstimation Application") + master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" + appCartMeanMaster = otb.Registry.CreateApplication("SARCartesianMeanEstimation") + appCartMeanMaster.SetParameterString("insar", os.path.join(burst_dir, burstDerampM)) + appCartMeanMaster.SetParameterString("indem", dem) + appCartMeanMaster.SetParameterString("indemproj", os.path.join(burst_dir, demProj_Master)) + appCartMeanMaster.SetParameterInt("indirectiondemc", appDEMProjectionMaster.GetParameterInt("directiontoscandemc")) + appCartMeanMaster.SetParameterInt("indirectiondeml", appDEMProjectionMaster.GetParameterInt("directiontoscandeml")) + appCartMeanMaster.SetParameterInt("mlran", 1) + appCartMeanMaster.SetParameterInt("mlazi", 1) + appCartMeanMaster.SetParameterString("ram", "4000") + appCartMeanMaster.SetParameterString("out", os.path.join(burst_dir, master_cartesian_mean)) + appCartMeanMaster.ExecuteAndWriteOutput() + + + print (directionDEMMaster) + + # Slave + print("\n Slave Ground chain \n", file=stdout_save) + logger.info("Slave Ground Application") + + gainSlave = [] + directionDEMSlave = [] + + for id_loop in range(0, len(validBurstMaster)): + #burstId = id_loop + firstBurst + burstId = validBurstMaster[id_loop] + burstId_slave = validBurstSlave[id_loop] + + print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) + logger.info("BurstId = {ind}".format(ind=burstId)) + + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + + burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" + + print("\n SARDEMProjection Application \n", file=stdout_save) + logger.info("SARDEMProjection Application") + demProj_Slave = "demProj" + "_burst" + str(burstId_slave) + "_Slave.tif" + appDEMProjectionSlave = otb.Registry.CreateApplication("SARDEMProjection") + appDEMProjectionSlave.SetParameterString("insar", os.path.join(burst_dir, burstDerampS)) + appDEMProjectionSlave.SetParameterString("indem", dem) + appDEMProjectionSlave.SetParameterString("withxyz", "true"); + appDEMProjectionSlave.SetParameterInt("nodata", -32768) + appDEMProjectionSlave.SetParameterString("out", os.path.join(burst_dir, demProj_Slave)) + appDEMProjectionSlave.SetParameterString("ram", "4000") + appDEMProjectionSlave.ExecuteAndWriteOutput() + + + ######################## DIn_SAR Chain ############################# + print("\n DIn_SAR chain \n", file=stdout_save) + logger.info("DIn_SAR chain") + + counter = 0 + list_of_Interferograms = [] + list_of_Grids = [] + + for id_loop in range(0, len(validBurstMaster)): + #burstId = id_loop + firstBurst + burstId = validBurstMaster[id_loop] + burstId_slave = validBurstSlave[id_loop] + + print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) + logger.info("BurstId = {ind}".format(ind=burstId)) + + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + + burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" + burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" + + burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" + + master_Image_ML = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + + + demProj_Master = "demProj" + "_burst" + str(burstId) +"_Master.tif" + demProj_Slave = "demProj" + "_burst" + str(burstId_slave) + "_Slave.tif" + master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" + + + ######## Step 1 : SARFineDeformationGrid ####### + ######## SARFineDeformationGrid Application (geo_grid step) ####### + print("\n SARFineDeformationGrid Application \n", file=stdout_save) + logger.info("SARFineDeformationGrid Application") + fine_grid = "fineDeformationGrid"+ "_burst" + str(burstId) + ".tif" + appFineDeformationGrid = otb.Registry.CreateApplication("SARFineDeformationGrid") + appFineDeformationGrid.SetParameterString("indem", dem) + appFineDeformationGrid.SetParameterString("insarmaster", os.path.join(burst_dir, burstDerampM)) + appFineDeformationGrid.SetParameterString("insarslave", os.path.join(burst_dir, burstDerampS)) + appFineDeformationGrid.SetParameterString("inmlmaster", os.path.join(burst_dir, master_Image_ML)) + appFineDeformationGrid.SetParameterString("inmlslave", os.path.join(burst_dir, slave_Image_ML)) + appFineDeformationGrid.SetParameterString("indemprojmaster", os.path.join(burst_dir, demProj_Master)) + appFineDeformationGrid.SetParameterString("indemprojslave", os.path.join(burst_dir, demProj_Slave)) + appFineDeformationGrid.SetParameterInt("mlran", ml_geoGrid_range) + appFineDeformationGrid.SetParameterInt("mlazi", ml_geoGrid_azimut) + appFineDeformationGrid.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appFineDeformationGrid.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold) + appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap) + appFineDeformationGrid.SetParameterString("advantage", advantage) + appFineDeformationGrid.SetParameterString("out", os.path.join(burst_dir, fine_grid)) + appFineDeformationGrid.SetParameterString("ram", "4000") + appFineDeformationGrid.ExecuteAndWriteOutput() + + list_of_Grids.append(os.path.join(burst_dir, fine_grid)) + + ####### Step 3 : SARCoRegistration + SARDeramp + SARRobustInterferogram ####### + ####### SARCoRegistration Application (changeo step) ####### + print("\n SARCoRegistration Application \n", file=stdout_save) + logger.info("SARCoRegistration Application") + slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated.tif" + appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") + appCoRegistration.SetParameterString("insarmaster", os.path.join(burst_dir, burstDerampM)) + appCoRegistration.SetParameterString("insarslave", os.path.join(burst_dir, burstDerampS)) + appCoRegistration.SetParameterString("ingrid", os.path.join(burst_dir, fine_grid)) + appCoRegistration.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appCoRegistration.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appCoRegistration.SetParameterFloat("doppler0", dop0Slave[burstId-firstBurst]) + appCoRegistration.SetParameterInt("sizetiles", 50) + appCoRegistration.SetParameterInt("margin", 7) + appCoRegistration.SetParameterInt("nbramps", 5121) #256*2*10+1 + appCoRegistration.SetParameterString("ram", "4000") + appCoRegistration.SetParameterString("out", os.path.join(burst_dir, slave_Image_CoRe)) + appCoRegistration.ExecuteAndWriteOutput() + + ######## SARDeramp Application (reramp mode) ####### + print("\n Deramping Application \n", file=stdout_save) + logger.info("Deramping Application") + # Slave (CoRegistrated) + burstCoReRerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated_reramp" + ".tif" + appDerampSlave = otb.Registry.CreateApplication("SARDeramp") + appDerampSlave.SetParameterString("in", os.path.join(burst_dir, slave_Image_CoRe)) + appDerampSlave.SetParameterString("out", os.path.join(burst_dir, burstCoReRerampS)) + appDerampSlave.SetParameterString("reramp", "true") + appDerampSlave.SetParameterString("shift", "true") + appDerampSlave.SetParameterString("ingrid", os.path.join(burst_dir, fine_grid)) + appDerampSlave.SetParameterString("inslave", os.path.join(burst_dir, burstDerampS)) + appDerampSlave.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appDerampSlave.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appDerampSlave.SetParameterString("ram", "4000") + appDerampSlave.ExecuteAndWriteOutput() + + ######## SARRobustInterferogram Application (interf step) ####### + print("\n SARRobustInterferogram Application \n", file=stdout_save) + logger.info("SARRobustInterferogram Application") + interferogram_path = os.path.join(burst_dir, "interferogram" + "_burst" + str(burstId) + ".tif") + + if esd_NbIter > 0 : + esd_dir = os.path.join(burst_dir, "esd") + + if not os.path.exists(esd_dir): + os.makedirs(esd_dir) + + interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") + + appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") + appInterferogram.SetParameterString("insarmaster", os.path.join(burst_dir, burstM)) + appInterferogram.SetParameterString("incoregistratedslave", os.path.join(burst_dir, burstCoReRerampS)) + appInterferogram.SetParameterString("insarslave", os.path.join(burst_dir, burstS)) + appInterferogram.SetParameterString("ingrid", os.path.join(burst_dir, fine_grid)) + appInterferogram.SetParameterString("incartmeanmaster", os.path.join(burst_dir, master_cartesian_mean)) + appInterferogram.SetParameterInt("mlran", 1) + appInterferogram.SetParameterInt("mlazi", 1) + appInterferogram.SetParameterInt("marginran", 1) + appInterferogram.SetParameterInt("marginazi", 1) + appInterferogram.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appInterferogram.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appInterferogram.SetParameterFloat("gain", gain_interfero) + appInterferogram.SetParameterString("ram", "4000") + appInterferogram.SetParameterString("out", interferogram_path) + appInterferogram.ExecuteAndWriteOutput() + + list_of_Interferograms.append(interferogram_path) + counter = counter + 1 + + + ######## Step 4 : ESD Loop (SARESD + SARGridOffset + Step 3) ####### + print ("Beginning of ESD Loop", file=stdout_save) + logger.info("Beginning of ESD Loop") + + if esd_NbIter > 0 : + + azimut_shift_esd = [] + azimut_shift_esd_global = [0.] * (lastBurst-firstBurst) + + for iter_esd in range(1, esd_NbIter+1) : + + print ("Iteration number = {number}".format(number=iter_esd), file=stdout_save) + logger.info("Iteration number = {number}".format(number=iter_esd)) + + # Clear all azimut shifts + azimut_shift_esd[:] = [] + + for id_loop in range(0, len(validBurstMaster)): + #burstId = id_loop + firstBurst + burstId = validBurstMaster[id_loop] + + print("\n BurstId = " + str(burstId) + "\n", file=stdout_save) + logger.info("BurstId = {ind}".format(ind=burstId)) + + # Paths + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + esd_dir = os.path.join(burst_dir, "esd") + esd_path = os.path.join(esd_dir, "esdOut" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") + + # Check number of burst index (not SARESD for lastBurst) + if burstId < lastBurst : + ######## SARESD Application ####### + print("\n ESD Application \n", file=stdout_save) + logger.info("Beginning of ESD Loop") + appESD = otb.Registry.CreateApplication("SARESD") + appESD.SetParameterString("ininterfup", list_of_Interferograms[burstId-firstBurst]) + appESD.SetParameterString("ininterflow", list_of_Interferograms[burstId-firstBurst+1]) + appESD.SetParameterString("insar", master_Image) + appESD.SetParameterInt("burstindex", burstId) + appESD.SetParameterFloat("threshold", 0.3) + appESD.SetParameterInt("mlazi", 1) + appESD.SetParameterString("ram", "4000") + appESD.SetParameterString("out", esd_path) + appESD.Execute() + + azimut_shift_esd.append(appESD.GetParameterFloat('azishift')) + + # Clear our list of interferograms + list_of_Interferograms[:] = [] + + for id_loop in range(0, len(validBurstMaster)): + #burstId = id_loop + firstBurst + burstId = validBurstMaster[id_loop] + burstId_slave = validBurstSlave[id_loop] + + # Paths + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + esd_dir = os.path.join(burst_dir, "esd") + gridOffset_path = os.path.join(esd_dir, "gridOffOut" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") + + burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" + burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" + burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" + master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" + + # Adjust azimut shift according to the burstId + aziShift = 0. + if burstId == (lastBurst - 1) : + # Only accumulation between iterations + azimut_shift_esd_global[burstId-firstBurst] += azimut_shift_esd[burstId-firstBurst] + aziShift = azimut_shift_esd_global[burstId-firstBurst] + elif burstId == lastBurst : + # Same as the lastBurst -1 + aziShift = azimut_shift_esd_global[burstId - 1 - firstBurst] + else : + # Accumulation of means between the current burstId and the next index + azimut_shift_esd_global[burstId-firstBurst] += ((azimut_shift_esd[burstId-firstBurst] + + azimut_shift_esd[burstId-firstBurst+1])/2) + aziShift = azimut_shift_esd_global[burstId-firstBurst] + + # Offset of deformation grid with the azimut shift + ######## SARGridOffset Application ####### + print("\n GridOffset Application \n", file=stdout_save) + logger.info("GridOffset Application") + appGridOff = otb.Registry.CreateApplication("SARGridOffset") + appGridOff.SetParameterString("ingrid", list_of_Grids[burstId-firstBurst]) + appGridOff.SetParameterFloat("offsetran", 0.) + appGridOff.SetParameterFloat("offsetazi", aziShift) + appGridOff.SetParameterString("ram", "4000") + appGridOff.SetParameterString("out", gridOffset_path) + appGridOff.Execute() + + ######## Step 3 : SARCoRegistration + SARDeramp + SARRobustInterferogram ####### + ######## SARCoRegistration Application (changeo step) ####### + print("\n SARCoRegistration Application \n", file=stdout_save) + logger.info("SARCoRegistration Application") + slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated.tif" + appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") + appCoRegistration.SetParameterString("insarmaster", os.path.join(burst_dir, burstDerampM)) + appCoRegistration.SetParameterString("insarslave", os.path.join(burst_dir, burstDerampS)) + appCoRegistration.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) + appCoRegistration.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appCoRegistration.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appCoRegistration.SetParameterFloat("doppler0", dop0Slave[burstId-firstBurst]) + appCoRegistration.SetParameterInt("sizetiles", 50) + appCoRegistration.SetParameterInt("margin", 7) + appCoRegistration.SetParameterInt("nbramps", 5121) #256*2*10+1 + appCoRegistration.SetParameterString("ram", "4000") + appCoRegistration.SetParameterString("out", os.path.join(esd_dir, slave_Image_CoRe)) + appCoRegistration.Execute() + + ######## SARDeramp Application (reramp mode) ####### + print("\n Deramping Application \n", file=stdout_save) + logger.info("Deramping Application") + # Slave (CoRegistrated) + burstCoReRerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated_reramp" + ".tif" + appDerampSlave = otb.Registry.CreateApplication("SARDeramp") + appDerampSlave.SetParameterInputImage("in", + appCoRegistration.GetParameterOutputImage("out")) + appDerampSlave.SetParameterString("out", os.path.join(esd_dir, burstCoReRerampS)) + appDerampSlave.SetParameterString("reramp", "true") + appDerampSlave.SetParameterString("shift", "true") + appDerampSlave.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) + appDerampSlave.SetParameterString("inslave", os.path.join(burst_dir, burstDerampS)) + appDerampSlave.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appDerampSlave.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appDerampSlave.SetParameterString("ram", "4000") + appDerampSlave.Execute() + + ######## SARRobustInterferogram Application (interf step) ####### + print("\n SARRobustInterferogram Application \n", file=stdout_save) + logger.info("SARRobustInterferogram Application") + interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(iter_esd) + ".tif") + + appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") + appInterferogram.SetParameterString("insarmaster", os.path.join(burst_dir, burstM)) + appInterferogram.SetParameterInputImage("incoregistratedslave", + appDerampSlave.GetParameterOutputImage("out")) + appInterferogram.SetParameterString("insarslave", os.path.join(burst_dir, burstS)) + appInterferogram.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) + appInterferogram.SetParameterString("incartmeanmaster", + os.path.join(burst_dir, master_cartesian_mean)) + appInterferogram.SetParameterInt("mlran", 1) + appInterferogram.SetParameterInt("mlazi", 1) + appInterferogram.SetParameterInt("marginran", 1) + appInterferogram.SetParameterInt("marginazi", 1) + appInterferogram.SetParameterInt("gridsteprange", geoGrid_gridstep_range) + appInterferogram.SetParameterInt("gridstepazimut", geoGrid_gridstep_azimut) + appInterferogram.SetParameterFloat("gain", gain_interfero) + appInterferogram.SetParameterString("ram", "4000") + appInterferogram.SetParameterString("out", interferogram_path) + appInterferogram.ExecuteAndWriteOutput() + + # Store the new interferogram paths + list_of_Interferograms.append(interferogram_path) + + # Check azimuth shift to end esd loop if automatic mode enabled + if esd_AutoMode : + absMax_aziShift = abs(max(azimut_shift_esd_global, key=abs)) + if absMax_aziShift < 0.01 : + break + + ######## Step 5 : SARConcatenateBursts ####### + ######### SARConcatenateBursts (for interferograms) ######### + print("\n SARConcatenateBursts Application (for interferograms) \n", file=stdout_save) + logger.info("SARConcatenateBursts Application") + interferogram_swath = "interferogram_swath.tif" + appCon = otb.Registry.CreateApplication("SARConcatenateBursts") + appCon.SetParameterStringList("il", list_of_Interferograms) + appCon.SetParameterString("insar", master_Image) + appCon.SetParameterInt("burstindex", firstBurst) + appCon.SetParameterString("out", os.path.join(output_dir, interferogram_swath)) + appCon.SetParameterString("ram", "4000") + appCon.ExecuteAndWriteOutput() + + + print("\n End of DiapOTB processing (S1 IW mode) \n", file=stdout_save) + logger.info("############ End of DiapOTB processing (S1 IW mode) ##############") diff --git a/python_src/processings/DInSar.py b/python_src/processings/DInSar.py new file mode 100644 index 0000000..0bb219f --- /dev/null +++ b/python_src/processings/DInSar.py @@ -0,0 +1,553 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +""" + The ``DINSAR`` chain + ====================== + + Use the module to launch DInSar processigns. + + main function + -------------- + + gridToInterferogram + +""" + +import os +import sys +from functools import partial + +import utils.DiapOTB_applications as diapOTBApp +import otbApplication as otb + + +# Coorection for interferograms thanks to esd shifts +def esd_correct_interferogram(validBurstMaster, validBurstSlave, master_Image_base, slave_Image_base, gridsteprange, gridstepazimut, gain, output_dir, master_dir, doppler0, iter_esd, list_of_Grids, azimut_shift_esd, azimut_shift_esd_global) : + """ + Correct interferograms with esd shifts. + + This function is used inside the ESD loop (called by gridToInterferogram_S1IW). + + :param validBurstMaster: Master burst to process + :param validBurstSlave: Slave burst to process + :param master_Image_base: Master image name + :param slave_Image_base: Slave image name + :param gridsteprange: Grid step for range dimension (into SLC/SAR geometry). + :param gridstepazimut: Grid step for azimut dimension (into SLC/SAR geometry). + :param gain: Gain to apply for amplitude estimation. + :param output_dir: Output directory + :param master_dir: Output directory (to store result on master image only) + :param doppler0: Doppler0 values for slave bursts + :param iter_esd: Current iteration number + :param list_of_Grids: List of deformation grid (one per burst) + :param azimut_shift_esd: List of azimut shift found by SARESD Application (one per burst) + :param azimut_shift_esd_global: List of azimut shift added for each iteration (one per burst) + + :type validBurstMaster: list + :type validBurstSlave: list + :type master_Image_base: string + :type slave_Image_base: string + :type gridsteprange: int + :type gridstepazimut: int + :type gain: float + :type output_dir: string + :type master_dir: string + :type doppler0: list + :type iter_esd: int + :type list_of_Grids: list + :type azimut_shift_esd: list + :type azimut_shift_esd_global: list + + """ + firstBurst = validBurstMaster[0] + lastBurst = validBurstMaster[len(validBurstMaster)-1] + + list_of_Interferograms = [] + + for id_loop in range(0, len(validBurstMaster)): + burstId = validBurstMaster[id_loop] + burstId_slave = validBurstSlave[id_loop] + + # Paths + Master_burst_dir = os.path.join(master_dir, "burst" + str(burstId)) + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + esd_dir = os.path.join(burst_dir, "esd") + gridOffset_path = os.path.join(esd_dir, "gridOffOut" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") + + burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" + burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" + burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" + master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" + + # Adjust azimut shift according to the burstId + aziShift = 0. + if burstId == (lastBurst - 1) : + # Only accumulation between iterations + azimut_shift_esd_global[burstId-firstBurst] += azimut_shift_esd[burstId-firstBurst] + aziShift = azimut_shift_esd_global[burstId-firstBurst] + elif burstId == lastBurst : + # Same as the lastBurst -1 + aziShift = azimut_shift_esd_global[burstId - 1 - firstBurst] + else : + # Accumulation of means between the current burstId and the next index + azimut_shift_esd_global[burstId-firstBurst] += ((azimut_shift_esd[burstId-firstBurst] + + azimut_shift_esd[burstId-firstBurst+1])/2) + aziShift = azimut_shift_esd_global[burstId-firstBurst] + + # Offset of deformation grid with the azimut shift + ## SARGridOffset Application ## + print("\n GridOffset Application \n") + appGridOff = otb.Registry.CreateApplication("SARGridOffset") + appGridOff.SetParameterString("ingrid", list_of_Grids[burstId-firstBurst]) + appGridOff.SetParameterFloat("offsetran", 0.) + appGridOff.SetParameterFloat("offsetazi", aziShift) + appGridOff.SetParameterString("ram", "4000") + appGridOff.SetParameterString("out", gridOffset_path) + appGridOff.Execute() + + + + ######## Step 3 : SARCoRegistration + SARDeramp + SARRobustInterferogram ####### + ## SARCoRegistration Application (changeo step) ## + print("\n SARCoRegistration Application \n") + slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated.tif" + coRe_path = os.path.join(esd_dir, slave_Image_CoRe) + appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") + appCoRegistration.SetParameterString("insarmaster", os.path.join(burst_dir, burstDerampM)) + appCoRegistration.SetParameterString("insarslave", os.path.join(burst_dir, burstDerampS)) + appCoRegistration.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) + appCoRegistration.SetParameterInt("gridstepazimut", gridstepazimut) + appCoRegistration.SetParameterInt("gridsteprange", gridsteprange) + appCoRegistration.SetParameterFloat("doppler0", doppler0[burstId-validBurstMaster[0]]) + appCoRegistration.SetParameterInt("sizetiles", 50) + appCoRegistration.SetParameterInt("margin", 7) + appCoRegistration.SetParameterInt("nbramps", 5121) #256*2*10+1 + appCoRegistration.SetParameterString("ram", "4000") + appCoRegistration.SetParameterString("out", coRe_path) + appCoRegistration.Execute() + + + ######## SARDeramp Application (reramp mode) ####### + print("\n Deramping Application \n") + # Slave (CoRegistrated) + burstCoReRerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated_reramp" + ".tif" + deramp_Path = os.path.join(esd_dir, burstCoReRerampS) + + appDerampSlave = otb.Registry.CreateApplication("SARDeramp") + appDerampSlave.SetParameterInputImage("in", appCoRegistration.GetParameterOutputImage("out")) + appDerampSlave.SetParameterString("out", deramp_Path) + appDerampSlave.SetParameterString("reramp", "true") + appDerampSlave.SetParameterString("shift", "true") + appDerampSlave.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) + appDerampSlave.SetParameterString("inslave", os.path.join(burst_dir, burstDerampS)) + appDerampSlave.SetParameterInt("gridsteprange", gridsteprange) + appDerampSlave.SetParameterInt("gridstepazimut", gridstepazimut) + appDerampSlave.SetParameterString("ram", "4000") + appDerampSlave.Execute() + + + ## SARRobustInterferogram Application (interf step) ## + print("\n SARRobustInterferogram Application \n") + interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(iter_esd) + ".tif") + + appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") + appInterferogram.SetParameterString("insarmaster", os.path.join(Master_burst_dir, burstM)) + appInterferogram.SetParameterInputImage("incoregistratedslave", appDerampSlave.GetParameterOutputImage("out")) + appInterferogram.SetParameterString("insarslave", os.path.join(burst_dir, burstS)) + appInterferogram.SetParameterInputImage("ingrid", appGridOff.GetParameterOutputImage("out")) + appInterferogram.SetParameterString("incartmeanmaster", os.path.join(Master_burst_dir, master_cartesian_mean)) + appInterferogram.SetParameterInt("mlran", 1) + appInterferogram.SetParameterInt("mlazi", 1) + appInterferogram.SetParameterInt("marginran", 1) + appInterferogram.SetParameterInt("marginazi", 1) + appInterferogram.SetParameterInt("gridsteprange", gridsteprange) + appInterferogram.SetParameterInt("gridstepazimut", gridstepazimut) + appInterferogram.SetParameterFloat("gain", gain) + appInterferogram.SetParameterString("ram", "4000") + appInterferogram.SetParameterString("out", interferogram_path) + appInterferogram.ExecuteAndWriteOutput() + + list_of_Interferograms.append(interferogram_path) + + return list_of_Interferograms + + +def gridToInterferogram(dem, master_Image, master_Image_base, slave_Image, slave_Image_base, master_dir, slave_dir, param, mode, output_dir): + """ + Main function for DInSAR processing + + This function applies several processing to provide the result interferogram + + :param dem: Input DEM + :param master_Image: Master image name (with path included) + :param master_Image_base: Master image name + :param slave_Image: Slave image name (with path included) + :param slave_Image_base: Slave image name + :param master_dir: Output directory (to store result on master image only) + :param slave_dir: Output directory (to store result on slave image only) + :param param: Dictionary to retrieve all parameters used into this processing chain + :param mode: Selected mode according to input product : S1IW or other (for Cosmo and S1 StripMap products) + :param output_dir: Output directory + + :type dem: string + :type master_Image: string + :type master_Image_base: string + :type slave_Image: string + :type slave_Image_base: string + :type master_dir: string + :type slave_dir: string + :type param: dict + :type mode: string + :type output_dir: string + + :return: list of grids and interferograms + :rtype: list, list + """ + if mode == 'S1_IW' : + return gridToInterferogram_S1IW(dem, master_Image, master_Image_base, slave_Image, slave_Image_base, master_dir, slave_dir, param, output_dir) + else : + return gridToInterferogram_Others(dem, master_Image, master_Image_base, slave_Image, slave_Image_base, master_dir, slave_dir, param, output_dir) + + + +def gridToInterferogram_Others(dem, master_Image, master_Image_base, slave_Image, slave_Image_base, master_dir, slave_dir, param, output_dir): + """ + Main function for DInSAR processing and other than S1IW mode (S1SM and Cosmo) + + This function applies several processing to provide the result interferogram + + :param dem: Input DEM + :param master_Image: Master image name (with path included) + :param master_Image_base: Master image name + :param slave_Image: Slave image name (with path included) + :param slave_Image_base: Slave image name + :param master_dir: Output directory (to store result on master image only) + :param slave_dir: Output directory (to store result on slave image only) + :param param: Dictionary to retrieve all parameters used into this processing chain + :param output_dir: Output directory + + :type dem: string + :type master_Image: string + :type master_Image_base: string + :type slave_Image: string + :type slave_Image_base: string + :type master_dir: string + :type slave_dir: string + :type param: dict + :type output_dir: string + + :return: list of grids and interferograms + :rtype: list, list + """ + # Get elements into param dictionnaries + ml_azimut = param['ml_azimut'] + ml_range = param['ml_range'] + ml_geoGrid_azimut = param['ml_geoGrid_azimut'] + ml_geoGrid_range = param['ml_geoGrid_range'] + geoGrid_gridstep_range = param['geoGrid_gridstep_range'] + geoGrid_gridstep_azimut = param['geoGrid_gridstep_azimut'] + geoGrid_threshold = param['geoGrid_threshold'] + geoGrid_gap = param['geoGrid_gap'] + doppler0 = param['doppler0'] + gain_interfero = param['gain_interfero'] + advantage = param['advantage'] + + print("doppler0 = " + str(doppler0)) + + nbramps = 257 + + # define applications functions with the previous parameters + fineGridDeformationApp = partial(diapOTBApp.fineGridDeformation, indem=dem, mlran=ml_geoGrid_range, + mlazi=ml_geoGrid_azimut, gridsteprange=geoGrid_gridstep_range, + gridstepazimut=geoGrid_gridstep_azimut, + threshold=geoGrid_threshold, gap=geoGrid_gap, advantage=advantage) + + coRegistrationApp = partial(diapOTBApp.coRegistration, gridsteprange=geoGrid_gridstep_range, gridstepazimut=geoGrid_gridstep_azimut, nbramps=nbramps) + + derampApp = partial(diapOTBApp.deramp, gridsteprange=geoGrid_gridstep_range, gridstepazimut=geoGrid_gridstep_azimut, reramp="true", shift="true") + + interferogramApp = partial(diapOTBApp.interferogram, gridsteprange=geoGrid_gridstep_range, gridstepazimut=geoGrid_gridstep_azimut, mlran=ml_range, mlazi=ml_azimut, gain=gain_interfero) + + + # Empty list + list_of_Grids = [] + list_of_Interferograms = [] + + + master_Image_ML = os.path.splitext(master_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + + + demProj_Master = "demProj_Master.tif" + demProj_Slave = "demProj_Slave.tif" + master_cartesian_mean = "CartMeanMaster.tif" + + ## SARFineDeformationGrid Application ## + print("\n SARFineDeformationGrid Application \n") + fine_grid = "fineDeformationGrid.tif" + fine_gird_path = os.path.join(output_dir, fine_grid) + fineGridDeformationApp(insarmaster=master_Image, + insarslave=slave_Image, + inmlmaster=os.path.join(master_dir, master_Image_ML), + inmlslave=os.path.join(output_dir, slave_Image_ML), + indemprojmaster=os.path.join(master_dir, demProj_Master), + indemprojslave=os.path.join(output_dir, demProj_Slave), + outPath=fine_gird_path) + + list_of_Grids.append(fine_gird_path) + + + ## SARCoRegistration Application (changeo step) ## + print("\n SARCoRegistration Application \n") + slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_coregistrated.tif" + coRe_path = os.path.join(output_dir, slave_Image_CoRe) + coRegistrationApp(insarmaster=master_Image, + insarslave=slave_Image, + ingrid=fine_gird_path, doppler0=doppler0, + outPath=coRe_path) + + + + ## SARRobustInterferogram Application (interf step) ## + print("\n SARRobustInterferogram Application \n") + interferogram_path = os.path.join(output_dir, "interferogram.tif") + + interferogramApp(insarmaster=master_Image, + incoregistratedslave=os.path.join(output_dir, slave_Image_CoRe), + insarslave=slave_Image, + ingrid=fine_gird_path, + incartmeanmaster=os.path.join(master_dir, master_cartesian_mean), + outPath=interferogram_path) + + list_of_Interferograms.append(interferogram_path) + + # Output lists + return list_of_Grids, list_of_Interferograms + + +def gridToInterferogram_S1IW(dem, master_Image, master_Image_base, slave_Image, slave_Image_base, master_dir, slave_dir, param, output_dir): + """ + Main function for DInSAR processing and S1IW mode + + This function applies several processing to provide the result interferogram (with concatenation of all interferograms) + + :param dem: Input DEM + :param master_Image: Master image name (with path included) + :param master_Image_base: Master image name + :param slave_Image: Slave image name (with path included) + :param slave_Image_base: Slave image name + :param master_dir: Output directory (to store result on master image only) + :param slave_dir: Output directory (to store result on slave image only) + :param param: Dictionary to retrieve all parameters used into this processing chain + :param output_dir: Output directory + + :type dem: string + :type master_Image: string + :type master_Image_base: string + :type slave_Image: string + :type slave_Image_base: string + :type master_dir: string + :type slave_dir: string + :type param: dict + :type output_dir: string + + :return: list of grids and interferograms + :rtype: list, list + """ + # Get elements into param dictionnaries + ml_azimut = param['ml_azimut'] + ml_range = param['ml_range'] + validBurstMaster = param['validBurstMaster'] + validBurstSlave = param['validBurstSlave'] + ml_geoGrid_azimut = param['ml_geoGrid_azimut'] + ml_geoGrid_range = param['ml_geoGrid_range'] + geoGrid_gridstep_range = param['geoGrid_gridstep_range'] + geoGrid_gridstep_azimut = param['geoGrid_gridstep_azimut'] + geoGrid_threshold = param['geoGrid_threshold'] + geoGrid_gap = param['geoGrid_gap'] + doppler0 = param['doppler0'] + gain_interfero = param['gain_interfero'] + advantage = param['advantage'] + esd_NbIter = param['esd_NbIter'] + esd_AutoMode = param['esd_AutoMode'] + + nbramps = 256*2*10+1 + + # define applications functions with the previous parameters + fineGridDeformationApp = partial(diapOTBApp.fineGridDeformation, indem=dem, mlran=ml_geoGrid_range, mlazi=ml_geoGrid_azimut, gridsteprange=geoGrid_gridstep_range, gridstepazimut=geoGrid_gridstep_azimut, + threshold=geoGrid_threshold, gap=geoGrid_gap, advantage=advantage) + + coRegistrationApp = partial(diapOTBApp.coRegistration, gridsteprange=geoGrid_gridstep_range, gridstepazimut=geoGrid_gridstep_azimut, nbramps=nbramps) + + derampApp = partial(diapOTBApp.deramp, gridsteprange=geoGrid_gridstep_range, gridstepazimut=geoGrid_gridstep_azimut, reramp="true", shift="true") + + interferogramApp = partial(diapOTBApp.interferogram, gridsteprange=geoGrid_gridstep_range, gridstepazimut=geoGrid_gridstep_azimut, mlran=1, mlazi=1, gain=gain_interfero) + + esdApp = partial(diapOTBApp.esd, insar=master_Image) + + esdCorrectInterferogram = partial(esd_correct_interferogram, validBurstMaster, validBurstSlave, master_Image_base, slave_Image_base, geoGrid_gridstep_range, geoGrid_gridstep_azimut, gain_interfero, output_dir, master_dir, doppler0) + + # Empty list + list_of_Grids = [] + list_of_Interferograms = [] + + + for id_loop in range(0, len(validBurstMaster)): + burstId = validBurstMaster[id_loop] + burstId_slave = validBurstSlave[id_loop] + + print("\n BurstId = " + str(burstId) + "\n") + + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + Master_burst_dir = os.path.join(master_dir, "burst" + str(burstId)) + + burstM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + ".tif" + burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + ".tif" + + burstDerampM = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_deramp" + ".tif" + + master_Image_ML = os.path.splitext(master_Image_base)[0] + "_burst" + str(burstId) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId_slave) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + + + demProj_Master = "demProj" + "_burst" + str(burstId) +"_Master.tif" + demProj_Slave = "demProj" + "_burst" + str(burstId_slave) + "_Slave.tif" + master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" + + counter = 0 + + ## SARFineDeformationGrid Application ## + print("\n SARFineDeformationGrid Application \n") + fine_grid = "fineDeformationGrid"+ "_burst" + str(burstId) + ".tif" + fine_gird_path = os.path.join(burst_dir, fine_grid) + fineGridDeformationApp(insarmaster=os.path.join(Master_burst_dir, burstDerampM), + insarslave=os.path.join(burst_dir, burstDerampS), + inmlmaster=os.path.join(Master_burst_dir, master_Image_ML), + inmlslave=os.path.join(burst_dir, slave_Image_ML), + indemprojmaster=os.path.join(Master_burst_dir, demProj_Master), + indemprojslave=os.path.join(burst_dir, demProj_Slave), + outPath=fine_gird_path) + + list_of_Grids.append(os.path.join(burst_dir, fine_grid)) + + + ## SARCoRegistration Application (changeo step) ## + print("\n SARCoRegistration Application \n") + slave_Image_CoRe = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated.tif" + coRe_path = os.path.join(burst_dir, slave_Image_CoRe) + coRegistrationApp(insarmaster=os.path.join(Master_burst_dir, burstDerampM), + insarslave=os.path.join(burst_dir, burstDerampS), + ingrid=fine_gird_path, doppler0=doppler0[burstId-validBurstMaster[0]], + outPath=coRe_path) + + + ## SARDeramp Application (reramp mode) ## + print("\n Deramping Application \n") + # Slave (CoRegistrated) + burstCoReRerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_coregistrated_reramp" + ".tif" + deramp_Path = os.path.join(burst_dir, burstCoReRerampS) + derampApp(inImg=os.path.join(burst_dir, slave_Image_CoRe), ingrid=fine_gird_path, + inslave=os.path.join(burst_dir, burstDerampS), outPath=deramp_Path) + + ## SARRobustInterferogram Application (interf step) ## + print("\n SARRobustInterferogram Application \n") + interferogram_path = os.path.join(burst_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") + if esd_NbIter > 0 : + esd_dir = os.path.join(burst_dir, "esd") + + if not os.path.exists(esd_dir): + os.makedirs(esd_dir) + + interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") + + interferogramApp(insarmaster=os.path.join(Master_burst_dir, burstM), + incoregistratedslave=deramp_Path, insarslave=os.path.join(burst_dir, burstS), + ingrid=fine_gird_path, + incartmeanmaster=os.path.join(Master_burst_dir, master_cartesian_mean), + outPath=interferogram_path) + + counter = counter + 1 + list_of_Interferograms.append(interferogram_path) + + ### ESD Loop (if required) #### + print ("Beginning of ESD Loop") + + firstBurst = validBurstMaster[0] + lastBurst = validBurstMaster[len(validBurstMaster)-1] + + if esd_NbIter > 0 : + + azimut_shift_esd = [] + azimut_shift_esd_global = [0.] * (lastBurst-firstBurst) + + for iter_esd in range(1, esd_NbIter+1) : + + print ("Iteration number = {number}".format(number=iter_esd)) + + # Clear all azimut shifts + azimut_shift_esd[:] = [] + + for id_loop in range(0, len(validBurstMaster)): + burstId = validBurstMaster[id_loop] + + print("\n BurstId = " + str(burstId) + "\n") + + # Paths + burst_dir = os.path.join(output_dir, "burst" + str(burstId)) + esd_dir = os.path.join(burst_dir, "esd") + esd_path = os.path.join(esd_dir, "esdOut" + "_burst" + str(burstId) + "_iter" + str(0) + ".tif") + + # Check number of burst index (not SARESD for lastBurst) + if burstId < lastBurst : + ## SARESD Application ## + print("\n ESD Application \n") + aziShiftEsd = esdApp(ininterfup=list_of_Interferograms[burstId-firstBurst], + ininterflow=list_of_Interferograms[burstId-firstBurst+1], + burstIndex=burstId, outPath=esd_path) + + + print("aziShift = " + str(aziShiftEsd)) + azimut_shift_esd.append(aziShiftEsd) + + # Correction of Interferogram with the esd shift + # new reference => new list of interferograms + list_of_Interferograms = esdCorrectInterferogram(iter_esd, list_of_Grids, azimut_shift_esd, + azimut_shift_esd_global) + + # Check azimuth shift to end esd loop if automatic mode enabled + if esd_AutoMode : + absMax_aziShift = abs(max(azimut_shift_esd_global, key=abs)) + if absMax_aziShift < 0.01 : + break + + + # Concatenate interferograms + ## SARConcatenateBursts (for interferograms) ## + print("\n SARConcatenateBursts Application (for interferograms) \n") + interferogram_swath = "interferogram_swath.tif" + diapOTBApp.concatenate(list_of_Interferograms, master_Image, firstBurst, os.path.join(output_dir, interferogram_swath)) + + # Output lists + return list_of_Grids, list_of_Interferograms diff --git a/python_src/processings/Ground.py b/python_src/processings/Ground.py new file mode 100644 index 0000000..10ed049 --- /dev/null +++ b/python_src/processings/Ground.py @@ -0,0 +1,186 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +""" + The ``Ground`` chain + ============================ + + Use the module to launch Ground chain. + + main function + -------------- + + demProjectionAndCartesianEstimation + +""" + +import os +import sys +from functools import partial + +import utils.DiapOTB_applications as diapOTBApp + + + +def demProjectionAndCartesianEstimation(sar_Image, sar_Image_base, dem, param, mode, output_dir): + """ + Main function for Ground chain + + This function applies several processing according to the selected mode : + For S1IW mode : For each burst, Burst Extraction + Deramping + Doppler0 estimation + MultiLook processing + For other mode : Doppler0 estimation + MultiLook processing + + :param sar_Image: SAR image name (with path included) + :param sar_Image_base: SAR image name + :param dem: Input DEM + :param param: Dictionary to retrieve all parameters used into this processing chain + :param mode: Selected mode according to input product : S1IW or other (for Cosmo and S1 StripMap products) + :param output_dir: Output directory + + :type sar_Image: string + :type sar_Image_base: string + :type dem: string + :type param: dict + :type mode: string + :type output_dir: string + """ + if mode == 'S1_IW' : + return demProjectionAndCartesianEstimation_S1IW(sar_Image, sar_Image_base, dem, param, output_dir) + else : + return demProjectionAndCartesianEstimation_Others(sar_Image, sar_Image_base, dem, param, output_dir) + + +def demProjectionAndCartesianEstimation_Others(sar_Image, sar_Image_base, dem, param, output_dir): + """ + Main function for Ground chain and other than S1IW mode (S1SM and Cosmo) + + This function applies several processing DEM Projection + Cartesian estimation + + :param sar_Image: SAR image name (with path included) + :param sar_Image_base: SAR image name + :param dem: Input DEM + :param param: Dictionary to retrieve all parameters used into this processing chain + :param output_dir: Output directory + + :type sar_Image: string + :type sar_Image_base: string + :type dem: string + :type param: dict + :type output_dir: string + """ + + for_slave_Image = False + + # Get elements into param dictionnaries + nodata = param['nodata'] + withxyz = param['withxyz'] + for_slave_Image = param['for_slave_Image'] + + + ## SARDEMProjection Application ## + print("\n SARDEMProjection Application \n") + extProj = "_Master.tif" + if for_slave_Image : + extProj = "_Slave.tif" + + demProj = "demProj" +extProj + + directiontoscandemc, directiontoscandeml, gain = diapOTBApp.demProjection(insar=sar_Image, indem=dem, withxyz=withxyz, outPath=os.path.join(output_dir, demProj), nodata=nodata) + + + if not for_slave_Image : + ## SARCartesianMeanEstimation Application ## + print("\n SARCartesianMeanEstimation Application \n") + master_cartesian_mean = "CartMeanMaster.tif" + + diapOTBApp.cartesianMeanEstimation(insar=sar_Image, indem=dem, + indemproj=os.path.join(output_dir, demProj), + indirectiondemc=directiontoscandemc, indirectiondeml=directiontoscandeml, + outPath=os.path.join(output_dir, master_cartesian_mean)) + + +def demProjectionAndCartesianEstimation_S1IW(sar_Image, sar_Image_base, dem, param, output_dir): + """ + Main function for Ground chain and S1IW mode + + This function applies several processing DEM Projection + Cartesian estimation + + :param sar_Image: SAR image name (with path included) + :param sar_Image_base: SAR image name + :param dem: Input DEM + :param param: Dictionary to retrieve all parameters used into this processing chain + :param output_dir: Output directory + + :type sar_Image: string + :type sar_Image_base: string + :type dem: string + :type param: dict + :type output_dir: string + """ + + for_slave_Image = False + + # Get elements into param dictionnaries + nodata = param['nodata'] + withxyz = param['withxyz'] + validBurstMaster = param['validBurstMaster'] + validBurstSlave = [] + if "validBurstSlave" in param : + validBurstSlave = param['validBurstSlave'] + for_slave_Image = True + + # Loop on burst + for id_loop in range(0, len(validBurstMaster)): + + burstId_dir = validBurstMaster[id_loop] + + burstId = burstId_dir + if for_slave_Image : + burstId = validBurstSlave[id_loop] + + print("\n BurstId = " + str(burstId_dir) + "\n") + + burst_dir = os.path.join(output_dir, "burst" + str(burstId_dir)) + + burstDeramp = os.path.splitext(sar_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + + + ## SARDEMProjection Application ## + print("\n SARDEMProjection Application \n") + extProj = "_Master.tif" + if for_slave_Image : + extProj = "_Slave.tif" + + demProj = "demProj" + "_burst" + str(burstId) +extProj + + directiontoscandemc, directiontoscandeml, gain = diapOTBApp.demProjection(insar=os.path.join(burst_dir, burstDeramp), indem=dem, withxyz=withxyz, outPath=os.path.join(burst_dir, demProj), nodata=nodata) + + + if not for_slave_Image : + ## SARCartesianMeanEstimation Application ## + print("\n SARCartesianMeanEstimation Application \n") + master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif" + + diapOTBApp.cartesianMeanEstimation(insar=os.path.join(burst_dir, burstDeramp), indem=dem, + indemproj=os.path.join(burst_dir, demProj), + indirectiondemc=directiontoscandemc, indirectiondeml=directiontoscandeml, + outPath=os.path.join(burst_dir, master_cartesian_mean)) diff --git a/python_src/processings/Metadata_Correction.py b/python_src/processings/Metadata_Correction.py new file mode 100644 index 0000000..05927e8 --- /dev/null +++ b/python_src/processings/Metadata_Correction.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +""" + The ``Metadata Correction`` chain + ================================= + + Use the module to launch Metadata Correction chain (only available for S1 SM or Cosmo. + + main function + -------------- + + metadataCorrection + +""" + +import os +import sys +from functools import partial + +import utils.DiapOTB_applications as diapOTBApp + + + +def fineMetadata(sar_Image, sar_Image_base, dem, param, output_dir): + """ + Main function for Metadata Correction chain + + This function applies several processing (only available for S1 SM and Cosmo products) : + Dem to Amplitude application + Correlation grid + Metadata Correction + + :param sar_Image: SAR image name (with path included) + :param sar_Image_base: SAR image name + :param dem: Input DEM + :param param: Dictionary to retrieve all parameters used into this processing chain + :param output_dir: Output directory + + :type sar_Image: string + :type sar_Image_base: string + :type dem: string + :type param: dict + :type output_dir: string + """ + + # Get elements into param dictionnaries + ml_azimut = param['ml_azimut'] + ml_range = param['ml_range'] + ml_gain = param['ml_gain'] + geoGrid_gridstep_range = param['geoGrid_gridstep_range'] + geoGrid_gridstep_azimut = param['geoGrid_gridstep_azimut'] + nodata = param['nodata'] + fine_metadata_file = param['fine_metadata_file'] + + ## SARDEMToAmplitude Application ## + amplitude_simu_image = os.path.splitext(sar_Image_base)[0] + "_simuSAR" + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + diapOTBApp.demToAmplitude(insar=sar_Image, dem=dem, mlran=ml_range, mlazi=ml_azimut, mlgain=ml_gain, + outPath=os.path.join(output_dir, amplitude_simu_image), nodata=nodata) + + + ## SARCorrelationGrid Application ## + correl_grid = "correl_simu" + "_gridstep" + str(geoGrid_gridstep_azimut) + str(geoGrid_gridstep_range) + ".tif" + sar_Image_ML = os.path.splitext(sar_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + + diapOTBApp.correlationGrid(inmaster=os.path.join(output_dir, sar_Image_ML), + inslave=os.path.join(output_dir, amplitude_simu_image), + mlran=ml_range, mlazi=ml_azimut, gridsteprange=geoGrid_gridstep_range, + gridstepazimut=geoGrid_gridstep_azimut, + outPath=os.path.join(output_dir, correl_grid)) + + ## SARFineMetadata Application ## + diapOTBApp.fineMetadata(insar=sar_Image, ingrid=os.path.join(output_dir, correl_grid), + stepmax=0.1, threshold=0.3, outPath=os.path.join(output_dir, fine_metadata_file)) + + diff --git a/python_src/processings/Pre_Processing.py b/python_src/processings/Pre_Processing.py new file mode 100644 index 0000000..04be2a7 --- /dev/null +++ b/python_src/processings/Pre_Processing.py @@ -0,0 +1,204 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +""" + The ``Pre_Processing`` chain + ============================ + + Use the module to launch Pre_Processing chain. + + main function + -------------- + + extractToMultilook + +""" + +import os +import sys +from functools import partial + +import utils.DiapOTB_applications as diapOTBApp + + + +def extractToMultilook(sar_Image, sar_Image_base, param, mode, output_dir): + """ + Main function for Pre_Processing chain + + This function applies several processing according to the selected mode : + For S1IW mode : For each burst, Burst Extraction + Deramping + Doppler0 estimation + MultiLook processing + For other mode : Doppler0 estimation + MultiLook processing + + :param sar_Image: SAR image name (with path included) + :param sar_Image_base: SAR image name + :param param: Dictionary to retrieve all parameters used into this processing chain + :param mode: Selected mode according to input product : S1IW or other (for Cosmo and S1 StripMap products) + :param output_dir: Output directory + + :type sar_Image: string + :type sar_Image_base: string + :type param: dict + :type mode: string + :type output_dir: string + + :return: Doppler0 values + :rtype: list + """ + if mode == 'S1_IW' : + return extractToMultilook_S1IW(sar_Image, sar_Image_base, param, output_dir) + else : + return extractToMultilook_Others(sar_Image, sar_Image_base, param, output_dir) + + +def extractToMultilook_Others(sar_Image, sar_Image_base, param, output_dir): + """ + Main function for Pre_Processing chain and other than S1IW mode (S1SM and Cosmo) + + This function applies several processing according to the selected mode : + For other mode : Doppler0 estimation + MultiLook processing + + :param sar_Image: SAR image name (with path included) + :param sar_Image_base: SAR image name + :param param: Dictionary to retrieve all parameters used into this processing chain + :param output_dir: Output directory + + :type sar_Image: string + :type sar_Image_base: string + :type param: dict + :type output_dir: string + + :return: Doppler0 value + :rtype: float + """ + + # Get elements into param dictionnaries + ml_azimut = param['ml_azimut'] + ml_range = param['ml_range'] + ml_gain = param['ml_gain'] + dop_file = param['dop_file'] + + # Define some applications functions with the previous parameters + multilookApp = partial(diapOTBApp.multilook, mlran=ml_range, mlazi=ml_azimut, mlgain=ml_gain) + + # Create an empty list + dop0 = 0 + + ## SARDoppler Application ## + print("\n Doppler Application \n") + dop0 = diapOTBApp.doppler0(insar=sar_Image, + outPath=os.path.join(output_dir, dop_file)) + + ## SARMultiLook Application ## + print("\n MultiLook Application \n") + sar_Image_ML = os.path.splitext(sar_Image_base)[0] + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + + multilookApp(inImg=sar_Image, + outPath=os.path.join(output_dir, sar_Image_ML)) + + # return the value of Doppler0 values + return dop0 + + + +def extractToMultilook_S1IW(sar_Image, sar_Image_base, param, output_dir): + """ + Main function for Pre_Processing chain and S1IW mode + + This function applies several processing according to the selected mode : + For S1IW mode : For each burst, Burst Extraction + Deramping + Doppler0 estimation + MultiLook processing + For other mode : Doppler0 estimation + MultiLook processing + + :param sar_Image: SAR image name (with path included) + :param sar_Image_base: SAR image name + :param param: Dictionary to retrieve all parameters used into this processing chain + :param output_dir: Output directory + + :type sar_Image: string + :type sar_Image_base: string + :type param: dict + :type output_dir: string + + :return: Doppler0 values + :rtype: list + """ + + for_slave_Image = False + + # Get elements into param dictionnaries + ml_azimut = param['ml_azimut'] + ml_range = param['ml_range'] + ml_gain = param['ml_gain'] + dop_file = param['dop_file'] + validBurstMaster = param['validBurstMaster'] + validBurstSlave = [] + if "validBurstSlave" in param : + validBurstSlave = param['validBurstSlave'] + for_slave_Image = True + + # Define some applications functions with the previous parameters + multilookApp = partial(diapOTBApp.multilook, mlran=ml_range, mlazi=ml_azimut, mlgain=ml_gain) + + derampApp = partial(diapOTBApp.deramp, reramp="false", shift="false", ingrid="", inslave="", + gridsteprange=0, gridstepazimut=0) + + # Create an empty list + dop0 = [] + + # Loop on bursts + for id_loop in range(0, len(validBurstMaster)): + + burstId_dir = validBurstMaster[id_loop] + + burstId = burstId_dir + if for_slave_Image : + burstId = validBurstSlave[id_loop] + + print("\n BurstId = " + str(burstId_dir) + "\n") + + burst_dir = os.path.join(output_dir, "burst" + str(burstId_dir)) + + ## SARBurstExtraction Application ## + print("\n Burst Extraction Application \n") + burstM = os.path.splitext(sar_Image_base)[0] + "_burst" + str(burstId) + ".tif" + + diapOTBApp.burstExtraction(sar_Image, burstId, os.path.join(burst_dir, burstM)) + + ## SARDeramp Application ## + print("\n Deramping Application \n") + burstDerampM = os.path.splitext(sar_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif" + derampApp(inImg=os.path.join(burst_dir, burstM), outPath=os.path.join(burst_dir, burstDerampM)) + + ## SARDoppler Application ## + print("\n Doppler Application \n") + dop0.append(diapOTBApp.doppler0(insar=os.path.join(burst_dir, burstDerampM), + outPath=os.path.join(output_dir, dop_file))) + + ## SARMultiLook Application ## + print("\n MultiLook Application \n") + sar_Image_ML = os.path.splitext(sar_Image_base)[0] + "_burst" + str(burstId) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif" + + multilookApp(inImg=os.path.join(burst_dir, burstDerampM), + outPath=os.path.join(burst_dir, sar_Image_ML)) + + # return the list of Doppler0 values + return dop0 diff --git a/python_src/processings/__init__.py b/python_src/processings/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python_src/utils/DiapOTB_applications.py b/python_src/utils/DiapOTB_applications.py new file mode 100644 index 0000000..8d57fec --- /dev/null +++ b/python_src/utils/DiapOTB_applications.py @@ -0,0 +1,584 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +""" + DiapOTB applications + ==================== + + Use the module to DiapOTB and OTB applications from Python + Applications are used with ExecuteAndWriteOutput as execute command => I/O on disk and not integrated + for a mode In-Memory. + + Applications list + ----------------- + + burstExtraction + doppler0 + multilook + demProjection + cartesianMeanEstimation + fineGridDeformation + coRegistration + deramp + interferogram + esd + concatenate + + +""" + +import os +import sys + +from utils import func_utils as fUtils +import otbApplication as otb + +##### Define function for OTB applications used into this processing ##### +def burstExtraction(inImg, burstIndex, outPath, allpixels="true", ram="4000"): + """ + Launch SARBurstExtraction application (from OTB). + + This application performs a burst extraction by keeping only lines and samples of a required burst. + Two modes are available for the output image : with all pixels and with only valid pixels. + + :param inImg: Input Sentinel1 IW SLC Image. + :param burstIndex: Index for the required Burst. + :param allpixels: Select the modes for output image. + :param outPath: Output burst. + + :type inImg: string + :type burstIndex: int + :type allpixels: bool + :type outPath: string + + :return: None + """ + appBurstExtraction = otb.Registry.CreateApplication("SARBurstExtraction") + appBurstExtraction.SetParameterString("in", inImg) + appBurstExtraction.SetParameterString("out", outPath) + appBurstExtraction.SetParameterInt("burstindex", burstIndex) + appBurstExtraction.SetParameterString("allpixels", allpixels) + appBurstExtraction.SetParameterString("ram", ram) + appBurstExtraction.ExecuteAndWriteOutput() + +def doppler0(insar, outPath, ram="4000"): + """ + Launch SARDoppler0 application (from DiapOTB). + + This application estimates the value of Doppler 0 for a SAR image. The algorithm is based on the CDE + method described into [Estimating the Doppler centroid of SAR data]. + + :param insar: Input SAR image. + :param outPath: Output file to store Doppler 0 value.. + + :type insar: string + :type outPath: string + + :return: Doppler0 value (azmiut dimension) + :rtype: float + """ + appDoppler0 = otb.Registry.CreateApplication("SARDoppler0") + appDoppler0.SetParameterString("insar", insar) + appDoppler0.SetParameterString("outfile", outPath) + appDoppler0.SetParameterString("ram", "4000") + appDoppler0.ExecuteAndWriteOutput() + + return appDoppler0.GetParameterFloat('doppler0') + +def multilook(inImg, mlran, mlazi, mlgain, outPath, ram="4000"): + """ + Launch SARMultiLook application (from DiapOTB). + + This application creates the multi-look image of a SLC product. + + :param inImg: Input image (complex here). + :param mlran: Averaging on distance). + :param mlazi: Averaging on azimut. + :param mlgain: Gain to apply on ML image. + :param outPath: Output ML Image. + + :type inImg: string + :type mlran: int + :type mlazi: int + :type mlgain: float + :type outPath: string + + :return: None + """ + appMultiLook = otb.Registry.CreateApplication("SARMultiLook") + appMultiLook.SetParameterString("incomplex", inImg) + appMultiLook.SetParameterString("out", outPath) + appMultiLook.SetParameterInt("mlran",mlran) + appMultiLook.SetParameterInt("mlazi",mlazi) + appMultiLook.SetParameterFloat("mlgain", mlgain) + appMultiLook.SetParameterString("ram", ram) + appMultiLook.ExecuteAndWriteOutput() + + +def demToAmplitude(insar, dem, mlran, mlazi, mlgain, outPath, nodata=-32768, ram="4000"): + """ + Launch SARDEMToAmplitude application (from DiapOTB). + + This application executes the Simu_SAR step with three internal applications : + SARDEMProjection, SARAmplitudeEstimation and SARMultiLook. + + :param insar: SAR Image to extract SAR geometry. + :param indem: DEM to extract DEM geometry. + :param mlran: Averaging on distance). + :param mlazi: Averaging on azimut. + :param mlgain: Gain to apply on ML image. + :param nodata: No Data values for the DEM. + :param outPath: Output Amplitude Image. + + :type insar: string + :type indem: string + :type mlran: int + :type mlazi: int + :type mlgain: float + :type nodata: int + :type outPath: string + + :return: None + """ + appDEMToAmplitude = otb.Registry.CreateApplication("SARDEMToAmplitude") + appDEMToAmplitude.SetParameterString("insar", insar) + appDEMToAmplitude.SetParameterString("indem", dem) + appDEMToAmplitude.SetParameterString("out", outPath) + appDEMToAmplitude.SetParameterInt("mlran", mlran) + appDEMToAmplitude.SetParameterInt("mlazi", mlazi) + appDEMToAmplitude.SetParameterFloat("mlgain", mlgain) + appDEMToAmplitude.SetParameterInt("nodata", nodata) + appDEMToAmplitude.SetParameterString("ram", ram) + appDEMToAmplitude.ExecuteAndWriteOutput() + + +def demProjection(insar, indem, withxyz, outPath, nodata=-32768, ram="4000"): + """ + Launch SARDEMProjection application (from DiapOTB). + + This application puts a DEM file into SAR geometry and estimates two additional coordinates. + In all for each point of the DEM input four components are calculated : + C (colunm into SAR image), L (line into SAR image), Z and Y. + + :param insar: SAR Image to extract SAR geometry. + :param indem: DEM to perform projection on. + :param withxyz: Set XYZ Cartesian components into projection. + :param nodata: No Data values for the DEM. + :param outPath: Output Vector Image with at least four components (C,L,Z,Y). + + :type insar: string + :type indem: string + :type withxyz: bool + :type nodata: int + :type outPath: string + + :return: Range/Azimut direction for DEM scan and Gain value + :rtype: int, int, float + """ + appDEMProjection = otb.Registry.CreateApplication("SARDEMProjection") + appDEMProjection.SetParameterString("insar", insar) + appDEMProjection.SetParameterString("indem", indem) + appDEMProjection.SetParameterString("withxyz", withxyz) + appDEMProjection.SetParameterInt("nodata", nodata) + appDEMProjection.SetParameterString("out", outPath) + appDEMProjection.SetParameterString("ram", ram) + appDEMProjection.ExecuteAndWriteOutput() + + return appDEMProjection.GetParameterInt('directiontoscandemc'), appDEMProjection.GetParameterInt('directiontoscandeml'), appDEMProjection.GetParameterFloat('gain') + +def cartesianMeanEstimation(insar,indem, indemproj, indirectiondemc, indirectiondeml, outPath, ram="4000"): + """ + Launch SARCartesianMeanEstimation application (from DiapOTB). + + This application estimates a simulated cartesian mean image thanks to a DEM file. + + :param insar: SAR Image to extract SAR geometry. + :param indem: DEM to extract DEM geometry. + :param indemproj: Input vector image for cartesian mean estimation. + :param indirectiondemc: Range direction for DEM scan + :param indirectiondeml: Azimut direction for DEM scan + :param outPath: Output cartesian (mean) Image for DEM Projection. + + :type insar: string + :type indem: string + :type indemproj: string + :type indirectiondemc: int + :type indirectiondeml: int + :type outPath: string + + :return: None + """ + appCartMean = otb.Registry.CreateApplication("SARCartesianMeanEstimation") + appCartMean.SetParameterString("insar", insar) + appCartMean.SetParameterString("indem", indem) + appCartMean.SetParameterString("indemproj", indemproj) + appCartMean.SetParameterInt("indirectiondemc", indirectiondemc) + appCartMean.SetParameterInt("indirectiondeml", indirectiondeml) + appCartMean.SetParameterInt("mlran", 1) + appCartMean.SetParameterInt("mlazi", 1) + appCartMean.SetParameterString("ram", ram) + appCartMean.SetParameterString("out", outPath) + appCartMean.ExecuteAndWriteOutput() + + +def correlationGrid(inmaster, inslave, mlran, mlazi, gridsteprange, gridstepazimut, outPath, ram="4000"): + """ + Launch SARCorrelationGrid application (from DiapOTB). + + This application computes correlation shifts between two images : + shift in range and shift in azimut. + The inputs of this application are MultiLooked images (real images) + + :param inmaster: Master Image (real image). + :param inslave: Slave Image (real image). + :param mlran: MultiLook factor on distance. + :param mlazi: MultiLook factor on azimut. + :param gridsteprange: Grid step for range dimension (into SLC/SAR geometry). + :param gridstepazimut: Grid step for azimut dimension (into SLC/SAR geometry). + :param outPath: Output Correlation Grid Vector Image (Shift_ran, Shift_azi). + + :type inmaster: string + :type inslave: string + :type mlran: int + :type mlazi: int + :type gridsteprange: int + :type gridstepazimut: int + :type outPath: string + + :return: None + """ + appCorGrid = otb.Registry.CreateApplication("SARCorrelationGrid") + appCorGrid.SetParameterString("inmaster", inmaster) + appCorGrid.SetParameterString("inslave", inslave) + appCorGrid.SetParameterInt("mlran", mlran) + appCorGrid.SetParameterInt("mlazi", mlazi) + appCorGrid.SetParameterInt("gridsteprange", gridsteprange) + appCorGrid.SetParameterInt("gridstepazimut", gridstepazimut) + appCorGrid.SetParameterString("out", outPath) + appCorGrid.SetParameterString("ram", ram) + appCorGrid.ExecuteAndWriteOutput() + + +def fineMetadata(insar, ingrid, stepmax, threshold, outPath, ram="4000"): + """ + Launch SARFineMetadata application (from DiapOTB). + + This application corrects two metadata : + the time of the first line and the slant near range. + + :param insar: SAR Image to extract SAR geometry. + :param ingrid: Input correlation grid. + :param stepmax: Max step for histogram. + :param threshold: Threshold on correlation rate. + :param outPath: Output file to store new metadata values. + + :type insar: string + :type ingrid: string + :type stepmax: float + :type threshold: float + :type outPath: string + + :return: None + """ + appFineMetadata = otb.Registry.CreateApplication("SARFineMetadata") + appFineMetadata.SetParameterString("insar", insar) + appFineMetadata.SetParameterString("ingrid", ingrid) + appFineMetadata.SetParameterFloat("stepmax", stepmax) + appFineMetadata.SetParameterFloat("threshold", threshold) + appFineMetadata.SetParameterString("outfile", outPath) + appFineMetadata.ExecuteAndWriteOutput() + + +def fineGridDeformation(indem, insarmaster, insarslave, inmlmaster, inmlslave, indemprojmaster, indemprojslave, mlran, mlazi, gridsteprange, gridstepazimut, threshold, gap, advantage, outPath, ram="4000"): + """ + Launch SARFineDeformationGrid application (from DiapOTB). + + This application executes the geo_grid step with + three internal applications : SARCorrelationGrid, SARDEMGrid and + SARCorrectionGrid. The aim is to obtain a fine deformation grid between + master and slave SAR images + + :param indem: DEM to extract DEM geometry. + :param insarmaster: SAR Master Image to extract SAR geometry. + :param insarslave: SAR Slave Image to extract SAR geometry. + :param inmlmaster: Master Image Multilooked (real image). + :param inmlslave: Slave Image Multilooked (real image). + :param indemprojmaster: Input vector of DEM projected into SAR Master geometry. + :param indemprojslave: Input vector of DEM projected into SAR Slave geometry. + :param mlran: MultiLook factor on distance. + :param mlazi: MultiLook factor on azimut. + :param gridsteprange: Grid step for range dimension (into SLC/SAR geometry). + :param gridstepazimut: Grid step for azimut dimension (into SLC/SAR geometry). + :param threshold: Threshold for correlation rate. + :param gap: Maximum difference between input grid values and mean values. + :param advantage: Give an advantage to DEM or Correlation Grid. + :param outPath: Output Deformation Grid Vector Image (Shift_ran, Shift_azi). + + :type indem: string + :type insarmaster: string + :type insarslave: string + :type inmlmaster: string + :type inmlslave: string + :type indemprojmaster: string + :type indemprojslave: string + :type mlran: int + :type mlazi: int + :type gridsteprange: int + :type gridstepazimut: int + :type threshold: float + :type gap: int + :type advantage: string + :type outPath: string + + :return: None + """ + appFineDeformationGrid = otb.Registry.CreateApplication("SARFineDeformationGrid") + appFineDeformationGrid.SetParameterString("indem", indem) + appFineDeformationGrid.SetParameterString("insarmaster", insarmaster) + appFineDeformationGrid.SetParameterString("insarslave", insarslave) + appFineDeformationGrid.SetParameterString("inmlmaster", inmlmaster) + appFineDeformationGrid.SetParameterString("inmlslave", inmlslave) + appFineDeformationGrid.SetParameterString("indemprojmaster", indemprojmaster) + appFineDeformationGrid.SetParameterString("indemprojslave", indemprojslave) + appFineDeformationGrid.SetParameterInt("mlran", mlran) + appFineDeformationGrid.SetParameterInt("mlazi", mlazi) + appFineDeformationGrid.SetParameterInt("gridsteprange", gridsteprange) + appFineDeformationGrid.SetParameterInt("gridstepazimut", gridstepazimut) + appFineDeformationGrid.SetParameterFloat("threshold", threshold) + appFineDeformationGrid.SetParameterFloat("gap", gap) + appFineDeformationGrid.SetParameterString("advantage", advantage) + appFineDeformationGrid.SetParameterString("out", outPath) + appFineDeformationGrid.SetParameterString("ram", ram) + appFineDeformationGrid.ExecuteAndWriteOutput() + + +def coRegistration(insarmaster, insarslave, ingrid, gridsteprange, gridstepazimut, doppler0, nbramps, outPath, ram="4000"): + """ + Launch SARCoRegistration application (from DiapOTB). + + This application does the coregistration between + two SAR images thanks to a deformation grid. + + :param insarmaster: SAR Master Image to extract SAR geometry. + :param insarslave: SAR Slave Image to extract SAR geometry. + :param ingrid: Input deformation grid. + :param gridsteprange: Grid step for range dimension (into SLC/SAR geometry). + :param gridstepazimut: Grid step for azimut dimension (into SLC/SAR geometry). + :param doppler0: Doppler 0 in azimut for SAR Slave Image. + :param nbramps : Number of ramps for coregistration + :param outPath: Coregistrated slave image into master geometry. + + :type insarmaster: string + :type insarslave: string + :type ingrid: string + :type gridsteprange: int + :type gridstepazimut: int + :type doppler0: int + :type nbramps: int + :type outPath: string + + :return: None + """ + appCoRegistration = otb.Registry.CreateApplication("SARCoRegistration") + appCoRegistration.SetParameterString("insarmaster", insarmaster) + appCoRegistration.SetParameterString("insarslave", insarslave) + appCoRegistration.SetParameterString("ingrid", ingrid) + appCoRegistration.SetParameterInt("gridsteprange", gridsteprange) + appCoRegistration.SetParameterInt("gridstepazimut", gridstepazimut) + appCoRegistration.SetParameterFloat("doppler0", doppler0) + appCoRegistration.SetParameterInt("sizetiles", 50) + appCoRegistration.SetParameterInt("margin", 7) + appCoRegistration.SetParameterInt("nbramps", nbramps) + appCoRegistration.SetParameterString("ram", ram) + appCoRegistration.SetParameterString("out", outPath) + appCoRegistration.ExecuteAndWriteOutput() + +def deramp(inImg, reramp, shift, ingrid, inslave, gridsteprange, gridstepazimut, outPath, ram="4000"): + """ + Launch SARDeramp application (from DiapOTB). + + This application does the deramping or reramping of S1 Iw burst. + + :param inImg: Input burst (from S1 Iw product). + :param reramp: Activate Reramping mode + :param shift: Activate Shift mode (To reramp/deramp accordingly a deformation grid) + :param ingrid: Input deformation grid. + :param insarslave: SAR Slave Image to extract SAR geometry. + :param gridsteprange: Grid step for range dimension (into SLC/SAR geometry). + :param gridstepazimut: Grid step for azimut dimension (into SLC/SAR geometry). + :param outPath: Coregistrated slave image into master geometry. + + :type inImg: string + :type reramp: bool + :type shift: bool + :type ingrid: string + :type inslave: string + :type gridsteprange: int + :type gridstepazimut: int + :type outPath: string + + :return: None + """ + appDerampSlave = otb.Registry.CreateApplication("SARDeramp") + appDerampSlave.SetParameterString("in", inImg) + appDerampSlave.SetParameterString("out", outPath) + if fUtils.str2bool(reramp) : + appDerampSlave.SetParameterString("reramp", reramp) + if fUtils.str2bool(shift) : + appDerampSlave.SetParameterString("shift", shift) + appDerampSlave.SetParameterString("ingrid", ingrid) + appDerampSlave.SetParameterInt("gridsteprange", gridsteprange) + appDerampSlave.SetParameterInt("gridstepazimut", gridstepazimut) + appDerampSlave.SetParameterString("inslave", inslave) + appDerampSlave.SetParameterString("ram", ram) + appDerampSlave.ExecuteAndWriteOutput() + +def interferogram(insarmaster, incoregistratedslave, insarslave, ingrid, incartmeanmaster, gridsteprange, gridstepazimut, mlran, mlazi, gain, outPath, ram="4000"): + """ + Launch SARRobustInterferogram application (from DiapOTB). + + This application estimates robust interferogram with flatenning and averaging. + + :param insarmaster: Input SAR Master image. + :param incoregistratedslave: Input SAR Slave image (Coregistrated image). + :param insarslave: Input SAR Slave image (Metadata). + :param ingrid: Input deformation grid. + :param incartmeanmaster: Input Cartesian Mean Master image. + :param gridsteprange: Grid step for range dimension (into SLC/SAR geometry). + :param gridstepazimut: Grid step for azimut dimension (into SLC/SAR geometry). + :param gain: Gain to apply for amplitude estimation. + :param outPath: Output interferogram (ML geometry). + + :type insarmaster: string + :type incoregistratedslave: string + :type insarslave: string + :type ingrid: string + :type incartmeanmaster: string + :type gridsteprange: int + :type gridstepazimut: int + :type gain: float + :type outPath: string + + :return: None + """ + appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram") + appInterferogram.SetParameterString("insarmaster", insarmaster) + appInterferogram.SetParameterString("incoregistratedslave", incoregistratedslave) + appInterferogram.SetParameterString("insarslave", insarslave) + appInterferogram.SetParameterString("ingrid", ingrid) + appInterferogram.SetParameterString("incartmeanmaster", incartmeanmaster) + appInterferogram.SetParameterInt("mlran", mlran) + appInterferogram.SetParameterInt("mlazi", mlazi) + appInterferogram.SetParameterInt("marginran", 1) + appInterferogram.SetParameterInt("marginazi", 1) + appInterferogram.SetParameterInt("gridsteprange", gridsteprange) + appInterferogram.SetParameterInt("gridstepazimut", gridstepazimut) + appInterferogram.SetParameterFloat("gain", gain) + appInterferogram.SetParameterString("ram", ram) + appInterferogram.SetParameterString("out", outPath) + appInterferogram.ExecuteAndWriteOutput() + +def esd(ininterfup, ininterflow, insar, burstIndex, outPath, ram="4000"): + """ + Launch SARESD application (from DiapOTB). + + This application applies esd processing on two consecutives bursts dimension. + + :param ininterfup: First interferogram. + :param ininterflow: Second interferogram. + :param insar: SAR Image to extract SAR geometry. + :param burstIndex: Index of first interferogram. + :param outPath: Output Tmp. + + :type ininterfup: string + :type ininterflow: string + :type insar: string + :type burstIndex: int + :type outPath: string + + :return: Azimut shift to correct phase jump. + :rtype: float + """ + appESD = otb.Registry.CreateApplication("SARESD") + appESD.SetParameterString("ininterfup", ininterfup) + appESD.SetParameterString("ininterflow", ininterflow) + appESD.SetParameterString("insar", insar) + appESD.SetParameterInt("burstindex", burstIndex) + appESD.SetParameterFloat("threshold", 0.3) + appESD.SetParameterInt("mlazi", 1) + appESD.SetParameterString("ram", ram) + appESD.SetParameterString("out", outPath) + appESD.Execute() + + return appESD.GetParameterFloat('azishift') + +def GridOffset(ingrid, offsetazi, outPath, ram="4000"): + """ + Launch SARGridOffset application (from DiapOTB). + + This application applies offsets on a deformation grid for each dimension. + + :param ingrid: Input Grid Vector Image (Shift_ran, Shift_azi). + :param offsetazi: Offset on azimut. + :param outPath: Output Grid (Tmp). + + :type ingrid: string + :type offsetazi: int + :type outPath: string + + :return: None + """ + appGridOff = otb.Registry.CreateApplication("SARGridOffset") + appGridOff.SetParameterString("ingrid", ingrid) + appGridOff.SetParameterFloat("offsetran", 0.) + appGridOff.SetParameterFloat("offsetazi", offsetazi) + appGridOff.SetParameterString("ram", ram) + appGridOff.SetParameterString("out", outPath) + appGridOff.ExecuteAndWriteOutput() + +def concatenate(list_of_Interferograms, insar, firstBurst, outPath, ram="4000"): + """ + Launch SARConcatenateBursts application (from DiapOTB). + + This application performs a burst concatenation and provides a SAR Deburst Image. + It reads the input image list (single bursts) and generates a whole SAR image with deburst operations. + + :param list_of_Interferograms: The list of bursts to concatenate. + :param insar: Raw Sentinel1 IW SLC image, or any extract of such made by OTB (geom file needed). + :param firstBurst: Index for the first required Burst. + :param outPath: The concatenated and debursted output image. + + :type list_of_Interferograms: list + :type insar: string + :type firstBurst: int + :type outPath: string + + :return: None + """ + appCon = otb.Registry.CreateApplication("SARConcatenateBursts") + appCon.SetParameterStringList("il", list_of_Interferograms) + appCon.SetParameterString("insar", insar) + appCon.SetParameterInt("burstindex", firstBurst) + appCon.SetParameterString("out", outPath) + appCon.SetParameterString("ram", ram) + appCon.ExecuteAndWriteOutput() + diff --git a/python_src/utils/__init__.py b/python_src/utils/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python_src/utils/func_utils.py b/python_src/utils/func_utils.py new file mode 100644 index 0000000..45149e6 --- /dev/null +++ b/python_src/utils/func_utils.py @@ -0,0 +1,235 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +""" +Hello func_utils +""" + + +# Imports +import logging +import json +from jsonschema import validate +import os +import sys +import argparse +import re + +import otbApplication as otb + +# Streamer to our log file +class StreamToLogger(object): + """ + Fake file-like stream object that redirects writes to a logger instance. + """ + def __init__(self, logger, log_level=logging.INFO): + self.logger = logger + self.log_level = log_level + + def write(self, buf): + for line in buf.rstrip().splitlines(): + self.logger.log(self.log_level, line.rstrip()) + + def flush(self): + for handler in self.logger.handlers: + handler.flush() + +# Global variables for logger and std_out +logger = logging.getLogger(__name__) +LogFormatter = logging.Formatter('%(filename)s :: %(levelname)s :: %(message)s') +s1 = StreamToLogger(logger, logging.INFO) +stdout_save = s1 + +### Functions for logger ### +def init_logger(): + + logger.setLevel(logging.INFO) + + # Create console handler with a high log level (warning level) + stream_handler = logging.StreamHandler() + stream_handler.setLevel(logging.WARNING) + + # Add Handlers + logger.addHandler(stream_handler) + + +def init_fileLog(output_dir): + # File handler for the logger + # Create file handler which logs even info messages (used as stdout redirection) + file_handler = logging.FileHandler(os.path.join(output_dir, 'info.log'), 'a') + file_handler.setLevel(logging.INFO) + file_handler.setFormatter(LogFormatter) + + # Add Handlers + logger.addHandler(file_handler) + + # Redirect stdout and stderr to logger + stdout_saveWrite = sys.stdout.write # Save stdout.write to print some info into the console + stdout_saveFlush = sys.stdout.flush # Save stdout.flush to print some info into the console + sys.stdout.write = s1.write # Replace stdout.write by our StreamToLogger + sys.stdout.flush = s1.flush # Replace stdout.flush by our StreamToLogger + #stdout_save = s1 # Different object + stdout_save.write = stdout_saveWrite # Restore stdout.write into stdout_save + stdout_save.flush = stdout_saveFlush # Restore stdout.write into stdout_save + + +def log(level, msg) : + logger.log(level, msg) + + +def printOnStd(msg): + # Print on stdout_save aka the original/true stdout + print(msg, file=stdout_save) + + + +### Functions for configFile (json) ### +def validate_json(json, schema): + # Read and Load the configuration file + try: + validate(json, schema) + except Exception as valid_err: + print("Invalid JSON: {}".format(valid_err)) + return False + else: + # Realise votre travail + print("Valid JSON") + return True + +def load_configfile(configfile, mode="Others"): + + # Empty dictionary + dataConfig = {} + + # Read and Load the configuration file + try: + with open(configfile, 'r') as f: + dataConfig = json.load(f) + + except Exception as err: + logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=configfile)) + quit() + + schema_json = "schema_S1SM.json" + if mode == "S1_IW": + schema_json = "schema_S1IW.json" + + # Load schema (into DiapOTB install) + diapOTB_install = os.getenv('DIAPOTB_INSTALL') + if diapOTB_install is not None and os.path.exists(diapOTB_install): + schemas_path = os.path.join(diapOTB_install, "json_schemas") + if os.path.exists(schemas_path): + schema_file = os.path.join(schemas_path, schema_json) + + try: + with open(schema_file, "r") as sch: + dataSchema = json.load(sch) + except Exception as err: + logger.critical("Impossible to read or load JSON configuration file : {err}. Check its path and content.".format(err=schema_S1SM)) + quit() + + # Check Json file + jsonIsValid = validate_json(dataConfig, dataSchema) + + if not jsonIsValid : + logger.critical("Error, provided config file does not match requirements") + quit() + + return dataConfig + + + +### Functions for metadata ### +def getImageKWL(Img): + # Retrieve some information about our input images with ReadImageInfo application + ReadImageInfo = otb.Registry.CreateApplication("ReadImageInfo") + ReadImageInfo.SetParameterString("in", Img) + ReadImageInfo.SetParameterString("keywordlist", "true") + ReadImageInfo.Execute() + + keyword = ReadImageInfo.GetParameterString("keyword") + keywordlist = ReadImageInfo.GetParameterString("keyword").split("\n") + keywordlist = filter(None, keywordlist) + dictKWL = { i.split(':')[0] : re.sub(r"[\n\t\s]*", "", i.split(':')[1]) for i in keywordlist } + + return dictKWL + + +def getDEMInformation(dem): + # Get information about DEM (spacing, size ..) + ReadDEMInfo = otb.Registry.CreateApplication("ReadImageInfo") + ReadDEMInfo.SetParameterString("in", dem) + ReadDEMInfo.SetParameterString("keywordlist", "true") + ReadDEMInfo.Execute() + + dictDEMInformation = {} + dictDEMInformation['spacingXDEM'] = ReadDEMInfo.GetParameterFloat("spacingx") + dictDEMInformation['estimatedGroundSpacingXDEM'] = ReadDEMInfo.GetParameterFloat("estimatedgroundspacingx") + dictDEMInformation['spacingYDEM'] = ReadDEMInfo.GetParameterFloat("spacingy") + dictDEMInformation['estimatedGroundSpacingYDEM'] = ReadDEMInfo.GetParameterFloat("estimatedgroundspacingy") + + return dictDEMInformation + + + +### Functions for some checks ### +def selectBurst(dictMaster, dictSlave, firstBurst, lastBurst, nbBurstSlave): + + key1Burst = "support_data.geom.bursts.burst[" + + # Initialize the output lists (empty lists) + validBurstMaster = [] + validBurstSlave = [] + + # Loop on Master bursts + for id_B in range(firstBurst, lastBurst+1): + keyBurstMaster = key1Burst + str(id_B) + "].azimuth_anx_time" + + # Get the anx time for the current burst (into Master Image) + anxMaster = float(dictMaster[keyBurstMaster]) + + # Loop on slave bursts to find the closest anx time + minDiff = 200 + id_B_save = id_B + for id_B_slave in range(0, nbBurstSlave): + keyBurstSlave = key1Burst + str(id_B_slave) + "].azimuth_anx_time" + + # Get anx time for slave burst + anxSlave = float(dictSlave[keyBurstSlave]) + + # Comparaison between master and slave + diff = abs(anxMaster - anxSlave) + + if minDiff > diff : + minDiff = diff + id_B_save = id_B_slave + + # Check if difference between the anx time is valid (must be inferior to 1) + if minDiff < 1. : + # Fill lists with master Burst_id and the selected slave burst_id + validBurstMaster.append(id_B) + validBurstSlave.append(id_B_save) + + return validBurstMaster, validBurstSlave + +def str2bool(v): + return v.lower() in ("yes", "true", "t", "1") -- GitLab