integratePS/unwrapStack: add --unw_method option

docs/workflows.md

@@ -50,18 +50,18 @@ To compute the amplitude dispersion of different pixels:
 ampdispersion.py -i coreg_stack/slcs_base.vrt -o ampDispersion/ampdispersion -m ampDispersion/mean
 ```
 
-To choose the PS pixels one can simply threshold the amplitude dispersion using a custom script. Here is an example using ISCE's imageMath funstion to threshold the PS pixels with a threshold of 0.4: 
+To choose the PS pixels one can simply threshold the amplitude dispersion using a custom script. Here is an example using ISCE's imageMath function to threshold the PS pixels with a threshold of 0.4: 
 
 ```
-imageMath.py -e="a<0.4" --a=ampDispersion/ampdispersion  -o ps_pixels -t byte
+imageMath.py -e="a<0.4" --a=ampDispersion/ampdispersion  -o ampDispersion/ps_pixels -t byte
 ```
 
 ### Integrate PS pixels to the wrapped phase series
 
 For PS pixels we extract the wrapped phase series of the SLCs through time and integrate to the maps of wrapped phase time-series from the previous step. As a result of this step the wrapped phase time-series of PS and DS pixels are obtained in which the phase series of DS pixels were estimated from full covariance analysis over local neighborhoods as explained above and for PS pixels the wrapped phase series was extracted from SLCs.
 
 ```
-integratePS.py -s coreg_stack/slcs_base.vrt -d adjusted_wrapped_DS/ -t Sequential/Datum_connection/EVD/tcorr.bin -p ampDispersion/ps_pixels -o PS_DS -c cor
+integratePS.py -s coreg_stack/slcs_base.vrt -d adjusted_wrapped_DS/ -t Sequential/Datum_connection/EVD/tcorr.bin -p ampDispersion/ps_pixels -o PS_DS --unwrap_method snaphu
 ```
 
 ## Phase 2: Wrapped time-series to deformation time-series
@@ -71,13 +71,15 @@ FRInGE currently is focused on estimating wrapped phase time-series for PS and D
 
 ### Unwrap the wrapped phase series
 
-To prepare unwrapping commands for each epoch of time-series run unwrapStacck.py:
+To prepare unwrapping commands for each epoch of time-series run `unwrapStack.py`. This will write the unwrapping command to a shell script _run_unwrap.sh_.
 
 ```
