ENH : New chain for IW processing (part 1 : Pre-processing)

python_src/diapOTB_S1IW.py

+#!/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
+try:
+    from ConfigParser import SafeConfigParser as Config2Parser
+except ImportError:
+    from configparser import ConfigParser
+import os
+import sys
+import argparse
+import re
+
+import otbApplication as otb
+
+logger = logging.getLogger(__name__)
+
+# ConfigParser functions
+# Get a map of sections 
+def ConfigSectionMap(config, section):
+    dict1 = {}
+    options = config.options(section)
+  
+    for option in options:
+        try:
+            dict1[option] = config.get(section, option)
+            if dict1[option] == -1:
+                DebugPrint("skip: %s" % option)
+        except:
+            print("exception on %s!" % option)
+            dict1[option] = None
+    
+    return dict1
+
+# Check the configuration file
+def ConfigCheck(config):
+    configOk = True
+    # List for configuration
+    listOfSections = ["Global", "Pre_Processing", "Ground", "DIn_SAR"]
+    globalList = ["Master_Image_Path", "Slave_Image_Path", "DEM_Path", "output_dir", "burst_index"]
+    preProcessingList = ["ML_range", "ML_azimut", "ML_gain", "doppler_file"]
+    GroundList = []
+    dInSARList = []
+
+    dictOfSections = {'Global': globalList, 'Pre_Processing': preProcessingList, 'Ground' : GroundList, 'DIn_SAR' : dInSARList}
+
+    # For each elt of listOfSections
+    for section in listOfSections :
+        # Check if section exists into config
+        if config.has_section(section) :
+            # For each elt of lists of dictOfSections
+            for option in dictOfSections[section] :
+                # Check if option exists into section of config
+                if not config.has_option(section, option) :
+                    print(option + " is missing for the section " + section + " of the configuration file")
+                    configOk = False
+
+        else :
+            print(section + " is missing into the configuration file")
+            configOk = False
+        
+
+    return configOk
+
+
+# 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 the configuration file
+    Config = ConfigParser()
+    Config.read(args.configfile)
+    configOk = ConfigCheck(Config)
+    
+    # If some values are missing, quit the application
+    if not configOk :
+        quit()
+        
+    # Get elements from configuration file
+    # Global
+    master_Image = ConfigSectionMap(Config, "Global")['master_image_path']
+    slave_Image =  ConfigSectionMap(Config, "Global")['slave_image_path']
+    dem =  ConfigSectionMap(Config, "Global")['dem_path']
+    output_dir = ConfigSectionMap(Config, "Global")['output_dir']
+    burst_index = ConfigSectionMap(Config, "Global")['burst_index']
+    
+    # Pre_Processing
+    ml_range = int(ConfigSectionMap(Config, "Pre_Processing")['ml_range'])
+    ml_azimut = int(ConfigSectionMap(Config, "Pre_Processing")['ml_azimut'])
+    ml_gain = float(ConfigSectionMap(Config, "Pre_Processing")['ml_gain'])
+    dop_file = ConfigSectionMap(Config, "Pre_Processing")['doppler_file']
+    # Ground
+    
+    # DIn_SAR
+    # geoGrid_gridstep_range = int(ConfigSectionMap(Config, "DIn_SAR")['gridstep_range'])
+    # geoGrid_gridstep_azimut = int(ConfigSectionMap(Config, "DIn_SAR")['gridstep_azimut'])
+    # geoGrid_threshold = float(ConfigSectionMap(Config, "DIn_SAR")['grid_threshold'])
+    # geoGrid_gap = float(ConfigSectionMap(Config, "DIn_SAR")['grid_gap'])
+    # ml_geoGrid_range = ml_range
+    # ml_geoGrid_azimut = ml_azimut
+    # gain_interfero = float(ConfigSectionMap(Config, "DIn_SAR")['interferogram_gain'])
+    # activateOrthoInterferogram = str2bool(ConfigSectionMap(Config, "DIn_SAR")['interferogram_ortho'])
+
+    # if (geoGrid_threshold < 0) or (geoGrid_threshold > 1) :
+    #     print("Wrong Threshold for fine deformation grid")
+
+
+    # Check if images 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 regex for Burst index
+    burstIndexOk = bool(re.match(r'^[\-0-9]+$', burst_index))
+    firstBurst = 0
+    lastBurst = 0
+
+    if burstIndexOk :
+        burstList = burst_index.split('-');
+        burstList = [int(i) for i in burstList]
+        
+        if len(burstList) != 2 :
+            burstIndexOk = False
+        else :
+            firstBurst = min(burstList)
+            lastBurst = max(burstList)
+
+    if not burstIndexOk :
+        print("Wrong Burst Index Format (for instance 0-5)")
+        quit()
+
+    
+    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'])])
+
+    if minNbBurst < firstBurst or minNbBurst < lastBurst :
+        print ("Wrong burst index (superior to number of bursts)")
+        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 = []
+    
+    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("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()
+        
+        ######## 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("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 #############################
+    
+
+
+    
+    ######################## DIn_SAR Chain #############################
+    # ######## SARDEMProjection Application #######
+    # print("\n SARDEMProjection Application \n")
+    # # 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.SetParameterEmpty("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.SetParameterEmpty("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")
+    # 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.SetParameterInt("patchsize", 20)
+    # appFineDeformationGrid.SetParameterInt("subpixel", 10)
+    # appFineDeformationGrid.SetParameterFloat("threshold", geoGrid_threshold)
+    # appFineDeformationGrid.SetParameterFloat("gap", geoGrid_gap)
+    # appFineDeformationGrid.SetParameterString("out", os.path.join(output_dir, fine_grid))
+    # appFineDeformationGrid.SetParameterString("ram", "4000")
+    # appFineDeformationGrid.ExecuteAndWriteOutput()
+
+    
+    # ######## SARCoRegistration Application (changeo step) #######
+    # print("\n SARCoRegistration Application \n")
+    # 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")
+    # # 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")
+    # interferogram = "interferogram.tif"
+    # interferogram_ortho = "interferogram_ortho.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))
+    # if activateOrthoInterferogram :
+    #     appInterferogram.SetParameterEmpty("ortho", True)
+    #     appInterferogram.SetParameterString("indem", dem)
+    #     appInterferogram.SetParameterString("outopt", os.path.join(output_dir, interferogram_ortho))
+    # appInterferogram.ExecuteAndWriteOutput()

python_src/ex_config_S1IW.ini

+[Global]
+Master_Image_Path = image_1.tif
+Slave_Image_Path = image_2.tif
+DEM_Path = ./DEM.hgt
+output_dir = ./output_diapOTB
+burst_index = 0-1
+
+[Pre_Processing]
+ML_range = 3
+ML_azimut = 3
+ML_gain = 0.1	
+doppler_file = dop0.txt
+
+[Ground]
+
+
+
+[DIn_SAR]
\ No newline at end of file