Merge pull request #87 from yunjunz/geometry

stripmapStack: generate multilooked geometry files

contrib/stack/stripmapStack/Stack.py

@@ -65,6 +65,8 @@ class config(object):
         self.f.write('master : ' + self.slcDir +'\n')
         self.f.write('dem : ' + self.dem +'\n')
         self.f.write('output : ' + self.geometryDir +'\n')
+        self.f.write('alks : ' + self.alks +'\n')
+        self.f.write('rlks : ' + self.rlks +'\n')
         if self.nativeDoppler:
             self.f.write('native : True\n')
         if self.useGPU:

contrib/stack/stripmapStack/topo.py

@@ -6,8 +6,11 @@ import numpy as np
 import shelve
 import os
 import datetime 
+import shutil
 from isceobj.Constants import SPEED_OF_LIGHT
 from isceobj.Util.Poly2D import Poly2D
+from mroipac.looks.Looks import Looks
+
 
 def createParser():
     '''
@@ -354,6 +357,86 @@ def runSimamp(outdir, hname='z.rdr'):
     simImage.renderHdr()
     hgtImage.finalizeImage()
     simImage.finalizeImage()
+    return
+
+
+def runMultilook(in_dir, out_dir, alks, rlks):
+    print('generate multilooked geometry files with alks={} and rlks={}'.format(alks, rlks))
+    from iscesys.Parsers.FileParserFactory import createFileParser
+    FP = createFileParser('xml')
+
+    if not os.path.isdir(out_dir):
+        os.makedirs(out_dir)
+        print('create directory: {}'.format(out_dir))
+
+    for fbase in ['hgt', 'incLocal', 'lat', 'lon', 'los', 'shadowMask', 'waterMask']:
+        fname = '{}.rdr'.format(fbase)
+        in_file = os.path.join(in_dir, fname)
+        out_file = os.path.join(out_dir, fname)
+
+        if os.path.isfile(in_file):
+            xmlProp = FP.parse(in_file+'.xml')[0]
+            if('image_type' in xmlProp and xmlProp['image_type'] == 'dem'):
+                inImage = isceobj.createDemImage()
+            else:
+                inImage = isceobj.createImage()
+
+            inImage.load(in_file+'.xml')
+            inImage.filename = in_file
+
+            lkObj = Looks()
+            lkObj.setDownLooks(alks)
+            lkObj.setAcrossLooks(rlks)
+            lkObj.setInputImage(inImage)
+            lkObj.setOutputFilename(out_file)
+            lkObj.looks()
+
+            # copy the full resolution xml/vrt file from ./merged/geom_master to ./geom_master
+            # to facilitate the number of looks extraction
+            # the file path inside .xml file is not, but should, updated
+            shutil.copy(in_file+'.xml', out_file+'.full.xml')
+            shutil.copy(in_file+'.vrt', out_file+'.full.vrt')
+
+    return out_dir
+
+
+def runMultilookGdal(in_dir, out_dir, alks, rlks):
+    print('generate multilooked geometry files with alks={} and rlks={}'.format(alks, rlks))
+    import gdal
+
+    # create 'geom_master' directory
+    if not os.path.isdir(out_dir):
+        os.makedirs(out_dir)
+        print('create directory: {}'.format(out_dir))
+
+    # multilook files one by one
+    for fbase in ['hgt', 'incLocal', 'lat', 'lon', 'los', 'shadowMask', 'waterMask']:
+        fname = '{}.rdr'.format(fbase)
+        in_file = os.path.join(in_dir, fname)
+        out_file = os.path.join(out_dir, fname)
+
+        if os.path.isfile(in_file):
+            ds = gdal.Open(in_file, gdal.GA_ReadOnly)
+            in_wid = ds.RasterXSize
+            in_len = ds.RasterYSize
+
+            out_wid = int(in_wid / rlks)
+            out_len = int(in_len / alks)
+            src_wid = out_wid * rlks
+            src_len = out_len * alks
+
+            cmd = 'gdal_translate -of ENVI -a_nodata 0 -outsize {ox} {oy} '.format(ox=out_wid, oy=out_len)
+            cmd += ' -srcwin 0 0 {sx} {sy} {fi} {fo} '.format(sx=src_wid, sy=src_len, fi=in_file, fo=out_file)
+            print(cmd)
+            os.system(cmd)
+
+            # copy the full resolution xml/vrt file from ./merged/geom_master to ./geom_master
+            # to facilitate the number of looks extraction
+            # the file path inside .xml file is not, but should, updated
+            shutil.copy(in_file+'.xml', out_file+'.full.xml')
+            shutil.copy(in_file+'.vrt', out_file+'.full.vrt')
+
+    return out_dir
 
 
 def extractInfo(frame, inps):
@@ -369,8 +452,8 @@ def extractInfo(frame, inps):
 
     info.lookSide = frame.instrument.platform.pointingDirection
     info.rangeFirstSample = frame.startingRange
-    info.numberRangeLooks = inps.rlks
-    info.numberAzimuthLooks = inps.alks
+    info.numberRangeLooks = 1 #inps.rlks
+    info.numberAzimuthLooks = 1 #inps.alks
 
     fsamp = frame.rangeSamplingRate
 
@@ -443,6 +526,14 @@ def main(iargs=None):
     runTopo(info,demImage,dop=doppler,nativedop=inps.nativedop, legendre=inps.legendre)
     runSimamp(os.path.dirname(info.heightFilename),os.path.basename(info.heightFilename))
 
+    # write multilooked geometry files in "geom_master" directory, same level as "Igrams"
+    if inps.rlks * inps.rlks > 1:
+        out_dir = os.path.join(os.path.dirname(os.path.dirname(info.outdir)), 'geom_master')
+        runMultilookGdal(in_dir=info.outdir, out_dir=out_dir, alks=inps.alks, rlks=inps.rlks)
+        #runMultilook(in_dir=info.outdir, out_dir=out_dir, alks=inps.alks, rlks=inps.rlks)
+
+    return
+
 
 if __name__ == '__main__':
     '''