-unwrapStacck.py -s slcs -m Sequential/miniStacks/ -d Sequential/Datum_connection/  -M 15 -o unwrapped
+unwrapStack.py -s slcs -m Sequential/miniStacks/ -d Sequential/Datum_connection/ -M 15 -u 'unwrap_fringe.py -m snaphu'
 ```
 
-By running the command above, the unwrapping commands will be written to shell script. Each row of the run file is then the unwrapping command for each epoch which is indpendent from other epochs and therefore unwrapping different epochs can be run in parallel if computing resourses such as clusters are available.  
+Alternatively one may use `integratePS.py --unw_method snaphu` in previous step, which will write a shell script _run_unwrap_ps_ds.sh_ (recommended).
+
+Each row of the run file is then the unwrapping command for each epoch which is indpendent from other epochs and therefore unwrapping different epochs can be run in parallel if computing resourses such as clusters are available.  
 
 TODO: 
 ingest FRInGE unwrapped phase series to MintPy ...

python/adjustMiniStacks.py

@@ -188,7 +188,7 @@ def adjust_acquisition_wrapped_vrt(miniStackSlc, adjustSlc, output, length, widt
       <SourceFilename>{miniStackSlc}</SourceFilename>
     </SimpleSource>
     <SimpleSource>
-      <SourceFilename >{adjustment}</SourceFilename>
+      <SourceFilename>{adjustment}</SourceFilename>
     </SimpleSource>
   </VRTRasterBand>
 </VRTDataset>'''

python/integratePS.py

@@ -23,6 +23,7 @@ def cmdLineParser(iargs = None):
 
     parser = argparse.ArgumentParser(description = 'integrate PS pixels into existing unwrapped DS filed',
                 formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+
     parser.add_argument('-s', '--slc_stack', type=str, dest='slcStack',
             required=True, help='slc stack dataset ')
 
@@ -39,7 +40,10 @@ def cmdLineParser(iargs = None):
             required=True, help='The output directory to store the wrapped phase series of DS and PS pixels')
 
     parser.add_argument('-c', '--coreg_slc_dir', type=str, dest='coregSlcDir',
-            required=True, help='The directory that contains the coregistered stack of SLCs')
+            required=False, help='The directory that contains the coregistered stack of SLCs')
+
+    parser.add_argument('-u', '--unw_method', '--unwrap_method', type=str, dest='unwrapMethod', choices=('snaphu','phass'),
+            help='phase unwrapping method. e.g., snaphu, phass. If enabled, a shell file "run_unwrap_ps_ds.sh" with unwrap command will be created.')
 
     return parser.parse_args()
 
@@ -93,22 +97,22 @@ def integratePS2DS(dsDS_band_i, dsDS_band_j, dsSlc, dsTcor, psDataset, outDatase
 
     x0 = 0
     xoff = ncols
-    
+
     for block in range(nblocks):
 
-        print("block: ", block)        
+        print("block: ", block)
 
         y0 = block*linesPerBlock
 
         if (y0 + linesPerBlock) > nrows:
-           yoff = nrows - y0
+            yoff = nrows - y0
         else:
-           yoff = linesPerBlock
+            yoff = linesPerBlock
 
         tcorr = dsTcor.ReadAsArray(x0, y0, xoff, yoff)
 
         psPixels = psDataset.ReadAsArray(x0, y0, xoff, yoff)
-        
+
         band_i_DS = dsDS_band_i.ReadAsArray(x0, y0, xoff, yoff)
         band_j_DS = dsDS_band_j.ReadAsArray(x0, y0, xoff, yoff)
 
@@ -119,10 +123,10 @@ def integratePS2DS(dsDS_band_i, dsDS_band_j, dsSlc, dsTcor, psDataset, outDatase
 
         # get the data for PS pixels
         ifgram_ds_ps[psPixels == 1] = ifgram[psPixels == 1]
- 
-        print(ifgram.dtype) 
+
+        print(ifgram.dtype)
         outDataset.GetRasterBand(1).WriteArray(ifgram_ds_ps, x0, y0)
-        
+
     return None
 
 
@@ -139,16 +143,16 @@ def getCoherence(dsTcor, psDataset, outDataset,
         y0 = block*linesPerBlock
 
         if (y0 + linesPerBlock) > nrows:
-           yoff = nrows - y0
+            yoff = nrows - y0
         else:
-           yoff = linesPerBlock
+            yoff = linesPerBlock
 
         tcorr = dsTcor.ReadAsArray(x0, y0, xoff, yoff)
 
         psPixels = psDataset.ReadAsArray(x0, y0, xoff, yoff)
-        
+
         tcorr[psPixels==1] = 0.95
-        
+
         outDataset.GetRasterBand(1).WriteArray(tcorr, x0, y0)
 
     return None
@@ -157,21 +161,21 @@ def getCoherence(dsTcor, psDataset, outDataset,
 def main(iargs=None):
 
     inps = cmdLineParser(iargs)
-   
-    # Open the SLC dataset to read 
+
+    # Open the SLC dataset to read
     dsSlc = gdal.Open(inps.slcStack, gdal.GA_ReadOnly)
-    
+
     # Open the tcorr dataset
     dsTcor = gdal.Open(inps.tcorrFile, gdal.GA_ReadOnly)
-    
+
     # Open the PS pixels dataset
     psDataset = gdal.Open(inps.psPixelsFile, gdal.GA_ReadOnly)
 
     nSlc = dsSlc.RasterCount
     nrows = dsSlc.RasterYSize
     ncols = dsSlc.RasterXSize
 
-    print("number of SLC: ", nSlc) 
+    print("number of SLC: ", nSlc)
     print("number of rows: ", nrows)
     print("number of columns: ", ncols)
 
@@ -182,30 +186,23 @@ def main(iargs=None):
     elif (nrows % (nblocks * linesPerBlock) != 0):
         nblocks += 1
 
-    print("nblocks: ", nblocks)        
+    print("nblocks: ", nblocks)
 
     if not os.path.exists(inps.outDir):
         os.makedirs(inps.outDir)
 
-    unwDir = os.path.join(inps.outDir, "unwrap")
-
-    if not os.path.exists(unwDir):
-        os.makedirs(unwDir)
-
     driver = gdal.GetDriverByName("ENVI")
     corName = os.path.join(inps.outDir, "tcorr_ds_ps.bin")
     corDataset = driver.Create(corName, ncols, nrows, 1, gdal.GDT_Float32)
     getCoherence(dsTcor, psDataset, corDataset,
-              ncols, nrows, nblocks, linesPerBlock)
-    corDataset.FlushCache()   
-
+                 ncols, nrows, nblocks, linesPerBlock)
+    corDataset.FlushCache()
 
     # extract the list of the dates
     dateList = []
     for band in range(nSlc):
         date_i =dsSlc.GetRasterBand(band + 1).GetMetadata("slc")["Date"]
         dateList.append(date_i)
-
     dateList.sort()
 
     # setup a network, compute coherence based on geometry, find min span tree pairs
@@ -220,7 +217,7 @@ def main(iargs=None):
     rangeSpacing = 1.02*29
     coherenceMatrix = networkObj.geometrical_coherence(rangeSpacing, theta, wvl, r)
     networkObj.min_span_tree(coherenceMatrix, 1)
-    
+
     networkObj.plot_network()
     '''
     networkObj.single_master()
@@ -229,15 +226,11 @@ def main(iargs=None):
     for pair in networkObj.pairsDates:
         print(pair)
     print(len(networkObj.pairsDates))
-
-    # an output script with commands to unwrap interferograms 
-    run_outname = "run_unwrap_ps_ds.sh"
-    runf= open(run_outname,'w')
 
     # loop over all pairs
     for pair in networkObj.pairsDates:
-        
-        date_i = pair.split('-')[0] 
+
+        date_i = pair.split('-')[0]
         date_j = pair.split('-')[1]
 
         band_i = dateList.index(date_i) + 1
@@ -254,22 +247,33 @@ def main(iargs=None):
         dsDS_band_j = gdal.Open(os.path.join(inps.dsStackDir, date_j + ".slc.vrt"), gdal.GA_ReadOnly)
 
         # integrate PS to DS for this pair and write to file block by block
-
         integratePS2DS(dsDS_band_i, dsDS_band_j, dsSlc, dsTcor, psDataset, outDataset,
                        nblocks, linesPerBlock, band_i, band_j)
-        
+
         # close the dataset
         outDataset.FlushCache()
         dsDS_band_i = None
         dsDS_band_j = None
 
-        # write the unwrapping command for this pair
-        unwrapName = os.path.join(unwDir, "{0}_{1}.unw".format(date_i, date_j))
-        cmd = "unwrap_fringe.py -m snaphu -i " + outName + " -c " + corName + " -o " + unwrapName
-        runf.write(cmd + "\n")
-
-    runf.close()
-
+    # write the unwrapping command for this pair
+    if inps.unwrapMethod is not None:
+        # prepare output dir
+        unwDir = os.path.join(inps.outDir, "unwrap")
+        if not os.path.exists(unwDir):
+            os.makedirs(unwDir)
+
+        # an output script with commands to unwrap interferograms
+        run_outname = "run_unwrap_ps_ds.sh"
+        runf= open(run_outname,'w')
+
+        for pair in networkObj.pairsDates:
+            date_i = pair.split('-')[0]
+            date_j = pair.split('-')[1]
+            intName = os.path.join(inps.outDir, "{0}_{1}.int".format(date_i, date_j))
+            unwName = os.path.join(unwDir, "{0}_{1}.unw".format(date_i, date_j))
+            cmd = "unwrap_fringe.py -m " + inps.unwrapMethod + " -i " + intName + " -c " + corName + " -o " + unwName
+            runf.write(cmd + "\n")
+        runf.close()
 
     dsSLC = None
     dsUnw = None
@@ -281,5 +285,3 @@ def main(iargs=None):
     Main driver.
     '''
     main()
-
-

python/unwrapStack.py

@@ -32,11 +32,13 @@ def cmdLineParser():
             required=True, help='size of each miniStack')
 
     parser.add_argument('-o', '--outDir', type=str, dest='outDir',
-            required=True, help='output directory')
+            required=False, help='output directory')
 
     parser.add_argument('-u', '--unwrapper_script', type=str, dest='unwrapperScript',
-            default="unwrap_fringe.py", help='unwraper script. e.g., unwrap.py')
-
+            default="unwrap_fringe.py", help='unwraper script. e.g., unwrap_fringe.py')
+
+    parser.add_argument('--unw_method','--unwrap_method', type=str, dest='unwrapMethod',
+            default='phass', choices=('snaphu','phass'), help='unwrap method.')
 
     return parser.parse_args()
 
@@ -173,7 +175,7 @@ def write_xml(outFile, length, width):
         print("coherence: ", coh) 
         print("adjusted output phase: ", output)
         print("*****************")
-        cmd = inps.unwrapperScript + " -m phass -i " + miniStackSlc + " -c " + coh + " -o " + output
+        cmd = inps.unwrapperScript + " -m " + inps.unwrapMethod + " -i " + miniStackSlc + " -c " + coh + " -o " + output
         runf.write(cmd + "\n") 
 
     for k in datumDict.keys():
@@ -185,7 +187,7 @@ def write_xml(outFile, length, width):
         print("coherence: ", coh)
         print("adjusted output phase: ", output)
         print("*****************")
-        cmd = inps.unwrapperScript + " -m phass -i " + datumSlc + " -c " + coh + " -o " + output
+        cmd = inps.unwrapperScript + " -m " + inps.unwrapMethod + " -i " + datumSlc + " -c " + coh + " -o " + output
         runf.write(cmd + "\n")
 
     runf.close()

src/evd/docs/README.md

@@ -52,8 +52,7 @@ For full details on the processing approach, see references.
 
 1. Ansari, H., De Zan, F., Adam, N., Goel, K., & Bamler, R. (2016, July). Sequential estimator for distributed scatterer interferometry. In Geoscience and Remote Sensing Symposium (IGARSS), 2016 IEEE International (pp. 6859-6862). IEEE.
 
-2. Ansari, H., De Zan, F., & Bamler, R. (2017). Sequential Estimator- A Proposal for High-Precision and Efficient Earth Deformation Monitoring with InSAR [FRINGE 2017 Presentation](http://fringe.esa.int/files/presentation324.pdf) 
-
+2. Ansari, H., De Zan, F., & Bamler, R. (2017). Sequential Estimator- A Proposal for High-Precision and Efficient Earth Deformation Monitoring with InSAR [FRINGE 2017 Presentation](http://fringe2017.esa.int/files/presentation324.pdf) 
 
 3. A. Monti Guarnieri and S. Tebaldini, “On the Exploitation of Target Statistics for SAR Interferometry Applications,” IEEE Trans. Geosci. Remote Sens., vol. 46, no. 11, pp. 3436–3443, Nov. 2008.