[issue_153] all remaining files merged with some changes from original

Author: basaks

pyrate/covariance.py

@@ -20,6 +20,7 @@
 the Matlab Pirate package.
 """
 from __future__ import print_function
+from os.path import basename, join
 import logging
 from numpy import array, where, isnan, real, imag, sqrt, meshgrid
 from numpy import zeros, vstack, ceil, mean, exp, reshape
@@ -31,6 +32,8 @@
 from pyrate import shared
 from pyrate.shared import PrereadIfg
 from pyrate.algorithm import master_slave_ids
+from pyrate import ifgconstants as ifc
+from pyrate import config as cf
 
 # distance division factor of 1000 converts to km and is needed to match
 # Matlab code output
@@ -64,17 +67,22 @@ def unique_points(points):  # pragma: no cover
     return vstack([array(u) for u in set(points)])
 
 
-def cvd(ifg_path, params, r_dist, calc_alpha=False):
+def cvd(ifg_path, params, r_dist, calc_alpha=False,
+        write_vals=False, save_acg=False):
     """
-    Calculate average covariance versus distance (autocorrelation) and its
-    best fitting exponential function
+    Calculate the 1D covariance function of an entire interferogram as the
+    radial average of its 2D autocorrelation.
 
     :param ifg_path: An interferogram.
         ifg: :py:class:`pyrate.shared.Ifg`.
-    :param: params: dict
-        dict of config params
+    :param params: dict
+        dictionary of configuration parameters
     :param calc_alpha: bool
-        whether you calculate alpha.
+        calculate alpha, the exponential length-scale of decay factor.
+    :param write_vals: bool
+        write maxvar and alpha values to interferogram metadata.
+    :param save_acg: bool
+        write acg and radial distance data to numpy array
     """
 
     if isinstance(ifg_path, str):  # used during MPI
@@ -91,15 +99,36 @@ def cvd(ifg_path, params, r_dist, calc_alpha=False):
     else:
         phase = ifg.phase_data
 
-    maxvar, alpha = cvd_from_phase(phase, ifg, r_dist, calc_alpha)
+    maxvar, alpha = cvd_from_phase(phase, ifg, r_dist, calc_alpha,
+                                   save_acg=save_acg)
+
+    if write_vals:
+        _add_metadata(ifg, maxvar, alpha)
 
     if isinstance(ifg_path, str):
         ifg.close()
 
     return maxvar, alpha
 
 
-def cvd_from_phase(phase, ifg, r_dist, calc_alpha):
+def _add_metadata(ifg, maxvar, alpha):
+    """convenience function for saving metadata to ifg"""
+    md = ifg.meta_data
+    md[ifc.PYRATE_MAXVAR] = str(maxvar) #.astype('str')
+    md[ifc.PYRATE_ALPHA] = str(alpha) #.astype('str')
+    ifg.write_modified_phase()
+
+
+def _save_cvd_data(acg, r_dist, ifg_path, outdir):
+    """ function to save numpy array of autocorrelation data to disk"""
+    data = np.column_stack((acg, r_dist))
+    data_file = join(outdir, 'cvd_data_{b}.npy'.format(
+        b=basename(ifg_path).split('.')[0]))
+    np.save(file=data_file, arr=data)
+
+
+def cvd_from_phase(phase, ifg, r_dist, calc_alpha, save_acg=False,
+                   params=None):
     """
     A convenience class reused in many places to compute cvd from phase data
     and a ifg class
@@ -110,6 +139,8 @@ def cvd_from_phase(phase, ifg, r_dist, calc_alpha):
     ifg: shared.Ifg class instance
     calc_alpha: bool, optional
         whether alpha is required
+    save_acg: bool, optional
+        write acg and radial distance data to numpy array
 
     Return
     ------
@@ -150,6 +181,10 @@ def cvd_from_phase(phase, ifg, r_dist, calc_alpha):
     indices_to_keep = r_dist < maxdist
     acg = acg[indices_to_keep]
 
+    # optionally save acg vs dist observations to disk
+    if save_acg:
+        _save_cvd_data(acg, r_dist, ifg.data_path, params[cf.TMPDIR])
+
     if calc_alpha:
         # bin width for collecting data
         bin_width = max(ifg.x_size, ifg.y_size) * 2 / DISTFACT  # km

pyrate/ref_phs_est.py

@@ -38,7 +38,7 @@ def estimate_ref_phase(ifgs, params, refpx, refpy):
         :ref_phs: reference phase correction
         :ifgs: reference phase data removed list of ifgs
     """
-    _check_ref_phs_ifgs(ifgs)
+    check_ref_phs_ifgs(ifgs)
 
     # set reference phase as the average of the whole image (recommended)
     if params[cf.REF_EST_METHOD] == 1:
@@ -135,12 +135,10 @@ def est_ref_phs_method1(phase_data, comp):
     """convenience function for ref phs estimate method 1 parallelisation"""
     ifgv = np.ravel(phase_data, order='F')
     ifgv[comp == 1] = np.nan
-    # reference phase
-    ref_ph = nanmedian(ifgv)
-    return ref_ph
+    return nanmedian(ifgv)
 
 
-def _check_ref_phs_ifgs(ifgs, preread_ifgs=None):
+def check_ref_phs_ifgs(ifgs, preread_ifgs=None):
     """
     Function to check that the ref phase status of all ifgs are the same
     """
@@ -149,6 +147,7 @@ def _check_ref_phs_ifgs(ifgs, preread_ifgs=None):
         raise ReferencePhaseError('Need to provide at least 2 ifgs')
 
     if preread_ifgs:  # check unless for mpi tests
+        # preread_ifgs[i].metadata contains ifg metadata
         flags = [ifc.PYRATE_REF_PHASE in preread_ifgs[i].metadata
                  for i in ifgs]
     else:

pyrate/scripts/run_pyrate.py

@@ -120,7 +120,7 @@ def create_ifg_dict(dest_tifs, params, tiles):
     mpiops.comm.barrier()
     preread_ifgs = OrderedDict(sorted(cp.load(open(preread_ifgs_file,
                                                    'rb')).items()))
-    log.info('finish converting phase_data to numpy '
+    log.info('Finished converting phase_data to numpy '
              'in process {}'.format(mpiops.rank))
     return preread_ifgs
 
@@ -270,14 +270,16 @@ def save_ref_pixel_blocks(grid, half_patch_size, ifg_paths, params):
 
 def orb_fit_calc(ifg_paths, params, preread_ifgs=None):
     """
