add dem_error.py --dem-err-file and save_qgis.py --zero-first option (#1007)

* dem_error: add `--dem-error-file` option to specify a customized file name for the estimated DEM error.

* save_qgis: add `--zf / --zero-first` option, same as `tsview.py`, to set displacement at 1st acquisition to zero.

Author: yunjunz

docs/QGIS.md

@@ -3,7 +3,7 @@ Displacement time-series can be exported as QGIS PS Time Series Viewer plugin co
 1. Download and install [QGIS](https://qgis.org/en/site/) if you have not done so.
 2. Download [ps-time-series](https://plugins.qgis.org/plugins/pstimeseries/) plugin as a *.zip file; and install it to QGIS through "Plugins -> Manage and install Plugins -> install from ZIP".
 3. On the left side of QGIS window, browse to the *.shp file generated by `save_qgis.py` and add to the layer.
-4. Open "View -> Layer Styling", then "Fill color -> Edit" and add the command below to color code the map based on velocity value.
+4. Open "View -> Panels -> Layer Styling", then on the right choose "Fill color -> Edit" and add the command below to color code the map based on velocity value.
 
 ```
 ramp_color('RdBu', scale_linear(VEL, -20, 20, 0, 1))
@@ -13,7 +13,7 @@ ramp_color('RdBu', scale_linear(VEL, -20, 20, 0, 1))
   <img width="1000" src="https://yunjunzhang.files.wordpress.com/2019/11/ps-time-series-viewer-1.png">
 </p>
 
-5. Select the [PS Time Series Viewer logo](https://gitlab.com/faunalia/ps-speed/blob/master/icons/logo.png) to activate the tool, and click/play on the map to display the time-series!
+5. Click the [PS Time Series Viewer logo](https://gitlab.com/faunalia/ps-speed/blob/master/icons/logo.png) to activate the tool, and click/play on the map to display the time-series!
 
 <p align="left">
   <img width="800" src="https://yunjunzhang.files.wordpress.com/2019/11/ps-time-series-viewer-2.png">

src/mintpy/cli/dem_error.py

@@ -41,15 +41,16 @@ def create_parser(subparsers=None):
     parser = arg_utils.create_argument_parser(
         name, synopsis=synopsis, description=synopsis, epilog=epilog, subparsers=subparsers)
 
-    parser.add_argument('timeseries_file',
-                        help='Timeseries file to be corrrected')
+    parser.add_argument('ts_file', help='Time-series HDF5 file to be corrrected.')
     parser.add_argument('-g', '--geometry', dest='geom_file',
                         help='geometry file including datasets:\n'+
                              'incidence angle\n'+
                              'slant range distance\n' +
                              'and/or 3D perpendicular baseline')
-    parser.add_argument('-o', '--outfile',
-                        help='Output file name for corrected time-series')
+    parser.add_argument('-o', '--outfile', dest='ts_cor_file',
+                        help='Output file name for corrected time-series (default: add suffix of "_demErr")')
+    parser.add_argument('--dem-err-file','--dem-error-file', dest='dem_err_file', default='demErr.h5',
+                        help='Output file name for the estimated DEM error (default: %(default)s).')
 
     defo_model = parser.add_argument_group('temporal deformation model')
     defo_model.add_argument('-t', '--template', dest='template_file',
@@ -69,7 +70,7 @@ def create_parser(subparsers=None):
                         help='Use phase velocity instead of phase for inversion constrain.')
     parser.add_argument('--update', dest='update_mode', action='store_true',
                         help='Enable update mode, and skip inversion if:\n'+
-                             '1) output timeseries file already exists, readable '+
+                             '1) output time-series file already exists, readable '+
                              'and newer than input interferograms file\n' +
                              '2) all configuration parameters are the same.')
     # computing
@@ -106,10 +107,18 @@ def cmd_line_parse(iargs=None):
     if inps.polyOrder < 1:
         raise argparse.ArgumentTypeError("Minimum polynomial order is 1")
 
-    # default: --output
-    if not inps.outfile:
-        fbase = os.path.splitext(inps.timeseries_file)[0]
-        inps.outfile = f'{fbase}_demErr.h5'
+    # check: --output / --dem-error-file option (must be HDF5 file)
+    for fname in [inps.ts_cor_file, inps.dem_err_file]:
+        if fname:
+            fext = os.path.splitext(fname)[1]
+            if fext not in ['.h5','.he5']:
+                msg = f'--output / --dem-err-file option ({fext}) supports HDF5 file only!'
+                raise ValueError(msg)
+
+    # default: --output option
+    if not inps.ts_cor_file:
+        fbase = os.path.splitext(inps.ts_file)[0]
+        inps.ts_cor_file = f'{fbase}_demErr.h5'
 
     return inps
 

src/mintpy/cli/save_qgis.py

@@ -37,6 +37,10 @@ def create_parser(subparsers=None):
     parser.add_argument('-B', '--geo-bbox', dest='geo_bbox', type=float, nargs=4, default=None,
                         metavar=('S','N','W','E'), help='bounding box in lat lon: South North West East')
 
+    # other options
+    parser.add_argument('--zf', '--zero-first', dest='zero_first', action='store_true',
+                        help='Set displacement at the first acquisition to zero.')
+
     return parser
 
 

src/mintpy/dem_error.py

@@ -35,28 +35,28 @@ def run_or_skip(inps):
     flag = 'skip'
 
     # check output file
-    if not os.path.isfile(inps.outfile):
+    if not os.path.isfile(inps.ts_cor_file):
         flag = 'run'
-        print(f'1) output file {inps.outfile} NOT found.')
+        print(f'1) output file {inps.ts_cor_file} NOT found.')
     else:
         # check if time-series file is partly written using file size
         # since time-series file is not compressed
-        with h5py.File(inps.outfile, 'r') as f:
+        with h5py.File(inps.ts_cor_file, 'r') as f:
             fsize_ref = f['timeseries'].size * 4
-        fsize = os.path.getsize(inps.outfile)
+        fsize = os.path.getsize(inps.ts_cor_file)
         if fsize <= fsize_ref:
             flag = 'run'
-            print(f'1) output file {inps.outfile} is NOT fully written.')
+            print(f'1) output file {inps.ts_cor_file} is NOT fully written.')
 
         else:
-            print(f'1) output file {inps.outfile} already exists.')
+            print(f'1) output file {inps.ts_cor_file} already exists.')
 
             # check modification time
-            infiles = [inps.timeseries_file]
+            infiles = [inps.ts_file]
             if inps.geom_file:
                 infiles.append(inps.geom_file)
             ti = max(os.path.getmtime(i) for i in infiles)
-            to = os.path.getmtime(inps.outfile)
+            to = os.path.getmtime(inps.ts_cor_file)
             if ti > to:
                 flag = 'run'
                 print(f'2) output file is NOT newer than input file: {infiles}.')
@@ -65,9 +65,9 @@ def run_or_skip(inps):
 
     # check configuration
     if flag == 'skip':
-        date_list_all = timeseries(inps.timeseries_file).get_date_list()
+        date_list_all = timeseries(inps.ts_file).get_date_list()
         inps.excludeDate = read_exclude_date(inps.excludeDate, date_list_all, print_msg=False)[1]
-        meta = readfile.read_attribute(inps.outfile)
+        meta = readfile.read_attribute(inps.ts_cor_file)
         if any(str(vars(inps)[key]) != meta.get(key_prefix+key, 'None') for key in config_keys):
             flag = 'run'
             print(f'3) NOT all key configuration parameters are the same:{config_keys}')
@@ -161,7 +161,7 @@ def get_design_matrix4defo(inps):
     model['periodic'] = inps.periodic
 
     # prepare SAR info
-    ts_obj = timeseries(inps.timeseries_file)
+    ts_obj = timeseries(inps.ts_file)
     date_list = ts_obj.get_date_list()
     seconds = ts_obj.get_metadata().get('CENTER_LINE_UTC', 0)
 
@@ -439,7 +439,7 @@ def correct_dem_error(inps):
     ## 1. input info
 
     # 1.1 read date info
-    ts_obj = timeseries(inps.timeseries_file)
+    ts_obj = timeseries(inps.ts_file)
     ts_obj.open()
     num_date = ts_obj.numDate
     length, width = ts_obj.length, ts_obj.width
@@ -465,20 +465,18 @@ def correct_dem_error(inps):
         meta[key_prefix+key] = str(vars(inps)[key])
 
     # 2.2 instantiate est. DEM error
-    dem_err_file = 'demErr.h5'
     meta['FILE_TYPE'] = 'dem'
     meta['UNIT'] = 'm'
     ds_name_dict = {'dem' : [np.float32, (length, width), None]}
-    writefile.layout_hdf5(dem_err_file, ds_name_dict, metadata=meta)
+    writefile.layout_hdf5(inps.dem_err_file, ds_name_dict, metadata=meta)
 
     # 2.3 instantiate corrected time-series
-    ts_cor_file = inps.outfile
     meta['FILE_TYPE'] = 'timeseries'
-    writefile.layout_hdf5(ts_cor_file, metadata=meta, ref_file=inps.timeseries_file)
+    writefile.layout_hdf5(inps.ts_cor_file, metadata=meta, ref_file=inps.ts_file)
 
     # 2.4 instantiate residual phase time-series
-    ts_res_file = os.path.join(os.path.dirname(inps.outfile), 'timeseriesResidual.h5')
-    writefile.layout_hdf5(ts_res_file, metadata=meta, ref_file=inps.timeseries_file)
+    ts_res_file = os.path.join(os.path.dirname(inps.ts_cor_file), 'timeseriesResidual.h5')
+    writefile.layout_hdf5(ts_res_file, metadata=meta, ref_file=inps.ts_file)
 
 
     ## 3. run the estimation and write to disk
@@ -503,7 +501,7 @@ def correct_dem_error(inps):
     # 3.2 prepare the input arguments for *_patch()
     data_kwargs = {
         'G_defo'         : G_defo,
-        'ts_file'        : inps.timeseries_file,
+        'ts_file'        : inps.ts_file,
         'geom_file'      : inps.geom_file,
         'date_flag'      : date_flag,
         'phase_velocity' : inps.phaseVelocity,
@@ -558,15 +556,15 @@ def correct_dem_error(inps):
         # DEM error - 2D
         block = [box[1], box[3], box[0], box[2]]
         writefile.write_hdf5_block(
-            dem_err_file,
+            inps.dem_err_file,
             data=delta_z,
             datasetName='dem',
             block=block)
 
         # corrected time-series - 3D
         block = [0, num_date, box[1], box[3], box[0], box[2]]
         writefile.write_hdf5_block(
-            ts_cor_file,
+            inps.ts_cor_file,
             data=ts_cor,
             datasetName='timeseries',
             block=block)
@@ -586,4 +584,4 @@ def correct_dem_error(inps):
     m, s = divmod(time.time()-start_time, 60)
     print(f'time used: {m:02.0f} mins {s:02.1f} secs.')
 
-    return dem_err_file, ts_cor_file, ts_res_file
+    return inps.dem_err_file, inps.ts_cor_file, ts_res_file

src/mintpy/save_qgis.py

@@ -49,20 +49,21 @@ def gather_files(ts_file, geom_file):
 
     #Can also add DEM Error here: corrected_DEM = DEM + DEM_error
     ts_dir = os.path.dirname(ts_file)
-    fDict = { 'TimeSeries' : ts_file,
-              'Velocity'   : os.path.join(ts_dir, vel_file),
-              'Coherence'  : os.path.join(ts_dir, coh_file),
-              'Mask'       : os.path.join(ts_dir, msk_file),
-              'Geometry'   : geom_file,
-            }
+    fDict = {
+        'TimeSeries' : ts_file,
+        'Velocity'   : os.path.join(ts_dir, vel_file),
+        'Coherence'  : os.path.join(ts_dir, coh_file),
+        'Mask'       : os.path.join(ts_dir, msk_file),
+        'Geometry'   : geom_file,
+    }
 
     #Check if the files exists.
     for fname in fDict.values():
         if not os.path.isfile(fname):
             raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), fname)
 
     for key, value in fDict.items():
-        print(f'{key:<12}: {value}')
+        print(f'{key:<10}: {value}')
     return fDict
 
 
@@ -86,13 +87,14 @@ def read_bounding_box(pix_box, geo_box, geom_file):
     return pix_box
 
 
-def write_shape_file(fDict, shp_file, box=None):
+def write_shape_file(fDict, shp_file, box=None, zero_first=False):
     '''Write time-series data to a shape file
 
-    Parameters: fDict    - dict, with value for path of data files
-                shp_file - str, output filename
-                box      - tuple of 4 int, in (x0, y0, x1, y1)
-    Returns:    shp_file - str, output filename
+    Parameters: fDict      - dict, with value for path of data files
+                shp_file   - str, output filename
+                box        - tuple of 4 int, in (x0, y0, x1, y1)
+                zero_first - bool, set displacement at 1st acquisition to zero
+    Returns:    shp_file   - str, output filename
     '''
 
     shpDriver = ogr.GetDriverByName("ESRI Shapefile")
@@ -161,6 +163,8 @@ def write_shape_file(fDict, shp_file, box=None):
     mask = readfile.read(fDict['Mask'], box=box)[0]
     nValid = np.sum(mask != 0)
     print('number of points with time-series:', nValid)
+    if zero_first:
+        print('set displacement at the first acquisition to zero.')
 
     lats, lons = ut.get_lat_lon(ts_obj.metadata, geom_file=fDict['Geometry'], box=box)
 
@@ -183,6 +187,9 @@ def write_shape_file(fDict, shp_file, box=None):
 
                         # read data for the line
                         ts = tsid['timeseries'][:, line, box[0]:box[2]].astype(np.float64)
+                        if zero_first:
+                            ts -= np.tile(ts[0, :], (ts.shape[0], 1))
+
                         coh = cohid['temporalCoherence'][line, box[0]:box[2]].astype(np.float64)
                         vel = velid['velocity'][line, box[0]:box[2]].astype(np.float64)
                         vel_std = velid['velocityStd'][line, box[0]:box[2]].astype(np.float64)
@@ -214,7 +221,7 @@ def write_shape_file(fDict, shp_file, box=None):
 
                             # update counter / progress bar
                             counter += 1
-                            prog_bar.update(counter, every=100, suffix=f'{counter}/{nValid}')
+                            prog_bar.update(counter, every=100, suffix=f'line {counter}/{nValid}')
                     prog_bar.close()
 
     print(f'finished writing to file: {shp_file}')
@@ -235,6 +242,6 @@ def save_qgis(inps):
     fDict = gather_files(inps.ts_file, inps.geom_file)
 
     # Write shape file
-    write_shape_file(fDict, inps.shp_file, box=box)
+    write_shape_file(fDict, inps.shp_file, box=box, zero_first=inps.zero_first)
 
     return