Merge pull request #64 from MIRSL/sentinel1_grd_to_RTC

Sentinel1 grd to rtc

README.md

@@ -623,7 +623,23 @@ between three files as follows:
     <property name="OUTPUT">20061231</property>
 </component>
 ```
-
+### rtcApp.xml 
+The inputs are Sentinel GRD zipfiles
+```xml
+<rtcApp>
+    <constant name="dir">/Users/data/sentinel1 </constant>
+    <component name="rtcApp">
+        <property name="posting">20</property>
+        <property name="sensor name">sentinel1</property>
+        <component name="master">
+        <property name="safe">$dir$/rtcApp/data/S1A_IW_GRDH_1SDV_20181221T225104_20181221T225129_025130_02C664_B46C.zip</property>
+        <property name="orbit directory">$dir$/orbits</property>
+        <property name="output directory">$dir$/rtcApp/output</property>
+        <property name="polarization">[VV, VH]</property>
+        </component>
+    </component>
+</rtcApp>
+```
 -----
 
 ## Component Configurability

applications/imageMath.py

@@ -293,6 +293,7 @@ def main(args, files):
 
     #######Determine number of input and output bands
     bandList = []
+    iMath['equations'] = []
     for ii,expr in enumerate(args.equation.split(';')):
 
         #####Now parse the equation to get the file names used
@@ -319,7 +320,11 @@ def main(args, files):
 
     ######Create input memmaps
     for ii,infile in enumerate(fileList):
+        if type(files) == list:
             fstr, files = parseInputFile(infile, files)
+        else:
+            fstr = getattr(files, infile)
+
         logger.debug('Input string for File %d: %s: %s'%(ii, infile, fstr))
 
         if len(fstr.split(';')) > 1:
@@ -341,6 +346,7 @@ def main(args, files):
         if bbox is not None:
             iMath['bboxes'].append(bbox)
 
+    if type(files) == list:
         if len(files):
             raise IOError('Unused input variables set:\n'+ ' '.join(files))
 

applications/rtcApp.py

@@ -155,7 +155,7 @@ NUMBER_RANGE_LOOKS = Application.Parameter('numberRangeLooks',
 )
 
 POSTING = Application.Parameter('posting',
-            public_name='azimuth looks',
+            public_name='posting',
             default = 20.0,
             type = float,
             mandatory = False,
@@ -363,6 +363,7 @@ class GRDSAR(Application):
         self.verifyDEM = RtcProc.createVerifyDEM(self)
         self.multilook = RtcProc.createLooks(self)
         self.runTopo  = RtcProc.createTopo(self)
+        self.runNormalize = RtcProc.createNormalize(self)
 #        self.runGeocode = RtcProc.createGeocode(self)
 
         return None
@@ -392,6 +393,9 @@ class GRDSAR(Application):
         ##Run topo for each bursts
         self.step('topo', func=self.runTopo)
 
+	    ##Run normalize to get gamma0
+        self.step('normalize', func=self.runNormalize)
+
         # Geocode
 #        self.step('geocode', func=self.runGeocode,
 #                args=(self.geocode_list, self.do_unwrap, self.geocode_bbox))
@@ -417,6 +421,9 @@ class GRDSAR(Application):
         ##Run topo for each burst
         self.runTopo()
 	
+	##Run normalize to get gamma0
+        self.runNormalize()
+
         ###Compute covariance
 #        self.runEstimateCovariance()
 

components/isceobj/RtcProc/Factories.py

@@ -46,5 +46,6 @@ createPreprocessor = _factory("runPreprocessor")
 createVerifyDEM = _factory("runVerifyDEM")
 createLooks = _factory("runLooks")
 createTopo = _factory("runTopo")
+createNormalize = _factory("runNormalize")
 #createGeocode = _factory("runGeocode")
 

components/isceobj/RtcProc/RtcProc.py

@@ -69,7 +69,7 @@ INC_FILENAME = Component.Parameter(
 GAMMA0_FILENAME = Component.Parameter(
         'gamma0FileName',
         public_name='Gamma0 backscatter file',
-        default = 'gamma0.rdr',
+        default = 'gamma0.img',
         type = str,
         mandatory = False,
         doc = 'Unmasked gamma0 backscatter file')

components/isceobj/RtcProc/runNormalize.py

@@ -1,4 +1,4 @@
-#
+#!/usr/bin/env python3
 # Author: Piyush Agram
 # Copyright 2016
 #
@@ -6,19 +6,23 @@
 import logging
 import isceobj
 import mroipac
+from .runTopo import filenameWithLooks
+from .runLooks import takeLooks
 import os
+import itertools
 import numpy as np
 from isceobj.Util.decorators import use_api
+from applications import imageMath
 
 logger = logging.getLogger('isce.grdsar.looks')
 
-
+class Dummy:
+    pass
 
 def runNormalize(self):
     '''
     Make sure that a DEM is available for processing the given data.
     '''
-   
     refPol = self._grd.polarizations[0]
     master = self._grd.loadProduct( os.path.join(self._grd.outputFolder, 'beta_{0}.xml'.format(refPol)))
 
@@ -26,17 +30,31 @@ def runNormalize(self):
     azlooks, rglooks = self._grd.getLooks( self.posting, master.groundRangePixelSize, master.azimuthPixelSize, self.numberAzimuthLooks, self.numberRangeLooks)
 
 
-    if (azlooks == 1) and (rglooks == 1):
-        return
-
-    slantRange = False
     for pol in self._grd.polarizations:
-        inname = os.path.join( self._grd.outputFolder, 'beta_{0}.img'.format(pol) )
+        if (azlooks == 1) and (rglooks == 1):
-        takeLooks(inname, azlooks, rglooks)
+            inname = os.path.join( self._grd.outputFolder, 'beta_{0}.img'.format(pol))
+        else:
+            inname = os.path.join( self._grd.outputFolder, filenameWithLooks('beta_{0}.img'.format(pol), azlooks, rglooks))
 
-        if not slantRange:
+        basefolder, output = os.path.split(self._grd.outputFolder)
-            inname = master.slantRangeImage.filename
+        incname = os.path.join(basefolder, self._grd.geometryFolder, self._grd.incFileName)
-            takeLooks(inname, azlooks, rglooks)
+        outname = os.path.join(self._grd.outputFolder, filenameWithLooks('gamma_{0}'.format(pol)+'.img', azlooks, rglooks))
-            slantRange = True
+        maskname = os.path.join(basefolder, self._grd.geometryFolder, self._grd.slMaskFileName)
+
+        args = imageMath.createNamespace()
+        args.equation = 'a*cos(b_0*PI/180.)/cos(b_1*PI/180.) * (c==0)'
+        args.dtype = np.float32
+        args.scheme = 'BIL'
+        args.out = outname
+        #args.debug = True
+
+        files = Dummy()
+        files.a = inname
+        files.b = incname
+        files.c = maskname
+
+
+
+        imageMath.main(args, files)
 
     return

components/isceobj/Sensor/GRD/Sentinel1.py

@@ -261,10 +261,10 @@ class Sentinel1(Component):
 
         self.validateUserInputs()
 
-        if self.xml.startswith('/vsizip'): #Read from zip file
+        if '.zip' in self.xml:
             try:
                 parts = self.xml.split(os.path.sep)
-                zipname = os.path.join(*(parts[2:-3]))
+                zipname = os.path.join('/',*(parts[:-3]))
                 fname = os.path.join(*(parts[-3:]))
 
                 with zipfile.ZipFile(zipname, 'r') as zf:
@@ -283,23 +283,22 @@ class Sentinel1(Component):
         self.populateMetadata()
         self.populateBbox()
 
-        ####Tru and locate an orbit file
+        ####Try and locate an orbit file
         if self.orbitFile is None:
             if self.orbitDir is not None:
                 self.orbitFile = self.findOrbitFile()
-
+                print('Found this orbitfile: %s' %self.orbitFile)
 
         ####Read in the orbits
-        if self.orbitFile:
+        if '_POEORB_' in self.orbitFile:
             orb = self.extractPreciseOrbit()
-        else:
+        elif '_RESORB_' in self.orbitFile:
             orb = self.extractOrbit()
         
         self.product.orbit.setOrbitSource('Header')
         for sv in orb:
             self.product.orbit.addStateVector(sv)
 
-
         self.populateIPFVersion()
         self.extractBetaLUT()
         self.extractNoiseLUT()
@@ -423,10 +422,11 @@ class Sentinel1(Component):
 
         nsp = "{http://www.esa.int/safe/sentinel-1.0}"
 
-        if self.manifest.startswith('/vsizip'):
+        if '.zip' in self.manifest:
+
             import zipfile    
             parts = self.manifest.split(os.path.sep)
-            zipname = os.path.join(*(parts[2:-2]))
+            zipname = os.path.join('/',*(parts[:-2]))
             fname = os.path.join(*(parts[-2:]))
 
             try:
@@ -462,23 +462,25 @@ class Sentinel1(Component):
 
         datefmt = "%Y%m%dT%H%M%S"
         types = ['POEORB', 'RESORB']
+        filelist = []
         match = []
-        timeStamp = self.product.sensingMid
+        timeStamp = self.product.sensingStart+(self.product.sensingStop - self.product.sensingStart)/2.
         
         for orbType in types:
             files = glob.glob( os.path.join(self.orbitDir, 'S1A_OPER_AUX_' + orbType + '_OPOD*'))
-            
+            filelist.extend(files)
             ###List all orbit files
-            for result in files:
+
+        for result in filelist:
             fields = result.split('_')
             taft = datetime.datetime.strptime(fields[-1][0:15], datefmt)
             tbef = datetime.datetime.strptime(fields[-2][1:16], datefmt)
+            print(taft, tbef)
                 
             #####Get all files that span the acquisition
             if (tbef <= timeStamp) and (taft >= timeStamp):
                 tmid = tbef + 0.5 * (taft - tbef)
                 match.append((result, abs((timeStamp-tmid).total_seconds())))
-
             #####Return the file with the image is aligned best to the middle of the file
             if len(match) != 0:
                 bestmatch = min(match, key = lambda x: x[1])
@@ -516,13 +518,7 @@ class Sentinel1(Component):
             vec.setVelocity(vel)
             frameOrbit.addStateVector(vec)
 
-
+        return frameOrbit
-        orbExt = OrbitExtender(planet=Planet(pname='Earth'))
-        orbExt.configure()
-        newOrb = orbExt.extendOrbit(frameOrbit)
-
-
-        return newOrb
             
     def extractPreciseOrbit(self):
         '''
@@ -534,11 +530,10 @@ class Sentinel1(Component):
             print("IOError: %s" % strerr)
             return
 
-        _xml_root = ElementTree(file=fp).getroot()
+        _xml_root = ElementTree.ElementTree(file=fp).getroot()
 
         node = _xml_root.find('Data_Block/List_of_OSVs')
 
-        print('Extracting orbit from Orbit File: ', self.orbitFile)
         orb = Orbit()
         orb.configure()
 
@@ -582,10 +577,10 @@ class Sentinel1(Component):
         if self.calibrationXml is None:
             raise Exception('No calibration file provided')
 
-        if self.calibrationXml.startswith('/vsizip'):
+        if '.zip' in self.calibrationXml:
             import zipfile
             parts = self.calibrationXml.split(os.path.sep)
-            zipname = os.path.join(*(parts[2:-4]))
+            zipname = os.path.join('/',*(parts[:-4]))
             fname = os.path.join(*(parts[-4:]))
             
             try:
@@ -723,7 +718,7 @@ class Sentinel1(Component):
 
         print('Extracting normalized image ....')
 
-        src = gdal.Open(self.tiff.strip(), gdal.GA_ReadOnly)
+        src = gdal.Open('/vsizip//'+self.tiff.strip(), gdal.GA_ReadOnly)
         band = src.GetRasterBand(1)
 
         if self.product.numberOfSamples != src.RasterXSize:

contrib/stack/stripmapStack/crossmul.py

@@ -33,6 +33,8 @@ def cmdLineParse(iargs = None):
 
 def run(imageSlc1, imageSlc2, resampName, azLooks, rgLooks):
     objSlc1 = isceobj.createSlcImage()
+    #right now imageSlc1 and 2 are just text files, need to open them as image
+
     IU.copyAttributes(imageSlc1, objSlc1)
     objSlc1.setAccessMode('read')
     objSlc1.createImage()
@@ -81,7 +83,6 @@ def run(imageSlc1, imageSlc2, resampName, azLooks, rgLooks):
 
 
 def main(iargs=None):
-
     inps = cmdLineParse(iargs)
 
     img1 = isceobj.createImage()
@@ -97,8 +98,7 @@ def main(iargs=None):
 
 if __name__ == '__main__':
 
+    main()
     '''
     Main driver.
     '''
-
-