Merge pull request #91 from yunjunz/water_mask

stripmapStack: generate water mask in radar coord

contrib/stack/stripmapStack/Stack.py

@@ -75,6 +75,17 @@ class config(object):
             self.f.write('useGPU : False\n')
         self.f.write('##########################'+'\n')
 
+    def createWaterMask(self, function):
+
+        self.f.write('##########################'+'\n')
+        self.f.write(function+'\n')
+        self.f.write('createWaterMask : '+'\n')
+        self.f.write('dem_file : ' + self.dem +'\n')
+        self.f.write('lat_file : ' + self.latFile +'\n')
+        self.f.write('lon_file : ' + self.lonFile +'\n')
+        self.f.write('output : ' + self.waterMaskFile + '\n')
+        self.f.write('##########################'+'\n')
+
     def geo2rdr(self, function):
 
         self.f.write('##########################'+'\n')
@@ -311,8 +322,7 @@ class run(object):
              self.runf.write(self.text_cmd+'stripmapWrapper.py -c '+ configName+'\n')
  
     def master_focus_split_geometry(self, stackMaster, config_prefix, split=False, focus=True, native=True):
-  	########
+        """focusing master and producing geometry files"""
-  	# focusing master and producing geometry files
         configName = os.path.join(self.configDir, config_prefix + stackMaster)
         configObj = config(configName)
         configObj.configure(self)
@@ -333,6 +343,14 @@ class run(object):
             configObj.outDir = configObj.slcDir
             configObj.shelve = os.path.join(configObj.slcDir, 'data')
             configObj.splitRangeSpectrum('[Function-{0}]'.format(counter))
+            counter += 1
+
+        # generate water mask in radar coordinates
+        configObj.latFile = os.path.join(self.workDir, 'geom_master/lat.rdr')
+        configObj.lonFile = os.path.join(self.workDir, 'geom_master/lon.rdr')
+        configObj.waterMaskFile = os.path.join(self.workDir, 'geom_master/waterMask.rdr')
+        configObj.createWaterMask('[Function-{0}]'.format(counter))
+        counter += 1
 
         configObj.finalize()
         del configObj

contrib/stack/stripmapStack/createWaterMask.py

@@ -2,30 +2,45 @@
 
 #Author: Heresh Fattahi
 
-import isce
+import os
-import isceobj
-from contrib.demUtils.SWBDStitcher import SWBDStitcher
-from iscesys.DataManager import createManager
 import argparse
 import configparser
-from numpy import round
+import numpy as np
+import isce
+import isceobj
+from iscesys.DataManager import createManager
+from contrib.demUtils.SWBDStitcher import SWBDStitcher
+
+
+EXAMPLE = """example:
+  createWaterMask.py -b 31 33 130 132
+  createWaterMask.py -b 31 33 130 132 -l lat.rdr -L lon.rdr -o waterMask.rdr
+  createWaterMask.py -d ../DEM/demLat_N31_N33_Lon_E130_E132.dem.wgs84 -l lat.rdr -L lon.rdr -o waterMask.rdr
+"""
 
 def createParser():
     '''
     Create command line parser.
     '''
     
-    parser = argparse.ArgumentParser( description='extracts the overlap geometry between master bursts')
+    parser = argparse.ArgumentParser(description='Create water body mask in geo and/or radar coordinates',
-   # parser.add_argument('-b', '--bbox', dest='bbox', type=str, default=None, 
+                                     formatter_class=argparse.RawTextHelpFormatter,
-    #                  help='Lat/Lon Bounding SNWE')
+                                     epilog=EXAMPLE)
-    parser.add_argument('-b', '--bbox', type = int, default = None, nargs = '+', dest = 'bbox', help = 'Defines the spatial region in the format south north west east.\
+    parser.add_argument('-b', '--bbox', dest='bbox', type=int, default=None, nargs=4, metavar=('S','N','W','E'), 
-                        The values should be integers from (-90,90) for latitudes and (0,360) or (-180,180) for longitudes.')
+                      help = 'Defines the spatial region in the format south north west east.\n'
+                             'The values should be integers from (-90,90) for latitudes '
+                             'and (0,360) or (-180,180) for longitudes.')
+    parser.add_argument('-d','--dem_file', dest='demName', type=str, default=None,
+                      help='DEM file in geo coordinates, i.e. demLat*.dem.wgs84.')
     parser.add_argument('-l', '--lat_file', dest='latName', type=str, default=None, 
                       help='pixel by pixel lat file in radar coordinate')
     parser.add_argument('-L', '--lon_file', dest='lonName', type=str, default=None, 
                       help='pixel by pixel lat file in radar coordinate')
