Merge pull request #274 from GeoscienceAustralia/develop

Release 0.4.2

.travis.yml

@@ -53,14 +53,11 @@ install:
   - pip install GDAL==$(gdal-config --version)
   - python setup.py install
   - rm -rf Py_Rate.egg-info  # remove the local egg
-  - export PYRATEPATH=$(pwd)
   - export PYTHONPATH=$PYRATEPATH:$PYTHONPATH
   - chmod 444 tests/test_data/small_test/tif/geo_070709-070813_unw.tif  # makes the file readonly, used in a test
 
 # command to run tests, e.g. python setup.py test
 script:
-  #  - python scripts/update_placeholder_paths.py
-  - mpirun -n 3 pytest tests/test_mpi.py
   - pytest --cov-config=.coveragerc --cov-report term-missing:skip-covered --cov=pyrate tests/
 
 
@@ -79,7 +76,7 @@ deploy:
   verbose: true
   on:
     branch: master
-    condition: $GDALVERSION="3.0.2" && $TRAVIS_PYTHON_VERSION=3.8.*
+    python: 3.8
   github_token: $GITHUB_TOKEN2
   local_dir: docs/_build/html
   project_name: PyRate

docs/history.rst

@@ -2,8 +2,50 @@
 
 Release History
 ===============
+
+0.4.2 (2020-06-26)
+------------------
+Added
++++++
+- Save full-res coherence files to disk in ``conv2tif`` step if ``cohmask = 1``.
+- Save multi-looked coherence files to disk in ``prepifg`` step if ``cohmask = 1``.
+- Additional ``DATA_TYPE`` geotiff header metadata for above coherence files.
+- ``conv2tif`` and ``prepifg`` output files have a tag applied to filename dependent
+  on data type, i.e. ``_ifg.tif``, ``_coh.tif``, ``_dem.tif``.
+- Metadata about used reference pixel is added to interferogram geotiff headers:
+  lat/lon and x/y values; mean and standard deviation of reference window samples.
+- Quicklook PNG and KML files are generated for the ``Stack Rate`` error map by default.
+
+Changed
++++++++
+- Bugfix: ensure ``prepifg`` treats input data files as `read only`.
+- Bugfix: fix the way that the reference phase is subtracted from interferograms
+  during ``process`` step.
+- Bugfix: manual entry of ``refx/y`` converted to type ``int``.
+- User supplies latitude and longitude values when specifying a reference pixel in
+  the config file. Pixel x/y values are calculated and used internally.
+- Move ``Stack Rate`` masking to a standalone function ``pyrate.core.stack.mask_rate``,
+  applied during the ``merge`` step and add unit tests.
+- Skip ``Stack Rate`` masking if threshold parameter ``maxsig = 0``.
+- Provide log message indicating the percentage of pixels masked by 
+  ``pyrate.core.stack.mask_rate``.
+- Refactor ``pyrate.core.stack`` module; expose two functions in documentation:
+  i) single pixel stacking algorithm, and
+  ii) loop function for processing full ifg array.
+- Refactor ``pyrate.merge`` script; remove duplicated code and create reusable
+  generic functions.
+- Colourmap used to render quicklook PNG images is calculated from min/max values of
+  the geotiff band.
+- Updated ``test`` and ``dev`` requirements.
+
+Removed
++++++++
+- Deprecate unused functions in ``pyrate.core.config`` and corresponding tests.
+- Static colourmap ``utils/colourmap.txt`` that was previously used to render
+  quicklook PNG images is removed. 
+
 0.4.1 (2020-05-19)
------------------------
+------------------
 Added
 +++++
 - Python 3.8 support.
@@ -31,7 +73,7 @@ Removed
 - Deprecate ``parallel = 2`` option; splitting image via rows for parallelisation.
 
 0.4.0 (2019-10-31)
------------------------
+------------------
 Added
 +++++
 - Python 3.7 support.
@@ -64,7 +106,7 @@ Removed
 - Unused tests for legacy api.
 
 0.3.0 (2019-07-26)
------------------------
+------------------
 Added
 +++++
 - ``utils/apt_install.sh`` script that lists Ubuntu/apt package requirements.

input_parameters.conf

@@ -4,24 +4,27 @@
 # Optional ON/OFF switches - ON = 1; OFF = 0
 
 # Coherence masking (PREPIFG)
-cohmask:   1
+cohmask:    0
 
 # Orbital error correction (PROCESS)
-orbfit:        1
+orbfit:     1
 
 # APS correction using spatio-temporal filter (PROCESS)
