Process Lutan-1 stripmap SLC images using stripmapApp.py (#852)

* Update estimateIono.py

* Update runDispersive.py

* Added Lutan-1 support

* Added revision for Lutan-1

* Added Lutan 1 in SConscript and CMakeLists.txt

* Fixed Lutan-1 name

* Bug fixes for Lutan-1 reader

* Added 128k on the Fortran code

* Added the ability to ingest annotation file orbit info

* Added warning for using orbit info from annotation file

* Update README.md

---------

Co-authored-by: bjmarfito <[email protected]>
Co-authored-by: Zhang Yunjun <[email protected]>

README.md

@@ -22,7 +22,7 @@ utility in a general sense for building other types of software packages.  In
 its InSAR aspect ISCE supports data from many space-borne satellites and one
 air-borne platform.  We continue to increase the number of sensors supported.
 At this time the sensors that are supported are the following: ALOS, ALOS2,
-COSMO_SKYMED, ENVISAT, ERS, KOMPSAT5, RADARSAT1, RADARSAT2, RISAT1, Sentinel1,
+COSMO_SKYMED, ENVISAT, ERS, KOMPSAT5, LUTAN1, RADARSAT1, RADARSAT2, RISAT1, Sentinel1,
 TERRASARX, UAVSAR and SAOCOM1A.
 
 ## Contents

applications/make_raw.py

@@ -226,7 +226,7 @@ def calculateSquint(self):
         self.frame.squintAngle = math.radians(self.squint)
 
     def make_raw(self):
-        from  isceobj.Image import createRawImage, createSlcImage
+        from isceobj.Image import createRawImage, createSlcImage
         self.activateInputPorts()
 
         # Parse the image metadata and extract the image

components/isceobj/Sensor/ALOS2.py

@@ -106,7 +106,7 @@ def __init__(self, name=''):
         self.leaderFile = None
         self.imageFile = None
 
-        #####Soecific doppler functions for ALOS2
+        #####Specific doppler functions for ALOS2
         self.doppler_coeff = None
         self.azfmrate_coeff = None
         self.lineDirection = None

components/isceobj/Sensor/CMakeLists.txt

@@ -41,6 +41,7 @@ set(installfiles
     UAVSAR_RPI.py
     UAVSAR_Stack.py
     SAOCOM_SLC.py
+    Lutan1.py
     )
 
 if(HDF5_FOUND)

components/isceobj/Sensor/Lutan1.py

@@ -0,0 +1,323 @@
+#!/usr/bin/python3
+
+# Reader for Lutan-1 SLC data
+# Used Sentinel1.py and ALOS.py as templates
+# Author: Bryan Marfito, EOS-RS
+
+
+import os
+import glob
+import numpy as np
+import xml.etree.ElementTree as ET
+import datetime
+import isce
+import isceobj
+from isceobj.Planet.Planet import Planet
+from iscesys.Component.Component import Component
+from isceobj.Sensor.Sensor import Sensor
+from isceobj.Scene.Frame import Frame
+from isceobj.Orbit.Orbit import StateVector, Orbit
+from isceobj.Planet.AstronomicalHandbook import Const
+from iscesys.DateTimeUtil.DateTimeUtil import DateTimeUtil as DTUtil
+from isceobj.Orbit.OrbitExtender import OrbitExtender
+from osgeo import gdal
+import warnings
+
+lookMap = { 'RIGHT' : -1,
+            'LEFT' : 1}
+
+# Antenna dimensions 9.8 x 3.4 m
+antennaLength = 9.8
+
+XML = Component.Parameter('xml',
+        public_name = 'xml',
+        default = None,
+        type = str,
+        doc = 'Input XML file')
+
+
+TIFF = Component.Parameter('tiff',
+                            public_name ='tiff',
+                            default = None,
+                            type=str,
+                            doc = 'Input image file')
+
+ORBIT_FILE = Component.Parameter('orbitFile',
+                            public_name ='orbitFile',
+                            default = None,
+                            type=str,
+                            doc = 'Orbit file')
+
+
+class Lutan1(Sensor):
+
+    "Class for Lutan-1 SLC data"
+
+    family = 'l1sm'
+    logging_name = 'isce.sensor.Lutan1'
+
+    parameter_list = (TIFF, ORBIT_FILE) + Sensor.parameter_list
+
+    def __init__(self, name = ''):
+        super(Lutan1,self).__init__(self.__class__.family, name=name)
+        self.frame = Frame()
+        self.frame.configure()
+        self._xml_root = None
+        self.doppler_coeff = None
+
+    def parse(self):
+        xmlFileName = self.tiff[:-4] + "meta.xml"
+        self.xml = xmlFileName
+
+        with open(self.xml, 'r') as fid:
+            xmlstr = fid.read()
+
+        self._xml_root = ET.fromstring(xmlstr)
+        self.populateMetadata()
+        fid.close()
+
+        if self.orbitFile:
+            # Check if orbit file exists or not
+            if os.path.isfile(self.orbitFile) == True:
+                orb = self.extractOrbit()
+                self.frame.orbit.setOrbitSource(os.path.basename(self.orbitFile))
+            else:
+                pass
+        else:
+            warnings.warn("WARNING! No orbit file found. Orbit information from the annotation file is used for processing.")
+            orb = self.extractOrbitFromAnnotation()
+            self.frame.orbit.setOrbitSource(os.path.basename(self.xml))
+            self.frame.orbit.setOrbitSource('Annotation')
+
+        for sv in orb:
+            self.frame.orbit.addStateVector(sv)
+
+    def convertToDateTime(self,string):
+        dt = datetime.datetime.strptime(string,"%Y-%m-%dT%H:%M:%S.%f")
+        return dt
+
+
+    def grab_from_xml(self, path):
+        try:
+            res = self._xml_root.find(path).text
+        except:
+            raise Exception('Tag= %s not found'%(path))
+
+        if res is None:
+            raise Exception('Tag = %s not found'%(path))
+
+        return res
+
+
+    def populateMetadata(self):
+        mission = self.grab_from_xml('generalHeader/mission')
+        polarization = self.grab_from_xml('productInfo/acquisitionInfo/polarisationMode')
+        frequency = float(self.grab_from_xml('instrument/radarParameters/centerFrequency'))
+        passDirection = self.grab_from_xml('productInfo/missionInfo/orbitDirection')
+        rangePixelSize = float(self.grab_from_xml('productInfo/imageDataInfo/imageRaster/columnSpacing'))
+        azimuthPixelSize = float(self.grab_from_xml('productInfo/imageDataInfo/imageRaster/rowSpacing'))
+        rangeSamplingRate = Const.c/(2.0*rangePixelSize)
+
+        prf = float(self.grab_from_xml('instrument/settings/settingRecord/PRF'))
+        lines = int(self.grab_from_xml('productInfo/imageDataInfo/imageRaster/numberOfRows'))
+        samples = int(self.grab_from_xml('productInfo/imageDataInfo/imageRaster/numberOfColumns'))
+
+        startingRange = float(self.grab_from_xml('productInfo/sceneInfo/rangeTime/firstPixel'))*Const.c/2.0
+        #slantRange = float(self.grab_from_xml('productSpecific/complexImageInfo/'))
+        incidenceAngle = float(self.grab_from_xml('productInfo/sceneInfo/sceneCenterCoord/incidenceAngle'))
+        dataStartTime = self.convertToDateTime(self.grab_from_xml('productInfo/sceneInfo/start/timeUTC'))
+        dataStopTime = self.convertToDateTime(self.grab_from_xml('productInfo/sceneInfo/stop/timeUTC'))
+        pulseLength = float(self.grab_from_xml('processing/processingParameter/rangeCompression/chirps/referenceChirp/pulseLength'))
+        pulseBandwidth = float(self.grab_from_xml('processing/processingParameter/rangeCompression/chirps/referenceChirp/pulseBandwidth'))
+        chirpSlope = pulseBandwidth/pulseLength
+
+        if self.grab_from_xml('processing/processingParameter/rangeCompression/chirps/referenceChirp/chirpSlope') == "DOWN":
+            chirpSlope = -1.0 * chirpSlope
+        else:
+            pass
+
+        # Check for satellite's look direction
+        if self.grab_from_xml('productInfo/acquisitionInfo/lookDirection') == "LEFT":
+            lookSide = lookMap['LEFT']
+            print("Look direction: LEFT")
+        else:
+            lookSide = lookMap['RIGHT']
+            print("Look direction: RIGHT")
+
+        processingFacility = self.grab_from_xml('productInfo/generationInfo/level1ProcessingFacility')
+
+        # Platform parameters
+        platform = self.frame.getInstrument().getPlatform()
+        platform.setPlanet(Planet(pname='Earth'))
+        platform.setMission(mission)
+        platform.setPointingDirection(lookSide)
+        platform.setAntennaLength(antennaLength)
+
+        # Instrument parameters
+        instrument = self.frame.getInstrument()
+        instrument.setRadarFrequency(frequency)
+        instrument.setPulseRepetitionFrequency(prf)
+        instrument.setPulseLength(pulseLength)
+        instrument.setChirpSlope(chirpSlope)
+        instrument.setIncidenceAngle(incidenceAngle)
+        instrument.setRangePixelSize(rangePixelSize)
+        instrument.setRangeSamplingRate(rangeSamplingRate)
+        instrument.setPulseLength(pulseLength)
+
+        # Frame parameters
+        self.frame.setSensingStart(dataStartTime)
+        self.frame.setSensingStop(dataStopTime)
+        self.frame.setProcessingFacility(processingFacility)
+
+        # Two-way travel time 
+        diffTime = DTUtil.timeDeltaToSeconds(dataStopTime - dataStartTime) / 2.0
+        sensingMid = dataStartTime + datetime.timedelta(microseconds=int(diffTime*1e6))
+        self.frame.setSensingMid(sensingMid)
+        self.frame.setPassDirection(passDirection)
+        self.frame.setPolarization(polarization)
+        self.frame.setStartingRange(startingRange)
+        self.frame.setFarRange(startingRange +  (samples - 1) * rangePixelSize)
+        self.frame.setNumberOfLines(lines)
+        self.frame.setNumberOfSamples(samples)
+
+        return
+
+
+    def extractOrbit(self):
+
+        '''
+        Extract orbit information from the orbit file
+        '''
+
+        try:
+            fp = open(self.orbitFile, 'r')
+        except IOError as strerr:
+            print("IOError: %s" % strerr)
+
+        _xml_root = ET.ElementTree(file=fp).getroot()
+        node = _xml_root.find('Data_Block/List_of_OSVs')
+
+        orb = Orbit()
+        orb.configure()
+
+        # I based the margin on the data that I have.
+        # Lutan-1 position and velocity sampling frequency is 1 Hz
+        margin = datetime.timedelta(seconds=1.0)
+        tstart = self.frame.getSensingStart() - margin
+        tend = self.frame.getSensingStop() + margin
+
+        for child in node:
+            timestamp = self.convertToDateTime(child.find('UTC').text)
+            if (timestamp >= tstart) and (timestamp <= tend):
+                pos = []
+                vel = []
+                for tag in ['VX', 'VY', 'VZ']:
+                    vel.append(float(child.find(tag).text))
+
+                for tag in ['X', 'Y', 'Z']:
+                    pos.append(float(child.find(tag).text))
+
+                vec = StateVector()
+                vec.setTime(timestamp)
+                vec.setPosition(pos)
+                vec.setVelocity(vel)
+                orb.addStateVector(vec)
+
+        fp.close()
+
+        return orb
+
+    def extractOrbitFromAnnotation(self):
+
+        '''
+        Extract orbit information from xml annotation
+        WARNING! Only use this method if orbit file is not available
+        '''
+
+        try:
+            fp = open(self.xml, 'r')
+        except IOError as strerr:
+            print("IOError: %s" % strerr)
+
+        _xml_root = ET.ElementTree(file=fp).getroot()
+        node = _xml_root.find('platform/orbit')
+        countNode = len(list(_xml_root.find('platform/orbit')))
+
+        frameOrbit = Orbit()
+        frameOrbit.setOrbitSource('Header')
+        margin = datetime.timedelta(seconds=1.0)
+        tstart = self.frame.getSensingStart() - margin
+        tend = self.frame.getSensingStop() + margin
+
+        for k in range(1,countNode):
+            timestamp = self.convertToDateTime(node.find('stateVec[{}]/timeUTC'.format(k)).text)
+            if (timestamp >= tstart) and (timestamp <= tend):
+                pos = [float(node.find('stateVec[{}]/posX'.format(k)).text), float(node.find('stateVec[{}]/posY'.format(k)).text), float(node.find('stateVec[{}]/posZ'.format(k)).text)]
+                vel = [float(node.find('stateVec[{}]/velX'.format(k)).text), float(node.find('stateVec[{}]/velY'.format(k)).text), float(node.find('stateVec[{}]/velZ'.format(k)).text)]
+
+                vec = StateVector()
+                vec.setTime(timestamp)
+                vec.setPosition(pos)
+                vec.setVelocity(vel)
+                frameOrbit.addStateVector(vec)
+
+        fp.close()
+        return frameOrbit
+
+    def extractImage(self):
+        self.parse()
+        width = self.frame.getNumberOfSamples()
+        lgth = self.frame.getNumberOfLines()
+        src = gdal.Open(self.tiff.strip(), gdal.GA_ReadOnly)
+
+        # Band 1 as real and band 2 as imaginary numbers
+        # Confirmed by Zhang Yunjun
+        band1 = src.GetRasterBand(1)
+        band2 = src.GetRasterBand(2)
+        cJ = np.complex64(1.0j)
+
+        fid = open(self.output, 'wb')
+        for ii in range(lgth):
+            # Combine the real and imaginary to make
+            # them in to complex numbers
+            real = band1.ReadAsArray(0,ii,width,1)
+            imag = band2.ReadAsArray(0,ii,width,1)
+            # Data becomes np.complex128 after combining them
+            data = real + (cJ * imag)
+            data.tofile(fid)
+
+        fid.close()
+        real = None
+        imag = None
+        src = None
+        band1 = None
+        band2 = None
+
+        ####
+        slcImage = isceobj.createSlcImage()
+        slcImage.setByteOrder('l')
+        slcImage.setFilename(self.output)
+        slcImage.setAccessMode('read')
+        slcImage.setWidth(self.frame.getNumberOfSamples())
+        slcImage.setLength(self.frame.getNumberOfLines())
+        slcImage.setXmin(0)
+        slcImage.setXmax(self.frame.getNumberOfSamples())
+        self.frame.setImage(slcImage)
+
+    def extractDoppler(self):
+        '''
+        Set doppler values to zero since the metadata doppler values are unreliable.
+        Also, the SLC images are zero doppler.
+        '''
+        dop = [0., 0., 0.]
+
+        ####For insarApp
+        quadratic = {}
+        quadratic['a'] = dop[0] / self.frame.getInstrument().getPulseRepetitionFrequency()
+        quadratic['b'] = 0.
+        quadratic['c'] = 0.
+
+        print("Average doppler: ", dop)
+        self.frame._dopplerVsPixel = dop
+
+        return quadratic

components/isceobj/Sensor/SConscript

@@ -28,7 +28,7 @@ listFiles = ['ALOS.py','CEOS.py','COSMO_SkyMed.py','COSMO_SkyMed_SLC.py',
              'UAVSAR_Polsar.py', 'ERS_EnviSAT.py', 'ICEYE_SLC.py',
 	     'ALOS2.py', 'ERS_SLC.py', 'ALOS_SLC.py', 'EnviSAT_SLC.py',
              'ERS_EnviSAT_SLC.py', 'SICD_RGZERO.py','UAVSAR_HDF5_SLC.py',
-             'SAOCOM_SLC.py','__init__.py']
+             'SAOCOM_SLC.py','__init__.py','Lutan1.py']
 
 helpList,installHelp = envSensor['HELP_BUILDER'](envSensor,'__init__.py',install)
 envSensor.Install(installHelp,helpList)