+    parser.add_argument('-o', '--output', dest='outfile', type=str,
+                      help='output filename of water mask in radar coordinates')
     return parser
 
+
 def cmdLineParse(iargs = None):
     '''
     Command line parser.
@@ -33,37 +48,69 @@ def cmdLineParse(iargs = None):
 
     parser = createParser()
     inps = parser.parse_args(args=iargs)
-    #inps.bbox = [int(round(val)) for val in inps.bbox.split()]
+
+    if not inps.bbox and not inps.demName:
+        parser.print_usage()
+        raise SystemExit('ERROR: no --bbox/--dem_file input, at least one is required.')
+
+    if not inps.outfile and (inps.latName and inps.lonName):
+        inps.outfile = os.path.join(os.path.dirname(inps.latName), 'waterMask.rdr')
+
     return inps
 
 
-def download_waterMask(inps):
+def dem2bbox(dem_file):
+    """Grab bbox from DEM file in geo coordinates"""
+    demImage = isceobj.createDemImage()
+    demImage.load(dem_file + '.xml')
+    demImage.setAccessMode('read')
+    N = demImage.getFirstLatitude()
+    W = demImage.getFirstLongitude()
+    S = N + demImage.getDeltaLatitude() * demImage.getLength()
+    E = W + demImage.getDeltaLongitude() * demImage.getWidth()
+    bbox = [np.floor(S).astype(int), np.ceil(N).astype(int),
+            np.floor(W).astype(int), np.ceil(E).astype(int)]
+    return bbox
+
+
+def download_waterMask(bbox, dem_file):
+    out_dir = os.getcwd()
+    # update out_dir and/or bbox if dem_file is input
+    if dem_file:
+        out_dir = os.path.dirname(dem_file)
+        if not bbox:
+            bbox = dem2bbox(dem_file)
 
     sw = createManager('wbd')
     sw.configure()
-    inps.waterBodyGeo = sw.defaultName(inps.bbox)
+    #inps.waterBodyGeo = sw.defaultName(inps.bbox)
+    sw.outputFile = os.path.join(out_dir, sw.defaultName(bbox))
     sw._noFilling = False
-    #sw._fillingValue = -1.0
+    sw._fillingValue = -1.0    #fill pixels without DEM data with value of -1, same as water body
-    sw._fillingValue = 0.0
+    #sw._fillingValue = 0.0
-    sw.stitch(inps.bbox[0:2],inps.bbox[2:])
+    sw.stitch(bbox[0:2], bbox[2:])
+    return sw.outputFile
 
-    return inps
 
-def geo2radar(inps):
+def geo2radar(geo_file, rdr_file, lat_file, lon_file):
-    inps.waterBodyRadar = inps.waterBodyGeo + '.rdr'
+    #inps.waterBodyRadar = inps.waterBodyGeo + '.rdr'
     sw = SWBDStitcher()
-    sw.toRadar(inps.waterBodyGeo, inps.latName, inps.lonName, inps.waterBodyRadar)
+    sw.toRadar(geo_file, lat_file, lon_file, rdr_file)
+    return rdr_file
 
 #looks.py -i watermask.msk -r 4 -a 14 -o 'waterMask.14alks_4rlks.msk'
 
 #imageMath.py -e='a*b' --a=filt_20100911_20101027.int --b=watermask.14alks_4rlks.msk -o filt_20100911_20101027_masked.int -t cfloat -s BIL
 
+
 def main(iargs=None):
 
     inps = cmdLineParse(iargs)
-  inps = download_waterMask(inps)    
+    geo_file = download_waterMask(inps.bbox, inps.demName)    
     if inps.latName and inps.lonName:
-     inps = geo2radar(inps)
+        geo2radar(geo_file, inps.outfile, inps.latName, inps.lonName)
+    return
+
 
 if __name__ == '__main__' :
   '''