-apsest:         0
+apsest:     0
 
 # Time series calculation (PROCESS)
-tscal:         1
+tscal:      1
+
+# Optional save of numpy array files for output products (MERGE)
+savenpy:    0
 
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # Multi-threading parameters used by stacking/timeseries/prepifg
 # gamma prepifg runs in parallel on a single machine if parallel = 1
 # parallel: 1 = parallel, 0 = serial
-parallel:  0
+parallel:   0
 # number of processes
-processes: 8
+processes:  8
 
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # Input/Output file locations
@@ -85,8 +88,8 @@ nan_conversion: 1
 # refnx/y: number of search grid points in x/y image dimensions
 # refchipsize: size of the data window at each search grid point
 # refminfrac: minimum fraction of valid (non-NaN) pixels in the data window
-refx:          
-refy:         
+refx:          150.941666654
+refy:          -34.218333314
 refnx:         5
 refny:         5
 refchipsize:   5
@@ -149,9 +152,9 @@ ts_pthr:       10
 #------------------------------------
 # Stacking calculation parameters
 
-# pthr: minimum number of coherent ifg connections for each pixel
-# nsig: n-sigma used as residuals threshold for iterative least squares stacking
-# maxsig: maximum residual used as a threshold for values in the rate map
+# pthr: threshold for minimum number of ifg observations for each pixel
+# nsig: threshold for iterative removal of observations
+# maxsig: maximum sigma (std dev) used as an output masking threshold applied in Merge step. 0 = OFF.
 pthr:          5
 nsig:          3
 maxsig:        1000

pyrate/configuration.py

@@ -13,13 +13,16 @@
 #   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
+"""
+This Python module contains utilities to validate user input parameters
+parsed in a PyRate configuration file.
+"""
 from configparser import ConfigParser
 from pathlib import Path, PurePath
-import re
 from pyrate.constants import NO_OF_PARALLEL_PROCESSES
 from pyrate.default_parameters import PYRATE_DEFAULT_CONFIGURATION
 from pyrate.core.algorithm import factorise_integer
-from pyrate.core.shared import extract_epochs_from_filename
+from pyrate.core.shared import extract_epochs_from_filename, InputTypes
 from pyrate.core.config import parse_namelist, ConfigException
 
 
@@ -44,17 +47,17 @@ def validate_parameter_value(input_name, input_value, min_value=None, max_value=
         if min_value is not None:
             if input_value < min_value:  # pragma: no cover
                 raise ValueError(
-                    "Invalid value for " + str(input_name) + " supplied: " + str(input_value) + ". Please provided a valid value greater than " + str(min_value) + ".")
+                    "Invalid value for " + str(input_name) + " supplied: " + str(input_value) + ". Provide a value greater than or equal to " + str(min_value) + ".")
     if input_value is not None:
         if max_value is not None:
             if input_value > max_value:  # pragma: no cover
                 raise ValueError(
-                    "Invalid value for " + str(input_name) + " supplied: " + str(input_value) + ". Please provided a valid value less than " + str(max_value) + ".")
+                    "Invalid value for " + str(input_name) + " supplied: " + str(input_value) + ". Provide a value less than or equal to " + str(max_value) + ".")
 
     if possible_values is not None:
         if input_value not in possible_values:  # pragma: no cover
             raise ValueError(
-                "Invalid value for " + str(input_name) + " supplied: " + str(input_value) + ". Please provided a valid value from with in: " + str(possible_values) + ".")
+                "Invalid value for " + str(input_name) + " supplied: " + str(input_value) + ". Provide a value from: " + str(possible_values) + ".")
     return True
 
 
@@ -74,7 +77,9 @@ def validate_file_list_values(file_list, no_of_epochs):
 
 
 class MultiplePaths:
-    def __init__(self, out_dir, file_name, ifglksx=1, ifgcropopt=1):
+    def __init__(self, out_dir: str, file_name: str, ifglksx: int = 1, ifgcropopt: int = 1,
+                 input_type: InputTypes = InputTypes.IFG):
+        self.input_type = input_type
         b = Path(file_name)
         if b.suffix == ".tif":
             self.unwrapped_path = None
@@ -83,7 +88,8 @@ def __init__(self, out_dir, file_name, ifglksx=1, ifgcropopt=1):
                 b.stem + '_' + str(ifglksx) + "rlks_" + str(ifgcropopt) + "cr.tif").as_posix()
         else:
             self.unwrapped_path = b.as_posix()