-    Orbital fit correction
+    MPI wrapper for orbital fit correction
 
     Parameters
     ----------
     ifg_paths: list
         list of ifg paths
     params: dict
         parameters dict corresponding to config file
+    preread_ifgs: dict, optional
+        dict containing information regarding ifgs
     """
     log.info('Calculating orbfit correction')
 
@@ -300,15 +302,15 @@ def orb_fit_calc(ifg_paths, params, preread_ifgs=None):
         # Here we do all the multilooking in one process, but in memory
         # can use multiple processes if we write data to disc during
         # remove_orbital_error step
-        # A performance comparison should be performed be performed for saving
-        # multilooked files on disc vs in memory single process multilooking
+        # A performance comparison should be made for saving multilooked
+        # files on disc vs in memory single process multilooking
         if mpiops.rank == MASTER_PROCESS:
             orbital.remove_orbital_error(ifg_paths, params, preread_ifgs)
     mpiops.comm.barrier()
-    log.info('Finished orbfit calculation in process {}'.format(mpiops.rank))
+    log.info('Finished orbfit calculation')
 
 
-def ref_phase_estimation(ifg_paths, params, refpx, refpy):
+def ref_phase_estimation(ifg_paths, params, refpx, refpy, preread_ifgs=None):
     """
     Reference phase estimation
 
@@ -322,18 +324,27 @@ def ref_phase_estimation(ifg_paths, params, refpx, refpy):
         reference pixel x-coordinate
     refpy: float
         reference pixel y-coordinate
+    preread_ifgs: dict
+        dict containing information regarding ifgs
     """
+    # perform some checks on existing ifgs
+    if preread_ifgs and mpiops.rank == MASTER_PROCESS:
+        ifg_paths = sorted(preread_ifgs.keys())
+        if mpiops.run_once(rpe.check_ref_phs_ifgs, ifg_paths, preread_ifgs):
+            return  # return if True condition returned
 
-    log.info('Estimating and removing reference phase')
     if params[cf.REF_EST_METHOD] == 1:
         # calculate phase sum for later use in ref phase method 1
         comp = phase_sum(ifg_paths, params)
+        log.info('Computing reference phase via method 1')
         process_ref_phs = ref_phs_method1(ifg_paths, comp)
     elif params[cf.REF_EST_METHOD] == 2:
+        log.info('Computing reference phase via method 2')
         process_ref_phs = ref_phs_method2(ifg_paths, params, refpx, refpy)
     else:
         raise ConfigException('Ref phase estimation method must be 1 or 2')
 
+    # Save reference phase numpy arrays to disk
     ref_phs_file = join(params[cf.TMPDIR], 'ref_phs.npy')
     if mpiops.rank == MASTER_PROCESS:
         ref_phs = np.zeros(len(ifg_paths), dtype=np.float64)
@@ -350,6 +361,7 @@ def ref_phase_estimation(ifg_paths, params, refpx, refpy):
         # send reference phase data to master process
         mpiops.comm.Send(process_ref_phs, dest=MASTER_PROCESS,
                          tag=mpiops.rank)
+    log.info('Completed reference phase estimation')
 
 
 def ref_phs_method2(ifg_paths, params, refpx, refpy):
@@ -443,7 +455,7 @@ def process_ifgs(ifg_paths, params, rows, cols):
     cols: int
         number of cols to break each ifg into
     """
-    if mpiops.size > 1: # turn of multiprocessing during mpi jobs
+    if mpiops.size > 1:  # turn of multiprocessing during mpi jobs
         params[cf.PARALLEL] = False
 
     tiles = mpiops.run_once(get_tiles, ifg_paths[0], rows, cols)
@@ -504,10 +516,10 @@ def linrate_calc(ifg_paths, params, vcmt, tiles, preread_ifgs):
     """
 
     process_tiles = mpiops.array_split(tiles)
-    log.info('Calculating linear rate')
+    log.info('Calculating linear rate map')
     output_dir = params[cf.TMPDIR]
     for t in process_tiles:
-        log.info('calculating lin rate of tile {}'.format(t.index))
+        log.info('Calculating linear rate of tile {}'.format(t.index))
         ifg_parts = [shared.IfgPart(p, t, preread_ifgs) for p in ifg_paths]
         mst_grid_n = np.load(os.path.join(output_dir,
                                           'mst_mat_{}.npy'.format(t.index)))
@@ -641,9 +653,6 @@ def timeseries_calc(ifg_paths, params, vcmt, tiles, preread_ifgs):
         dict containing ifg characteristics for efficient computing
 
     """
-    if not params[cf.TIME_SERIES_CAL]:
-        log.info('Time series calculation not required.')
-        return
     process_tiles = mpiops.array_split(tiles)
     log.info('Calculating time series')
     output_dir = params[cf.TMPDIR]