From 6277960ccdb1106844d65e73fc1cebcc82361701 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ga=C3=ABlle=20USSEGLIO?= <gaelle.usseglio@cnes.fr>
Date: Tue, 19 Feb 2019 12:57:00 +0000
Subject: [PATCH] ENH : S1_IW CoRegistration Chain
---
python_src/coRegistation_S1IW.py | 709 +++++++++++++++++++++++++++++++
1 file changed, 709 insertions(+)
create mode 100644 python_src/coRegistation_S1IW.py
diff --git a/python_src/coRegistation_S1IW.py b/python_src/coRegistation_S1IW.py
new file mode 100644
index 0000000..1c7b6b9
--- /dev/null
+++ b/python_src/coRegistation_S1IW.py
@@ -0,0 +1,709 @@
+#!/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
+
+logger = logging.getLogger(__name__)
+
+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")
+
+
+# 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)
+
+ # Read and Load the configuration file
+ try:
+ with open(args.configfile, 'r') as f:
+ dataConfig = json.load(f)
+
+ except Exception as err:
+ print("Impossible to read or load JSON configuration file")
+ 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:
+ print("Impossible to read or load JSON schema file")
+ quit()
+
+ # Check Json file
+ jsonIsValid = validate_json(dataConfig, dataSchema)
+
+ if not jsonIsValid :
+ 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']
+ activateOrthoInterferogram = dict_DInSAR['parameter']['Interferogram_ortho']
+
+ esd_NbIter = 0
+
+ if 'ESD_iter' in dict_DInSAR['parameter']:
+ esd_NbIter = dict_DInSAR['parameter']['ESD_iter']
+
+ print("esd_NbIter" + str(esd_NbIter))
+
+ # Check Threshold
+ if (geoGrid_threshold < 0) or (geoGrid_threshold > 1) :
+ print("Wrong Threshold for fine deformation grid")
+ geoGrid_threshold = 0.3
+
+ # Check if images/dem exist
+ if not os.path.exists(master_Image) :
+ print(master_Image + " does not exists")
+ quit()
+ if not os.path.exists(slave_Image) :
+ print(slave_Image + " does not exists")
+ quit()
+ if not os.path.exists(dem) :
+ print(dem + " does not exists")
+ 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")
+
+ # 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 :
+ print ("Upgrade your geom file")
+ quit()
+
+ # 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:
+ print("Wrong burst index")
+ quit()
+
+ if minNbBurst < firstBurst or minNbBurst < lastBurst or lastBurst < 0 or firstBurst < 0 :
+ print ("Wrong burst index")
+ quit()
+
+ # Create directory for each burst
+ for burstId in range(firstBurst, lastBurst+1):
+ if not os.path.exists(os.path.join(output_dir, "burst" + str(burstId))):
+ os.makedirs(os.path.join(output_dir, "burst" + str(burstId)))
+
+ ####################### 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")
+
+ dop0Master = []
+ list_of_Deramps_Master = []
+
+ for burstId in range(firstBurst, lastBurst+1):
+
+ print("\n BurstId = " + str(burstId) + "\n")
+
+ burst_dir = os.path.join(output_dir, "burst" + str(burstId))
+
+ ######## SARBurstExtraction Application #######
+ print("\n Burst Extraction Application \n")
+ 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")
+ 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()
+
+ list_of_Deramps_Master.append(os.path.join(burst_dir, burstDerampM))
+
+ ######## SARDoppler Application #######
+ print("\n Doppler Application \n")
+ 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")
+ 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")
+
+ dop0Slave = []
+
+ for burstId in range(firstBurst, lastBurst+1):
+
+ print("\n BurstId = " + str(burstId) + "\n")
+
+ burst_dir = os.path.join(output_dir, "burst" + str(burstId))
+
+ ######## SARBurstExtraction Application #######
+ print("\n Burst Extraction Application \n")
+ burstS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + ".tif"
+ appBurstExtractionSlave = otb.Registry.CreateApplication("SARBurstExtraction")
+ appBurstExtractionSlave.SetParameterString("in", slave_Image)
+ appBurstExtractionSlave.SetParameterString("out", os.path.join(burst_dir, burstS))
+ appBurstExtractionSlave.SetParameterInt("burstindex", burstId)
+ #appBurstExtractionSlave.SetParameterString("allpixels", "true")
+ appBurstExtractionSlave.SetParameterString("ram", "4000")
+ appBurstExtractionSlave.ExecuteAndWriteOutput()
+
+ ######## SARDeramp Application #######
+ # print("\n Deramping Application \n")
+ burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_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")
+ 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) + "\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")
+ slave_Image_ML = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_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")
+
+ gainMaster = []
+ directionDEMMaster = []
+
+ for burstId in range(firstBurst, lastBurst+1):
+
+ print("\n BurstId = " + str(burstId) + "\n")
+
+ 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")
+ 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")
+ 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")
+
+ gainSlave = []
+ directionDEMSlave = []
+
+ for burstId in range(firstBurst, lastBurst+1):
+
+ print("\n BurstId = " + str(burstId) + "\n")
+
+ burst_dir = os.path.join(output_dir, "burst" + str(burstId))
+
+ burstDerampS = os.path.splitext(slave_Image_base)[0] + "_burst" + str(burstId) + "_deramp" + ".tif"
+
+ demProj_Slave = "demProj" + "_burst" + str(burstId) + "_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")
+
+ counter = 0
+ list_of_Interferograms = []
+ list_of_Grids = []
+ list_of_Coregistrated_Slave = []
+
+ for burstId in range(firstBurst, lastBurst+1):
+
+ print("\n BurstId = " + str(burstId) + "\n")
+
+ 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) + ".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) + "_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) + "_ml" + str(ml_azimut) + str(ml_range) + ".tif"
+
+
+ demProj_Master = "demProj" + "_burst" + str(burstId) +"_Master.tif"
+ demProj_Slave = "demProj" + "_burst" + str(burstId) + "_Slave.tif"
+ master_cartesian_mean = "CartMeanMaster" + "_burst" + str(burstId) + ".tif"
+
+
+ ######## Step 1 : SARFineDeformationGrid #######
+ ######## SARFineDeformationGrid Application (geo_grid step) #######
+ print("\n SARFineDeformationGrid Application \n")
+ 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.SetParameterInt("patchsize", 20)
+ appFineDeformationGrid.SetParameterInt("subpixel", 10)
+ appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold)
+ appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap)
+ 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")
+ 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()
+
+ list_of_Coregistrated_Slave.append(os.path.join(burst_dir, slave_Image_CoRe))
+
+ ######## 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"
+ 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")
+ interferogram_path = os.path.join(burst_dir, "interferogram" + "_burst" + str(burstId) + ".tif")
+ interferogram_ortho_path = os.path.join(burst_dir, "interferogram_ortho" + "_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")
+ interferogram_ortho_path = os.path.join(esd_dir, "interferogram_ortho" + "_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)
+ if activateOrthoInterferogram :
+ appInterferogram.SetParameterString("ortho", "true")
+ appInterferogram.SetParameterString("indem", dem)
+ appInterferogram.SetParameterString("outopt", interferogram_ortho_path)
+ appInterferogram.ExecuteAndWriteOutput()
+
+ list_of_Interferograms.append(interferogram_path)
+ counter = counter + 1
+
+
+ ######## Step 4 : ESD Loop (SARESD + SARGridOffset + Step 3) #######
+ if esd_NbIter > 0 :
+
+ azimut_shift_esd = []
+ azimut_shift_esd_global = [0.] * (lastBurst-firstBurst)
+
+ for iter_esd in range(1, esd_NbIter+1) :
+
+ # Clear all azimut shifts
+ azimut_shift_esd[:] = []
+
+ for burstId in range(firstBurst, lastBurst+1):
+
+ 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")
+ 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[:] = []
+ list_of_Coregistrated_Slave[:] = []
+
+ for burstId in range(firstBurst, lastBurst+1):
+
+ # 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) + ".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) + "_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"
+ 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.ExecuteAndWriteOutput()
+
+ list_of_Coregistrated_Slave.append(os.path.join(esd_dir, slave_Image_CoRe))
+
+ ######## 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"
+ appDerampSlave = otb.Registry.CreateApplication("SARDeramp")
+ appDerampSlave.SetParameterString("in", os.path.join(esd_dir, slave_Image_CoRe))
+ 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")
+ interferogram_path = os.path.join(esd_dir, "interferogram" + "_burst" + str(burstId) + "_iter" + str(iter_esd) + ".tif")
+ interferogram_ortho_path = os.path.join(esd_dir, "interferogram_ortho" + "_burst" + str(burstId) + "_iter" + str(iter_esd) + ".tif")
+
+ appInterferogram = otb.Registry.CreateApplication("SARRobustInterferogram")
+ appInterferogram.SetParameterString("insarmaster", os.path.join(burst_dir, burstM))
+ appInterferogram.SetParameterComplexInputImage("incoregistratedslave",
+ appDerampSlave.GetParameterComplexOutputImage("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)
+ if activateOrthoInterferogram :
+ appInterferogram.SetParameterString("ortho", "true")
+ appInterferogram.SetParameterString("indem", dem)
+ appInterferogram.SetParameterString("outopt", interferogram_ortho_path)
+ appInterferogram.ExecuteAndWriteOutput()
+
+ # Store the new interferogram paths
+ list_of_Interferograms.append(interferogram_path)
+
+
+ ######## Step 5 : SARConcatenateBursts #######
+ ######### SARConcatenateBursts (for coregistrated and deramped images) #########
+ print("\n SARConcatenateBursts Application (for coregistrated and deramped images) \n")
+ deramp_master = "master_deramp.tif"
+ appCon1 = otb.Registry.CreateApplication("SARConcatenateBursts")
+ appCon1.SetParameterStringList("il", list_of_Deramps_Master)
+ appCon1.SetParameterString("insar", master_Image)
+ appCon1.SetParameterInt("burstindex", firstBurst)
+ appCon1.SetParameterString("out", os.path.join(output_dir, deramp_master))
+ appCon1.SetParameterString("ram", "4000")
+ appCon1.ExecuteAndWriteOutput()
+
+
+ deramp_slave = "slave_coregistrated_deramp.tif"
+ appCon2 = otb.Registry.CreateApplication("SARConcatenateBursts")
+ appCon2.SetParameterStringList("il", list_of_Coregistrated_Slave)
+ appCon2.SetParameterString("insar", master_Image)
+ appCon2.SetParameterInt("burstindex", firstBurst)
+ appCon2.SetParameterString("out", os.path.join(output_dir, deramp_slave))
+ appCon2.SetParameterString("ram", "4000")
+ appCon2.ExecuteAndWriteOutput()
--
GitLab