-            converted_path = Path(out_dir).joinpath(b.stem + '_' + b.suffix[1:]).with_suffix('.tif')
+            converted_path = Path(out_dir).joinpath(
+                b.stem.split('.')[0] + '_' + b.suffix[1:] + input_type.value).with_suffix('.tif')
             self.sampled_path = converted_path.with_name(
                 converted_path.stem + '_' + str(ifglksx) + "rlks_" + str(ifgcropopt) + "cr.tif").as_posix()
         self.converted_path = converted_path.as_posix()
@@ -171,20 +177,21 @@ def __init__(self, config_file_path):
         if self.cohfilelist is not None:
             # if self.processor != 0:  # not roipac
             validate_file_list_values(self.cohfilelist, 1)
-            self.coherence_file_paths = self.__get_files_from_attr('cohfilelist')
+            self.coherence_file_paths = self.__get_files_from_attr('cohfilelist', input_type=InputTypes.COH)
 
-        self.header_file_paths = self.__get_files_from_attr('hdrfilelist')
+        self.header_file_paths = self.__get_files_from_attr('hdrfilelist', input_type=InputTypes.HEADER)
 
         self.interferogram_files = self.__get_files_from_attr('ifgfilelist')
 
-        self.dem_file = MultiplePaths(self.outdir, self.demfile, self.ifglksx, self.ifgcropopt)
+        self.dem_file = MultiplePaths(self.outdir, self.demfile, self.ifglksx, self.ifgcropopt,
+                                      input_type=InputTypes.DEM)
 
         # backward compatibility for string paths
         for key in self.__dict__:
             if isinstance(self.__dict__[key], PurePath):
                 self.__dict__[key] = str(self.__dict__[key])
 
-    def __get_files_from_attr(self, attr):
+    def __get_files_from_attr(self, attr, input_type=InputTypes.IFG):
         val = self.__getattribute__(attr)
         files = parse_namelist(val)
-        return [MultiplePaths(self.outdir, p, self.ifglksx, self.ifgcropopt) for p in files]
+        return [MultiplePaths(self.outdir, p, self.ifglksx, self.ifgcropopt, input_type=input_type) for p in files]

pyrate/constants.py

@@ -27,7 +27,6 @@
 PROCESS = 'process'
 MERGE = 'merge'
 
-REF_COLOR_MAP_PATH = os.path.join(PYRATEPATH, "utils", "colourmap.txt")
 # distance division factor of 1000 converts to km and is needed to match legacy output
 DISTFACT = 1000
 # mappings for metadata in header for interferogram

pyrate/conv2tif.py

@@ -22,9 +22,11 @@
 from typing import Tuple, List
 from joblib import Parallel, delayed
 import numpy as np
+from pathlib import Path
 
 from pyrate.core.prepifg_helper import PreprocessError
 from pyrate.core import shared, mpiops, config as cf, gamma, roipac
+from pyrate.core import ifgconstants as ifc
 from pyrate.core.logger import pyratelogger as log
 from pyrate.configuration import MultiplePaths
 from pyrate.core.shared import mpi_vs_multiprocess_logging
@@ -103,7 +105,9 @@ def _geotiff_multiprocessing(unw_path: MultiplePaths, params: dict) -> Tuple[str
             header = roipac.roipac_header(unw_path.unwrapped_path, params)
         else:
             raise PreprocessError('Processor must be ROI_PAC (0) or GAMMA (1)')
+        header[ifc.INPUT_TYPE] = unw_path.input_type
         shared.write_fullres_geotiff(header, unw_path.unwrapped_path, dest, nodata=params[cf.NO_DATA_VALUE])
+        Path(dest).chmod(0o444)  # readonly output
         return dest, True
     else:
         log.warning(f"Full-res geotiff already exists in {dest}! Returning existing geotiff!")

pyrate/core/aps.py

@@ -19,32 +19,32 @@
 signals.
 """
 # pylint: disable=invalid-name, too-many-locals, too-many-arguments
-import logging
 import os
 from copy import deepcopy
 from collections import OrderedDict
 import numpy as np
 from numpy import isnan
 from scipy.fftpack import fft2, ifft2, fftshift, ifftshift
 from scipy.interpolate import griddata
+from pyrate.core.logger import pyratelogger as log
 
 from pyrate.core import shared, ifgconstants as ifc, mpiops, config as cf
 from pyrate.core.covariance import cvd_from_phase, RDist
 from pyrate.core.algorithm import get_epochs
 from pyrate.core.shared import Ifg
 from pyrate.core.timeseries import time_series
-from pyrate.merge import _assemble_tiles
-
-log = logging.getLogger(__name__)
+from pyrate.merge import assemble_tiles
 
 
 def wrap_spatio_temporal_filter(ifg_paths, params, tiles, preread_ifgs):
     """
     A wrapper for the spatio-temporal filter so it can be tested.
     See docstring for spatio_temporal_filter.
     """
-    if not params[cf.APSEST]:
-        log.info('APS correction not required.')
+    if params[cf.APSEST]:
+        log.info('Doing APS spatio-temporal filtering')
+    else:
+        log.info('APS spatio-temporal filtering not required')
         return
 
     # perform some checks on existing ifgs
@@ -101,16 +101,15 @@ def _calc_svd_time_series(ifg_paths, params, preread_ifgs, tiles):
     new_params[cf.TIME_SERIES_METHOD] = 2  # use SVD method
 
     process_tiles = mpiops.array_split(tiles)
-    output_dir = params[cf.TMPDIR]
 
     nvels = None
     for t in process_tiles:
         log.debug('Calculating time series for tile {} during APS '
                  'correction'.format(t.index))
         ifg_parts = [shared.IfgPart(p, t, preread_ifgs, params) for p in ifg_paths]
-        mst_tile = np.load(os.path.join(output_dir, 'mst_mat_{}.npy'.format(t.index)))
+        mst_tile = np.load(os.path.join(params[cf.TMPDIR], 'mst_mat_{}.npy'.format(t.index)))
         tsincr = time_series(ifg_parts, new_params, vcmt=None, mst=mst_tile)[0]
-        np.save(file=os.path.join(output_dir, 'tsincr_aps_{}.npy'.format(t.index)), arr=tsincr)
+        np.save(file=os.path.join(params[cf.TMPDIR], 'tsincr_aps_{}.npy'.format(t.index)), arr=tsincr)
         nvels = tsincr.shape[2]
 
     nvels = mpiops.comm.bcast(nvels, root=0)
@@ -125,11 +124,14 @@ def _assemble_tsincr(ifg_paths, params, preread_ifgs, tiles, nvels):
     """
     Helper function to reconstruct time series images from tiles
     """
+    # pre-allocate dest 3D array
     shape = preread_ifgs[ifg_paths[0]].shape + (nvels,)
     tsincr_g = np.empty(shape=shape, dtype=np.float32)
+    # shape of one 2D time-slice array
+    s = preread_ifgs[ifg_paths[0]].shape
+    # loop over the time slices and assemble dest 3D array
     for i in range(nvels):
-        for n, t in enumerate(tiles):
-            _assemble_tiles(i, n, t, tsincr_g[:, :, i], params[cf.TMPDIR], 'tsincr_aps')
+        tsincr_g[:, :, i] = assemble_tiles(s, params[cf.TMPDIR], tiles, out_type='tsincr_aps', index=i)
 
     return tsincr_g
 
@@ -184,7 +186,7 @@ def spatial_low_pass_filter(ts_lp, ifg, params):
     :return: ts_hp: filtered time series data of shape (ifg.shape, n_epochs)
     :rtype: ndarray
     """
-    log.info('Applying APS spatial low-pass filter')
+    log.info('Applying spatial low-pass filter')
     if params[cf.SLPF_NANFILL] == 0:
         ts_lp[np.isnan(ts_lp)] = 0  # need it here for cvd and fft
     else:
@@ -281,7 +283,7 @@ def temporal_low_pass_filter(tsincr, epochlist, params):
     :return: tsfilt_incr: filtered time series data, shape (ifg.shape, nepochs)
     :rtype: ndarray
     """
-    log.info('Applying APS temporal low-pass filter')
+    log.info('Applying temporal low-pass filter')
     nanmat = ~isnan(tsincr)
     tsfilt_incr = np.empty_like(tsincr, dtype=np.float32) * np.nan
     intv = np.diff(epochlist.spans)  # time interval for the neighboring epoch
@@ -298,7 +300,7 @@ def temporal_low_pass_filter(tsincr, epochlist, params):
         func = mean_filter
 
     _tlpfilter(cols, cutoff, nanmat, rows, span, threshold, tsfilt_incr, tsincr, func)
-    log.debug("Finished applying temporal low pass filter")
+    log.debug("Finished applying temporal low-pass filter")
     return tsfilt_incr
 
 # Throwaway function to define Gaussian filter weights