diff --git a/.travis.yml b/.travis.yml index 12c95c37a..7375824fa 100644 --- a/.travis.yml +++ b/.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 diff --git a/docs/history.rst b/docs/history.rst index ef726373f..e8201cc80 100644 --- a/docs/history.rst +++ b/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. diff --git a/input_parameters.conf b/input_parameters.conf index 5514bc12f..c10f55628 100644 --- a/input_parameters.conf +++ b/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 diff --git a/pyrate/configuration.py b/pyrate/configuration.py index 9927da350..66d0c4c3c 100644 --- a/pyrate/configuration.py +++ b/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] diff --git a/pyrate/constants.py b/pyrate/constants.py index 7309ca12c..6937de552 100644 --- a/pyrate/constants.py +++ b/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 diff --git a/pyrate/conv2tif.py b/pyrate/conv2tif.py index db988f6dc..b7572683f 100644 --- a/pyrate/conv2tif.py +++ b/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!") diff --git a/pyrate/core/aps.py b/pyrate/core/aps.py index 84ee106c0..2ffd46309 100644 --- a/pyrate/core/aps.py +++ b/pyrate/core/aps.py @@ -19,7 +19,6 @@ signals. """ # pylint: disable=invalid-name, too-many-locals, too-many-arguments -import logging import os from copy import deepcopy from collections import OrderedDict @@ -27,15 +26,14 @@ 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): @@ -43,8 +41,10 @@ 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 diff --git a/pyrate/core/config.py b/pyrate/core/config.py index b73d66934..09741be1b 100644 --- a/pyrate/core/config.py +++ b/pyrate/core/config.py @@ -16,8 +16,7 @@ """ This Python module contains utilities to parse PyRate configuration files. It also includes numerous general constants relating to options -in configuration files. Examples of PyRate configuration files are -provided in the configs/ directory +in configuration files. """ # coding: utf-8 # pylint: disable=invalid-name @@ -295,7 +294,7 @@ INT_KEYS = [APS_CORRECTION, APS_METHOD] -def get_config_params(path: str, validate: bool=True, step: str=CONV2TIF) -> Dict: +def get_config_params(path: str) -> Dict: """ Reads the parameters file provided by the user and converts it into a dictionary. @@ -317,13 +316,13 @@ def get_config_params(path: str, validate: bool=True, step: str=CONV2TIF) -> Dic # create expanded line line = line[:pos] + os.environ['HOME'] + line[(pos+1):] txt += line - params = _parse_conf_file(txt, validate, step) + params = _parse_conf_file(txt) params[TMPDIR] = os.path.join(os.path.abspath(params[OUT_DIR]), 'tmpdir') return params -def _parse_conf_file(content, validate: bool=True, step: str=CONV2TIF) -> Dict: +def _parse_conf_file(content) -> Dict: """ Converts the parameters from their text form into a dictionary. @@ -358,7 +357,7 @@ def _is_valid(line): if not parameters: raise ConfigException('Cannot parse any parameters from config file') - return _parse_pars(parameters, validate, step) + return _parse_pars(parameters) def _handle_extra_parameters(params): @@ -383,7 +382,7 @@ def _handle_extra_parameters(params): return params -def _parse_pars(pars, validate: bool=True, step: str=CONV2TIF) -> Dict: +def _parse_pars(pars) -> Dict: """ Takes dictionary of parameters, converting values to required type and providing defaults for missing values. @@ -414,8 +413,6 @@ def _parse_pars(pars, validate: bool=True, step: str=CONV2TIF) -> Dict: if pars.get(p) is None: pars[p] = pars[OBS_DIR] - if validate: - validate_parameters(pars, step) return pars @@ -433,6 +430,7 @@ def parse_namelist(nml): lines = [line.rstrip() for line in f_in] return filter(None, lines) + def write_config_file(params, output_conf_file): """ Takes a param object and writes the config file. Reverse of get_conf_params. @@ -450,6 +448,7 @@ def write_config_file(params, output_conf_file): else: f.write(''.join([k, ':\t', '', '\n'])) + def transform_params(params): """ Returns subset of all parameters for cropping and multilooking. @@ -547,41 +546,6 @@ def get_dest_paths(base_paths, crop, params, looks): return [os.path.join(params[OUT_DIR], p) for p in dest_mlooked_ifgs] -def get_ifg_paths(config_file, step=CONV2TIF): - """ - Read the configuration file, extract interferogram file list and determine - input and output interferogram path names. - - :param str config_file: Configuration file path - - :return: base_unw_paths: List of unwrapped inteferograms - :return: dest_paths: List of multi-looked and cropped geotifs - :return: params: Dictionary corresponding to the config file - :rtype: list - :rtype: list - :rtype: dict - """ - params = get_config_params(config_file, step=step) - ifg_file_list = params.get(IFG_FILE_LIST) - - xlks, _, crop = transform_params(params) - - # base_unw_paths need to be geotiffed by conv2tif and multilooked by prepifg - base_unw_paths = original_ifg_paths(ifg_file_list, params[OBS_DIR]) - - # dest_paths are tifs that have been coherence masked (if enabled), - # cropped and multilooked - - if "tif" in base_unw_paths[0].split(".")[1]: - - dest_paths = get_dest_paths(base_unw_paths, crop, params, xlks) - for i, dest_path in enumerate(dest_paths): - dest_paths[i] = dest_path.replace("_tif", "") - else: - dest_paths = get_dest_paths(base_unw_paths, crop, params, xlks) - - return base_unw_paths, dest_paths, params - # ==== PARAMETER VALIDATION ==== # _PARAM_VALIDATION = { @@ -818,832 +782,7 @@ def get_ifg_paths(config_file, step=CONV2TIF): f"'{REF_EST_METHOD}': must select option 1 or 2." ), } -"""dict: basic validation functions for reference pixel search parameters.""" - -def convert_geographic_coordinate_to_pixel_value(refpx, refpy, transform): - """ - Converts a lat/long coordinate to a pixel coordinate given the - geotransform of the image. - - Args: - refpx: The longitude of the coordinate. - refpx: The latitude of the coordinate. - transform: The geotransform array of the image. - - Returns: - Tuple of refpx, refpy in pixel values. - """ - # transform = ifg.dataset.GetGeoTransform() - - xOrigin = transform[0] - yOrigin = transform[3] - pixelWidth = transform[1] - pixelHeight = -transform[5] - - refpx = int((refpx - xOrigin) / pixelWidth) - refpy = int((yOrigin - refpy) / pixelHeight) - - return int(refpx), int(refpy) - - -def validate_parameters(pars: Dict, step: str=CONV2TIF): - """ - Main validation function. Calls validation subfunctions and gathers - some required variables for performing validation. - - Args: - pars: The parameters dictionary. - step: The current step of the PyRate workflow. Determines what - parameters to validate and what resources are available in - terms of geotiffs. - - Raises: - ConfigException: If errors occur during parameter validation. - """ - is_GAMMA = pars[PROCESSOR] == GAMMA - ifl = pars[IFG_FILE_LIST] - - # TODO: Call bounds checking functions based on the current step - # Call basic bounds checking functions for all parameters. - validate_compulsory_parameters(pars) - validate_optional_parameters(pars) - - # Validate that provided filenames contain correct epochs and that - # the files exist. - if is_GAMMA: - validate_epochs(ifl, SIXTEEN_DIGIT_EPOCH_PAIR) - validate_epochs(pars[HDR_FILE_LIST], EIGHT_DIGIT_EPOCH) - validate_gamma_headers(ifl, pars[HDR_FILE_LIST]) - else: - validate_epochs(ifl, TWELVE_DIGIT_EPOCH_PAIR) - - if step == CONV2TIF: - # Check that unwrapped interferograms exist. - validate_ifgs(ifl) - - elif step == PREPIFG: - # Check full res geotiffs exist before continuing. - validate_tifs_exist(ifl, pars[OUT_DIR]) - - # Check the minimum number of epochs. - n_epochs = 0 - with open(ifl, "r") as f: - list_of_epoches = [] - for line in f.readlines(): - - PTN = re.compile(r'\d{8}') # match 8 digits for the dates - - # if ROI_PAC - if 0 == pars["processor"]: - PTN = re.compile(r'\d{6}') # match 8 digits for the dates - - epochs = PTN.findall(line.strip()) - list_of_epoches.extend(epochs) - - list_of_epoches = set(list_of_epoches) - n_epochs = len(list_of_epoches) - validate_minimum_epochs(n_epochs, MINIMUM_NUMBER_EPOCHS) - - # validate - with open(ifl, "r") as f: - # validate params for each geotiff - for line in f.readlines(): - line = line.strip('\n') - if line.endswith('.tif'): - tif_file_path = Path(pars[OUT_DIR]).joinpath(Path(line.strip()).name) - else: - p = Path(line) - base, ext = p.stem, p.suffix[1:] - tif_file_path = Path(pars[OUT_DIR]).joinpath(base + "_" + ext + ".tif") - - if not tif_file_path.exists(): - raise Exception("GeoTiff: " + tif_file_path.as_posix() + " not found.") - - raster = os.path.join(tif_file_path) - gtif = gdal.Open(raster) - - latitudes = [] - longitudes = [] - - for line in gdal.Info(gtif).split('\n'): - if "Size is" in line: - x_size, y_size = line.split("Size is")[1].split(",") - x_size, y_size = int(x_size.strip()), int(y_size.strip()) - - for line_tag in ["Upper Left", "Lower Left", "Upper Right", "Lower Right"]: - if line_tag in line: - latitude, longitude = line.split(")")[0].split("(")[1].split(",") - latitudes.append(float(latitude.strip())) - longitudes.append(float(longitude.strip())) - - # validate multi-look parameters - if pars["ifglksx"] < 0 or pars["ifglksx"] > x_size: - raise Exception("Value of ifglksx: "+str(pars["ifglksx"])+" out of bounds: [0,"+str(x_size)+"]") - if pars["ifglksy"] < 0 or pars["ifglksy"] > x_size: - raise Exception("Value of ifglksy: "+str(pars["ifglksy"])+" out of bounds: [0,"+str(x_size)+"]") - - # validate crop parameters - if pars["ifgxfirst"] < min(latitudes) or pars["ifgxfirst"] > max(latitudes): - raise Exception( - "ifgxfirst: " + str(pars["ifgxfirst"]) + " not with in range {" + str( - min(latitudes)) + "," + str(max(latitudes)) + "}" - ) - - if pars["ifgxlast"] < min(latitudes) or pars["ifgxlast"] > max(latitudes): - raise Exception( - "ifgxlast: " + str(pars["ifgxlast"]) + " not with in range {" + str(min(latitudes)) + "," + str( - max(latitudes)) + "}" - ) - - if pars["ifgyfirst"] < min(longitudes) or pars["ifgyfirst"] > max(longitudes): - raise Exception( - "ifgyfirst: " + str(pars["ifgyfirst"]) + " not with in range {" + str( - min(longitudes)) + "," + str(max(longitudes)) + "}" - ) - - if pars["ifgylast"] < min(longitudes) or pars["ifgylast"] > max(longitudes): - raise Exception( - "ifgylast: " + str(pars["ifgylast"]) + " not with in range {" + str( - min(longitudes)) + "," + str(max(longitudes)) + "}" - ) - - del gtif # manually close raster - - # Check coherence masking if enabled - if pars[COH_MASK]: - validate_epochs(pars[COH_FILE_LIST], SIXTEEN_DIGIT_EPOCH_PAIR) - validate_coherence_files(ifl, pars) - - elif step == PROCESS: - - tmp_ifg_list = os.path.join(pars[OUT_DIR], "tmpIfgList") - with open(tmp_ifg_list, "w") as fileHandler: - for p in Path(pars[OUT_DIR]).rglob("*rlks_*cr.tif"): - if "dem" not in str(p): - fileHandler.write(p.name + "\n") - # Check that cropped/subsampled tifs exist. - validate_prepifg_tifs_exist(tmp_ifg_list, pars[OUT_DIR], pars) - - # Check the minimum number of epochs. - n_epochs, max_span = _get_temporal_info(tmp_ifg_list, pars[OBS_DIR]) - validate_minimum_epochs(n_epochs, MINIMUM_NUMBER_EPOCHS) - - # Get spatial information from tifs. - - extents, n_cols, n_rows, transform = _get_prepifg_info(tmp_ifg_list, pars[OBS_DIR], pars) - - # test if refx/y already set to default value of -1 - if pars[REFX] != -1 and pars[REFY] != -1: - - # Convert refx/refy from lat/long to pixel and validate... - if (pars[REFX] <= 180) and (pars[REFX] >= -180) and (pars[REFY] >= -90) and (pars[REFY] <= 90): - - pars[REFX], pars[REFY] = convert_geographic_coordinate_to_pixel_value(pars[REFX], pars[REFY], transform) - _logger.debug("converted pars[REFX], pars[REFY] to: " + str(pars[REFX]) + " " + str(pars[REFY])) - - if pars[REFX] < 0 or pars[REFX] > n_cols: - _logger.info("converted pars[REFX] out of range") - pars[REFX] = -1 - if pars[REFY] < 0 or pars[REFY] > n_rows: - _logger.info("converted pars[REFY] out of range") - pars[REFY] = -1 - - # otherwise we need to search for the pixel so validate the search parameters. - else: - _logger.info("given pars[REFX] or pars[REFY] out of range") - pars[REFX] = -1 - pars[REFY] = -1 - # TODO fix tests when below validation is included - # validate_reference_pixel_search_windows(n_cols, n_rows, pars) - - validate_multilook_parameters(n_cols, n_rows, ORBITAL_FIT_LOOKS_X, ORBITAL_FIT_LOOKS_Y, pars) - validate_slpf_cutoff(extents, pars) - validate_epoch_thresholds(n_epochs, pars) - validate_epoch_cutoff(max_span, TLPF_CUTOFF, pars) - validate_obs_thresholds(tmp_ifg_list, pars) - - elif step == MERGE: - validate_prepifg_tifs_exist(ifl, pars[OBS_DIR], pars) - - -def _raise_errors(errors: List[str]) -> bool: - """ - Convenience function for raising an exception with errors. - """ - if errors: - errors.insert(0, "invalid parameters") - raise ConfigException('\n'.join(errors)) - return True - -def validate_compulsory_parameters(pars: Dict) -> Optional[bool]: - """ - Calls the validators for compulsory (always used) parameters. - - Args: - pars: The parameters dictionary. - - Returns: - True if validation is successful. - - Raises: - ConfigException: If validation fails. - """ - errors = [] - for k in pars.keys(): - validator = _PARAM_VALIDATION.get(k) - if validator is None: - # _logger.debug(f"No validator implemented for '{k}'.") - continue - if not validator[0](pars[k]): - errors.append(validator[1]) - - return _raise_errors(errors) - -def validate_optional_parameters(pars: Dict): - """ - Calls the validators for optional parameters. - - Args: - pars: The parameters dictionary. - - Returns: - True if validation successful. - - Raises: - ConfigException: If validation fails. - """ - def validate(on: bool, validators: Dict, pars: Dict) -> List[str]: - """ - Convenience method for calling validators. - - Args: - on: Determines whether to call the validators. - validators: A dictionary of validator functions. - pars: Parameters dictionary. - - Returns: - A list of errors. - """ - errors = [] - if on: - for k, validator in validators.items(): - if not validator[0](pars[k]): - errors.append(validator[1]) - return errors - - errors = [] - - errors.extend(validate(pars[COH_MASK], _COHERENCE_VALIDATION, pars)) - errors.extend(validate(pars[APSEST], _APSEST_VALIDATION, pars)) - errors.extend(validate(pars[TIME_SERIES_CAL], _TIME_SERIES_VALIDATION, pars)) - errors.extend(validate(pars[ORBITAL_FIT], _ORBITAL_FIT_VALIDATION, pars)) - errors.extend(validate(pars[PROCESSOR] == GAMMA, _GAMMA_VALIDATION, pars)) - errors.extend(validate(pars[IFG_CROP_OPT] == 3, _CUSTOM_CROP_VALIDATION, pars)) - # errors.extend(validate(pars[REFX] > 0 and pars[REFY] > 0, _REFERENCE_PIXEL_VALIDATION, pars)) - - return _raise_errors(errors) - -def validate_minimum_epochs(n_epochs: int, min_epochs: int) -> Optional[bool]: - """ - Validates the minimum number of epochs required for PyRate to produce - good results. - - Args: - n_epochs: The number of unique epochs in the collection of interferograms - provided as input. - min_epochs: The minimum number of epochs PyRate requires. - - Returns: - True if there are enough epochs to satisfy minimum epochs. - - Raises: - ConfigException: If there are not enough epochs to satisfy min epochs. - """ - errors = [] - if n_epochs < min_epochs: - errors.append(f"'{IFG_FILE_LIST}': total number of epochs given is {n_epochs}." - f" {min_epochs} or more unique epochs are required by PyRate.") - _raise_errors(errors) - -def validate_epochs(file_list: str, pattern: str) -> Optional[bool]: - """ - Validate that user provided file names contain the correct pattern of - epochs. - - Args: - file_list: Path to the list of files to validate. - pattern: Regex string for finding the epoch(s). - - Returns: - True if all names in file list contain the epoch pattern *once*. - - Raises: - ConfigException: If not all names in the file list don't contain - the epoch pattern or contain it more than once. - """ - errors = [] - PTN = re.compile(pattern) - filenames = parse_namelist(file_list) - for fn in filenames: - epochs = PTN.findall(fn) - if not epochs: - errors.append(f"'{file_list}': {fn} does not contain an epoch of " - f"format {pattern}.") - if len(epochs) > 1: - errors.append(f"'{file_list}': {fn} does contains more than " - f"one epoch of {pattern}. There must be only one " - f"epoch in the filename.") - - return _raise_errors(errors) - -def validate_epoch_cutoff(max_span: float, cutoff: str, pars: Dict) -> Optional[bool]: - """ - Validate cutoff parameters that rely on the data timespan. - - Args: - max_span: The maximum temporal span of the provided data in years. - cutoff: The key of the cutoff parameter. - pars: The parameters dictionary. - - Returns: - True if the cutoff is less than the maximum data timespan. - - Raises: - ConfigException: If the cutoff is greater than the max data timespan. - """ - errors = [] - if pars[cutoff] > max_span: - errors.append("'{cutoff}': must be less than max time span of " - "data in years ({max_span}).") - return _raise_errors(errors) - -def validate_prepifg_tifs_exist(ifg_file_list: str, obs_dir: str, pars: Dict) -> Optional[bool]: - """ - Validates that cropped and multilooked interferograms exist in geotiff - format. - - Args: - ifg_file_list: Path to file containing interfergram file names. - obs_dir: Path to observations directory. - pars: Parameters dictionary. - - Returns: - True if all interferograms exist in geotiff format. - - Raises: - ConfigException: If not all intergerograms exist in geotiff format. - """ - - errors = [] - base_paths = [os.path.join(obs_dir, ifg) for ifg in parse_namelist(ifg_file_list)] - ifg_paths = get_dest_paths(base_paths, pars[IFG_CROP_OPT], pars, pars[IFG_LKSX]) - for i, ifg_path in enumerate(ifg_paths): - ifg_paths[i] = ifg_path.replace("_tif", "") - - for path in ifg_paths: - if not os.path.exists(path): - fname = os.path.split(path)[1] - errors.append(f"'{IFG_FILE_LIST}': interferogram '{fname}' is " - f"required as a cropped and subsampled geotiff but " - f"could not be found. Make sure 'prepifg' has been " - f"run and ensure the '{IFG_LKSX}' and '{IFG_CROP_OPT}' " - f"parameters have not been changed since 'prepifg' was run.") - - return _raise_errors(errors) - - -def validate_tifs_exist(ifg_file_list: str, obs_dir: str) -> Optional[bool]: - """ - Validates that provided interferograms exist in geotiff format. - - Args: - ifg_file_list: Path to file containing interfergram file names. - obs_dir: Path to observations directory. - - Returns: - True if all interferograms exist in geotiff format. - - Raises: - ConfigException: If not all intergerograms exist in geotiff format. - """ - from pyrate.core.shared import output_tiff_filename - - errors = [] - ifgs = parse_namelist(ifg_file_list) - ifg_paths = [os.path.join(obs_dir, ifg) for ifg in ifgs] - gtiff_paths = [output_tiff_filename(f, obs_dir) for f in ifg_paths] - for gtp in gtiff_paths: - if not os.path.exists(gtp): - fname = os.path.split(gtp)[1] - errors.append(f"'{IFG_FILE_LIST}': interferogram '{fname}' is " - "required in geotiff format but no geotiff file " - "could be found.") - - return _raise_errors(errors) - - -def validate_ifgs(ifg_file_list: str) -> Optional[bool]: - """ - Validates that provided interferograms exist. - - Args: - ifg_file_list: Path to file containing interferogram file names. - obs_dir: Path to observations directory. - - Returns: - True if all interferograms exist. - - Raises: - ConfigException: If not all interferograms exist. - """ - errors = [] - ifgs = parse_namelist(ifg_file_list) - for path in ifgs: - if not os.path.exists(path): - fname = os.path.split(path)[1] - errors.append(f"'{IFG_FILE_LIST}': interferogram '{fname}' does not exist.") - - return _raise_errors(errors) - - -def validate_obs_thresholds(ifg_file_list: str, pars: Dict) -> Optional[bool]: - """ - Validates parameters that specify an observations threshold. - - Args: - ifg_file_list: Path to the file containing interferogram file names. - pars: Parameters dictionary. - - Returns: - True if there are enough interferograms to satisfy all observation thresholds. - - Raises: - ConfigException: If there not enough interferograms to satisfy all observation - thresholds. - """ - def validate(n, p, k): - thresh = p[k] - if thresh > n: - return [f"'{k}': not enough interferograms have been specified " - f"({n}) to satisfy threshold ({thresh})."] - return [] - - errors = [] - n_ifgs = len(list(parse_namelist(ifg_file_list))) - errors.extend(validate(n_ifgs, pars, TIME_SERIES_PTHRESH)) - if pars[APSEST]: - errors.extend(validate(n_ifgs, pars, TLPF_PTHR)) - - return _raise_errors(errors) - -def validate_epoch_thresholds(n_epochs: int, pars: Dict) -> Optional[bool]: - """ - Validates threshold paramters that rely on the number of epochs - available. - - Args: - n_epochs: The number of unique epochs in the collection of interferograms - provided as input. - pars: Parameters dictionary. - - Returns: - True if there are enough epochs to satisfy all epoch thresholds. - - Raises: - ConfigException: If there are not enough epochs to satisfy all epoch thresholds. - """ - errors = [] - thresh = pars[LR_PTHRESH] - if n_epochs < thresh: - errors.append(f"'{LR_PTHRESH}': not enough epochs have been specified " - f"({n_epochs}) to satisfy threshold ({thresh}).") - - return _raise_errors(errors) - -def validate_crop_parameters(min_extents: Tuple[float, float, float, float], - pars: Dict) -> Optional[bool]: - """ - Validates IFG crop parameters by ensuring user provided crop coordinates - are within the *minimum* (intersecting) bounds of the interferogram stack. - - Args: - min_extents: The minimum extents of the interferogram stack. The - crop coordinates must be within these extents. - pars: Parameters dictionary. - - Returns: - True if validation is successful. - - Raises: - ConfigException: If validation fails. - """ - errors = [] - min_xmin, min_ymin, min_xmax, min_ymax = min_extents - x_dim_string = f"(xmin: {min_xmin}, xmax: {min_xmax})" - y_dim_string = f"(ymin: {min_ymin}, ymax: {min_ymax})" - - # Check crop coordinates within scene. - def _validate_crop_coord(var_name, dim_min, dim_max, dim_string): - if not dim_min < pars[var_name] < dim_max: - return [f"'{var_name}': crop coordinate ({pars[var_name]}) " - f"is outside bounds of scene {dim_string}."] - - return [] - - errors.extend(_validate_crop_coord(IFG_XFIRST, min_xmin, min_xmax, x_dim_string)) - errors.extend(_validate_crop_coord(IFG_YFIRST, min_ymin, min_ymax, y_dim_string)) - errors.extend(_validate_crop_coord(IFG_XLAST, min_xmin, min_xmax, x_dim_string)) - errors.extend(_validate_crop_coord(IFG_YLAST, min_ymin, min_ymax, y_dim_string)) - - return _raise_errors(errors) - - -def validate_slpf_cutoff(extents: Tuple[float, float, float, float], - pars: Dict) -> Optional[bool]: - """ - Validate SLPF_CUTOFF is within the bounds of the scene being processed - (after prepifg cropping/subsampling has occurred). - - Args: - extents : Tuple of (xmin, xmax, ymin, ymax) describing the extents - of the scene, after cropping & multisampling, in degrees. - pars: Parameters dictionary. - - Returns: - True if validation is successful. - - Raises: - ConfigException: If validation fails. - """ - xmin, ymin, xmax, ymax = extents - # Check SLPF_CUTOFF within scene *in kilometeres*. - DEG_TO_KM = 111.32 # km per degree - x_extent = abs(xmin - xmax) - y_extent = abs(ymin - ymax) - x_extent_km = x_extent * DEG_TO_KM - y_extent_km = y_extent * DEG_TO_KM - if pars[SLPF_CUTOFF] > max(x_extent_km, y_extent_km): - errors = [f"'{SLPF_CUTOFF}': cutoff is out of bounds, must be " - "less than max scene bound (in km) " - f"({max(x_extent_km, y_extent_km)})."] - else: - errors = [] - - return _raise_errors(errors) - - - -def validate_multilook_parameters(cols: int, rows: int, - xlks_name: str, ylks_name: str, - pars: Dict) -> Optional[bool]: - """ - Validate multilook parameters by ensuring resulting resolution will be - at least 1 pixel and less than the current resolution. - - Args: - cols: Number of pixel columns (X) in the scene. - rows: Number of pixel rows (Y) in the scene. - xlks_name: Key for looking up the X looks factor. - ylks_name: Key for looking up the Y looks factor. - pars: Parameters dictionary. - - Returns: - True if validation is successful. - - Raises: - ConfigException: If validation fails. - """ - errors = [] - - def _validate_mlook(var_name, dim_val, dim_str): - if pars[var_name] > dim_val: - return [f"'{var_name}': the multilook factor ({pars[var_name]}) " - f"must be less than the number of pixels on the {dim_str} " - f"axis ({dim_val})."] - return [] - - errors.extend(_validate_mlook(xlks_name, cols, 'x')) - errors.extend(_validate_mlook(ylks_name, rows, 'y')) - return _raise_errors(errors) - -def validate_reference_pixel_search_windows(n_cols: int, n_rows: int, - pars: Dict) -> Optional[bool]: - """ - Validates that the reference pixel search windows provided by user - fit within the scene being processed. - - Args: - n_cols: Number of pixel columns (X) in the scene. - n_rows: Number of pixel rows (Y) in the scene. - - Returns: - True if the scene can accomodate the search windows (no overlap). - - Raises: - ConfigException: If the scene cannot accomodate the search windows. - """ - from math import floor - - errors = [] - refnx = pars[REFNX] - refny = pars[REFNY] - chip_size = pars[REF_CHIP_SIZE] - - x_windows = floor(n_cols/chip_size) - if refnx > x_windows: - errors.append(f"'{REFNX}' & '{REF_CHIP_SIZE}': search windows do not " - f"fit in scene on X axis (number of columns: {n_cols}). " - f"Reduce {REF_CHIP_SIZE} or {REFNX} so that {REFNX} " - f"is less than or equal to (columns / {REF_CHIP_SIZE}).") - y_windows = floor(n_rows/chip_size) - if refny > y_windows: - errors.append(f"'{REFNY}' & '{REF_CHIP_SIZE}': search windows do not " - f"fit in scene on Y axis (number of rows: {n_rows}). " - f"Reduce {REF_CHIP_SIZE} or {REFNY} so that {REFNY} " - f"is less than or equal to (rows / {REF_CHIP_SIZE}).") - - return _raise_errors(errors) - - -def validate_gamma_headers(ifg_file_list: str, slc_file_list: str) -> Optional[bool]: - """ - Validates that a pair of GAMMA headers exist for each provided - GAMMA interferogram. - - Args: - ifg_file_list: Path to the file listing filenames of interferograms. - slc_file_list: Path to the file listing filenames of GAMMA headers. - slc_dir: Path to directory containing GAMMA headers. - - Returns: - True if there are exactly 2 headers for each interogram (one for - each epoch). - - Raises: - ConfigException: If there are 0 or more than 2 matching headers - for an interferogram. - """ - from pyrate.core.gamma import get_header_paths - errors = [] - - for ifg in parse_namelist(ifg_file_list): - headers = get_header_paths(ifg, slc_file_list) - if len(headers) < 2: - errors.append(f"'{SLC_DIR}': Headers not found for interferogram '{ifg}'.") - - return _raise_errors(errors) - - -def validate_coherence_files(ifg_file_list: str, pars: Dict) -> Optional[bool]: - """ - Validates that there is a matching coherence file for each provided - interferogram. - - Args: - ifg_file_list: Path to file listing interferogram names. - pars: The parameters dictionary. - - Returns: - True if there is exactly 1 coherence file for each interferogram. - - Raises: - ConfigException: If there are 0 or more than 1 matching coherence - files for an interferogram. - """ - errors = [] - - for ifg in parse_namelist(ifg_file_list): - paths = coherence_paths_for(ifg, pars) - if not paths: - errors.append(f"'{COH_FILE_DIR}': no coherence files found for intergerogram '{ifg}'.") - return _raise_errors(errors) - - -def _get_temporal_info(ifg_file_list: str, obs_dir: str) -> Tuple: - """ - Retrieves the number of unique epochs and maximum time span of the - dataset. - - Args: - ifg_file_list: Path to file containing list of interferogram file names. - obs_dir: Path to observations directory. - - Returns: - Tuple containing the number of unique epochs and the maximum timespan. - """ - from pyrate.core.algorithm import get_epochs - from pyrate.core.shared import Ifg, output_tiff_filename - - ifg_paths = [os.path.join(obs_dir, ifg) for ifg in parse_namelist(ifg_file_list)] - rasters = [Ifg(output_tiff_filename(f, obs_dir)) for f in ifg_paths] - - for r in rasters: - if not r.is_open: - r.open(readonly=True) - - # extents = xmin, ymin, xmax, ymax - epoch_list = get_epochs(rasters)[0] - n_epochs = len(epoch_list.dates) - max_span = max(epoch_list.spans) - - for r in rasters: - if not r.is_open: - r.close() - - return n_epochs, max_span - - -def _get_prepifg_info(ifg_file_list: str, obs_dir: str, pars: Dict) -> Tuple: - """ - Retrives spatial information from prepifg interferograms (images that - have been cropped and multilooked). - - Args: - ifg_file_list: Path to file containing list of interferogram file names. - obs_dir: Path to observations directory. - - Returns: - """ - from pyrate.core.shared import Ifg - base_paths = [os.path.join(obs_dir, ifg) for ifg in parse_namelist(ifg_file_list)] - ifg_paths = get_dest_paths(base_paths, pars[IFG_CROP_OPT], pars, pars[IFG_LKSX]) - - for i, ifg_path in enumerate(ifg_paths): - ifg_paths[i] = ifg_path.replace("_tif","") - - # This function assumes the stack of interferograms now have the same - # bounds and resolution after going through prepifg. - raster = Ifg(ifg_paths[0]) - if not raster.is_open: - raster.open(readonly=True) - - n_cols = raster.ncols - n_rows = raster.nrows - extents = raster.x_first, raster.y_first, raster.x_last, raster.y_last - transform = raster.dataset.GetGeoTransform() - - raster.close() - - return extents, n_cols, n_rows, transform - - -def _get_fullres_info(ifg_file_list: str, out_dir: str, crop_opts: Tuple) -> Tuple: - """ - Retrieves spatial information from the provided interferograms. - Requires the interferograms to exist in geotiff format. - - Args: - ifg_file_list: Path to file containing list of interferogram file names. - out_dir: Path to observations directory. - - Returns: - Tuple containing extents (xmin, ymin, xmax, ymax), number of pixel - columns, number of pixel rows, number of unique epochs and maximum - time span of the data. - """ - from pyrate.core.prepifg_helper import _min_bounds, _get_extents - from pyrate.core.shared import Ifg, output_tiff_filename - - ifg_paths = parse_namelist(ifg_file_list) - rasters = [Ifg(output_tiff_filename(f, out_dir)) for f in ifg_paths] - - for r in rasters: - if not r.is_open: - r.open(readonly=True) - - # extents = xmin, ymin, xmax, ymax - min_extents = _min_bounds(rasters) - post_crop_extents = _get_extents(rasters, crop_opts[0], user_exts=crop_opts[1]) - x_step = rasters[0].x_step - y_step = rasters[0].y_step - # Get the pixel bounds. Ifg/Raster objects do have 'ncols'/'nrows' - # properties, but we'll calculate it off the extents we got above - # because these take into account the chosen cropping option (until - # the stack of interferograms is cropped it's not known what the - # pixel dimensions will be). - n_cols = abs(int(abs(post_crop_extents[0] - post_crop_extents[2]) / x_step)) - n_rows = abs(int(abs(post_crop_extents[1] - post_crop_extents[3]) / y_step)) - - for r in rasters: - r.close() - - return min_extents, n_cols, n_rows - -def _crop_opts(params: Dict) -> Tuple: - """ - Convenience function for getting crop options from parameters. - """ - from pyrate.core.prepifg_helper import CustomExts - - crop_opt = params[IFG_CROP_OPT] - if crop_opt == 3: - xfirst = params[IFG_XFIRST] - yfirst = params[IFG_YFIRST] - xlast = params[IFG_XLAST] - ylast = params[IFG_YLAST] - return crop_opt, CustomExts(xfirst, yfirst, xlast, ylast) - return crop_opt, None class ConfigException(Exception): """ diff --git a/pyrate/core/gdal_python.py b/pyrate/core/gdal_python.py index 4e9c10485..0d85f4d39 100644 --- a/pyrate/core/gdal_python.py +++ b/pyrate/core/gdal_python.py @@ -276,9 +276,9 @@ def _gdalwarp_width_and_height(max_x, max_y, min_x, min_y, geo_trans): def crop_resample_average( - input_tif, extents: Union[List, Tuple], new_res, output_file, thresh, - out_driver_type='GTiff', - match_pyrate=False, hdr=None, coherence_path=None, coherence_thresh=None): + input_tif, extents: Union[List, Tuple], new_res, output_file, thresh, hdr, out_driver_type='GTiff', + match_pyrate=False, coherence_path=None, coherence_thresh=None + ): """ Crop, resample, and average a geotiff image. @@ -325,11 +325,8 @@ def crop_resample_average( gt = dst_ds.GetGeoTransform() wkt = dst_ds.GetProjection() - # TEST HERE IF EXISTING FILE HAS PYRATE METADATA. IF NOT ADD HERE - if ifc.DATA_TYPE not in dst_ds.GetMetadata() and hdr is not None: - md = shared.collate_metadata(hdr) - else: - md = dst_ds.GetMetadata() + # insert metadata from the header + md = shared.collate_metadata(hdr) # update metadata for output @@ -342,15 +339,12 @@ def crop_resample_average( md.update({ifc.DATA_TYPE: ifc.COHERENCE}) elif (v == ifc.ORIG) and (coherence_path is None): md.update({ifc.DATA_TYPE: ifc.MULTILOOKED}) + elif v == ifc.COH: + md.update({ifc.DATA_TYPE: ifc.MULTILOOKED_COH}) elif v == ifc.DEM: md.update({ifc.DATA_TYPE: ifc.MLOOKED_DEM}) elif v == ifc.INCIDENCE: md.update({ifc.DATA_TYPE: ifc.MLOOKED_INC}) - elif (v == ifc.COHERENCE) and (coherence_path is None): - # during orbital fit multilooking - pass - elif v == ifc.MULTILOOKED: - pass else: raise TypeError(f'Data Type metadata {v} not recognised') diff --git a/pyrate/core/ifgconstants.py b/pyrate/core/ifgconstants.py index 9558ee583..a5ece2b91 100644 --- a/pyrate/core/ifgconstants.py +++ b/pyrate/core/ifgconstants.py @@ -41,7 +41,9 @@ PYRATE_ALPHA = 'CVD_ALPHA' COHERENCE = 'COHERENCE_MASKED_MULTILOOKED_IFG' MULTILOOKED = 'MULTILOOKED_IFG' +MULTILOOKED_COH = 'MULTILOOKED_COH' ORIG = 'ORIGINAL_IFG' +COH = 'ORIGINAL_COH' DEM = 'ORIGINAL_DEM' MLOOKED_DEM = 'MULTILOOKED_DEM' INCIDENCE = 'INCIDENCE_ANGLE_MAP' @@ -60,6 +62,15 @@ NAN_CONVERTED = 'CONVERTED' DATA_TYPE = 'DATA_TYPE' DATA_UNITS = 'DATA_UNITS' +INPUT_TYPE = 'INPUT_TYPE' + +PYRATE_REFPIX_X = 'REF_PIX_X' +PYRATE_REFPIX_Y = 'REF_PIX_Y' +PYRATE_REFPIX_LAT = 'REF_PIX_LAT' +PYRATE_REFPIX_LON = 'REF_PIX_LON' +PYRATE_MEAN_REF_AREA = 'REF_AREA_MEAN' +PYRATE_STDDEV_REF_AREA = 'REF_AREA_STDDEV' +SEQUENCE_POSITION = 'SEQUENCE_POSITION' DAYS_PER_YEAR = 365.25 # span of year, not a calendar year YEARS_PER_DAY = 1 / DAYS_PER_YEAR diff --git a/pyrate/core/orbital.py b/pyrate/core/orbital.py index 4c8e13287..4bf38479d 100644 --- a/pyrate/core/orbital.py +++ b/pyrate/core/orbital.py @@ -68,20 +68,13 @@ PART_CUBIC = cf.PART_CUBIC -def remove_orbital_error(ifgs: Iterable, params: dict, preread_ifgs=None) -> None: +def remove_orbital_error(ifgs: Iterable, params: dict, headers: List[dict], preread_ifgs=None) -> None: """ Wrapper function for PyRate orbital error removal functionality. NB: the ifg data is modified in situ, rather than create intermediate files. The network method assumes the given ifgs have already been reduced to a minimum spanning tree network. - - :param list ifgs: List of interferograms class objects - :param dict params: Dictionary containing configuration parameters - :param dict preread_ifgs: Dictionary containing information specifically - for MPI jobs (optional) - - :return: None - interferogram phase data is updated and saved to disk """ ifg_paths = [i.data_path for i in ifgs] if isinstance(ifgs[0], Ifg) else ifgs @@ -96,7 +89,9 @@ def remove_orbital_error(ifgs: Iterable, params: dict, preread_ifgs=None) -> Non xlooks=params[cf.ORBITAL_FIT_LOOKS_X], ylooks=params[cf.ORBITAL_FIT_LOOKS_Y], thresh=params[cf.NO_DATA_AVERAGING_THRESHOLD], - write_to_disc=False) + write_to_disc=False, + headers=headers + ) mlooked = [Ifg(m[1]) for m in mlooked_dataset] for m in mlooked: diff --git a/pyrate/core/prepifg_helper.py b/pyrate/core/prepifg_helper.py index 994244d19..7435607b1 100644 --- a/pyrate/core/prepifg_helper.py +++ b/pyrate/core/prepifg_helper.py @@ -24,13 +24,13 @@ from math import modf from numbers import Number from decimal import Decimal - +from typing import List, Tuple, Union from numpy import array, nan, isnan, nanmean, float32, zeros, sum as nsum from osgeo import gdal from pyrate.core.gdal_python import crop_resample_average from pyrate.core import config as cf -from pyrate.core.shared import output_tiff_filename, dem_or_ifg +from pyrate.core.shared import output_tiff_filename, dem_or_ifg, Ifg, DEM CustomExts = namedtuple('CustExtents', ['xfirst', 'yfirst', 'xlast', 'ylast']) @@ -45,7 +45,8 @@ GRID_TOL = 1e-6 -def get_analysis_extent(crop_opt, rasters, xlooks, ylooks, user_exts): +def get_analysis_extent(crop_opt: int, rasters: List[Union[Ifg, DEM]], xlooks: int, ylooks: int, + user_exts: Tuple[float, float, float, float]) -> Tuple[float, float, float, float]: """ Function checks prepifg parameters and returns extents/bounding box. @@ -144,7 +145,7 @@ def _get_extents(ifgs, crop_opt, user_exts=None): return extents -def prepare_ifg(raster_path, xlooks, ylooks, exts, thresh, crop_opt, write_to_disk=True, out_path=None, header=None, +def prepare_ifg(raster_path, xlooks, ylooks, exts, thresh, crop_opt, header, write_to_disk=True, out_path=None, coherence_path=None, coherence_thresh=None): """ Open, resample, crop and optionally save to disk an interferogram or DEM. @@ -202,7 +203,7 @@ def prepare_ifg(raster_path, xlooks, ylooks, exts, thresh, crop_opt, write_to_di # TODO: crop options 0 = no cropping? get rid of same size -def prepare_ifgs(raster_data_paths, crop_opt, xlooks, ylooks, thresh=0.5, user_exts=None, write_to_disc=True, +def prepare_ifgs(raster_data_paths, crop_opt, xlooks, ylooks, headers, thresh=0.5, user_exts=None, write_to_disc=True, out_path=None): """ Wrapper function to prepare a sequence of interferogram files for @@ -223,13 +224,13 @@ def prepare_ifgs(raster_data_paths, crop_opt, xlooks, ylooks, thresh=0.5, user_e :return: resampled_data: output cropped and resampled image :rtype: ndarray :return: out_ds: destination gdal dataset object - :rtype: gdal.Dataset + :rtype: List[gdal.Dataset] """ # use metadata check to check whether it's a dem or ifg rasters = [dem_or_ifg(r) for r in raster_data_paths] exts = get_analysis_extent(crop_opt, rasters, xlooks, ylooks, user_exts) - - return [prepare_ifg(d, xlooks, ylooks, exts, thresh, crop_opt, write_to_disc, out_path) for d in raster_data_paths] + return [prepare_ifg(d, xlooks, ylooks, exts, thresh, crop_opt, h, write_to_disc, out_path) for d, h + in zip(raster_data_paths, headers)] # TODO: Not being used. Remove in future? @@ -267,7 +268,7 @@ def _resample(data, xscale, yscale, thresh): return dest -def _min_bounds(ifgs): +def _min_bounds(ifgs: List[Ifg]) -> Tuple[float, float, float, float]: """ Returns bounds for overlapping area of the given interferograms. """ @@ -279,7 +280,7 @@ def _min_bounds(ifgs): return xmin, ymin, xmax, ymax -def _max_bounds(ifgs): +def _max_bounds(ifgs: List[Ifg]) -> Tuple[float, float, float, float]: """ Returns bounds for the total area covered by the given interferograms. """ @@ -291,7 +292,7 @@ def _max_bounds(ifgs): return xmin, ymin, xmax, ymax -def _get_same_bounds(ifgs): +def _get_same_bounds(ifgs: List[Ifg]) -> Tuple[float, float, float, float]: """ Check and return bounding box for ALREADY_SAME_SIZE option. """ diff --git a/pyrate/core/ref_phs_est.py b/pyrate/core/ref_phs_est.py index bbf94225f..622bcfeaa 100644 --- a/pyrate/core/ref_phs_est.py +++ b/pyrate/core/ref_phs_est.py @@ -29,6 +29,7 @@ MASTER_PROCESS = 0 + def est_ref_phase_method2(ifg_paths, params, refpx, refpy): """ Reference phase estimation using method 2. Reference phase is the @@ -61,8 +62,7 @@ def _inner(ifg_paths): if params[cf.PARALLEL]: ref_phs = Parallel(n_jobs=params[cf.PROCESSES], verbose=joblib_log_level(cf.LOG_LEVEL))( - delayed(_est_ref_phs_method2)(p, half_chip_size, - refpx, refpy, thresh) + delayed(_est_ref_phs_method2)(p, half_chip_size, refpx, refpy, thresh) for p in phase_data) for n, ifg in enumerate(ifgs): @@ -70,20 +70,20 @@ def _inner(ifg_paths): else: ref_phs = np.zeros(len(ifgs)) for n, ifg in enumerate(ifgs): - ref_phs[n] = \ - _est_ref_phs_method2(phase_data[n], half_chip_size, - refpx, refpy, thresh) + ref_phs[n] = _est_ref_phs_method2(phase_data[n], half_chip_size, refpx, refpy, thresh) ifg.phase_data -= ref_phs[n] for ifg in ifgs: _update_phase_metadata(ifg) ifg.close() + return ref_phs process_ifgs_paths = mpiops.array_split(ifg_paths) ref_phs = _inner(process_ifgs_paths) return ref_phs + def _est_ref_phs_method2(phase_data, half_chip_size, refpx, refpy, thresh): """ Convenience function for ref phs estimate method 2 parallelisation @@ -144,10 +144,9 @@ def _inner(proc_ifgs, phase_data_sum): comp = np.ravel(comp, order='F') if params[cf.PARALLEL]: - phase_data = [i.phase_data for i in proc_ifgs] log.info("Calculating ref phase using multiprocessing") ref_phs = Parallel(n_jobs=params[cf.PROCESSES], verbose=joblib_log_level(cf.LOG_LEVEL))( - delayed(_est_ref_phs_method1)(p, comp) for p in phase_data + delayed(_est_ref_phs_method1)(p.phase_data, comp) for p in proc_ifgs ) for n, ifg in enumerate(proc_ifgs): ifg.phase_data -= ref_phs[n] @@ -183,6 +182,7 @@ def _est_ref_phs_method1(phase_data, comp): def _update_phase_metadata(ifg): ifg.meta_data[ifc.PYRATE_REF_PHASE] = ifc.REF_PHASE_REMOVED ifg.write_modified_phase() + log.debug(f"Reference phase corrected for {ifg.data_path}") class ReferencePhaseError(Exception): diff --git a/pyrate/core/refpixel.py b/pyrate/core/refpixel.py index 8539d33a6..d25fbafba 100644 --- a/pyrate/core/refpixel.py +++ b/pyrate/core/refpixel.py @@ -25,10 +25,103 @@ from numpy import isnan, std, mean, sum as nsum from joblib import Parallel, delayed +from osgeo import gdal + import pyrate.core.config as cf +from pyrate.core import ifgconstants as ifc +from pyrate.core import mpiops from pyrate.core.shared import Ifg from pyrate.core.shared import joblib_log_level from pyrate.core.logger import pyratelogger as log +from pyrate.core import prepifg_helper + + +def update_refpix_metadata(ifg_paths, refx, refy, transform, params): + """ + Function that adds metadata about the chosen reference pixel to each interferogram. + """ + + pyrate_refpix_lon, pyrate_refpix_lat = mpiops.run_once(convert_pixel_value_to_geographic_coordinate, + refx, refy, transform) + + process_ifgs_paths = mpiops.array_split(ifg_paths) + + for ifg_file in process_ifgs_paths: + log.debug("Updating metadata for: "+ifg_file) + ifg = Ifg(ifg_file) + log.debug("Open dataset") + ifg.open(readonly=True) + log.debug("Set no data value") + ifg.nodata_value = params["noDataValue"] + log.debug("Update no data values in dataset") + ifg.convert_to_nans() + log.debug("Convert mm") + ifg.convert_to_mm() + half_patch_size = params["refchipsize"] // 2 + x, y = refx, refy + log.debug("Extract reference pixel windows") + data = ifg.phase_data[y - half_patch_size: y + half_patch_size + 1, + x - half_patch_size: x + half_patch_size + 1] + log.debug("Calculate standard deviation for reference window") + stddev_ref_area = np.nanstd(data) + log.debug("Calculate mean for reference window") + mean_ref_area = np.nanmean(data) + ifg.add_metadata(**{ + ifc.PYRATE_REFPIX_X: str(refx), + ifc.PYRATE_REFPIX_Y: str(refy), + ifc.PYRATE_REFPIX_LAT: str(pyrate_refpix_lat), + ifc.PYRATE_REFPIX_LON: str(pyrate_refpix_lon), + ifc.PYRATE_MEAN_REF_AREA: str(mean_ref_area), + ifc.PYRATE_STDDEV_REF_AREA: str(stddev_ref_area) + }) + + ifg.close() + + +def convert_pixel_value_to_geographic_coordinate(refx, refy, transform): + """ + Converts a pixel coordinate to a latitude/longitude coordinate given the + geotransform of the image. + Args: + refx: The pixel x coordinate. + refy: The pixel ye coordinate. + transform: The geotransform array of the image. + Returns: + Tuple of lon, lat geographic coordinate. + """ + + xOrigin = transform[0] + yOrigin = transform[3] + pixelWidth = transform[1] + pixelHeight = -transform[5] + + lon = refx*pixelWidth + xOrigin + lat = yOrigin - refy*pixelHeight + + return lon, lat + + +def convert_geographic_coordinate_to_pixel_value(lon, lat, transform): + """ + Converts a latitude/longitude coordinate to a pixel coordinate given the + geotransform of the image. + Args: + lon: Pixel longitude. + lat: Pixel latitude. + transform: The geotransform array of the image. + Returns: + Tuple of refx, refy pixel coordinates. + """ + + xOrigin = transform[0] + yOrigin = transform[3] + pixelWidth = transform[1] + pixelHeight = -transform[5] + + refx = round((lon - xOrigin) / pixelWidth) + refy = round((yOrigin - lat) / pixelHeight) + + return refx, refy # TODO: move error checking to config step (for fail fast) @@ -66,8 +159,8 @@ def ref_pixel(ifgs, params): mean_sds.append(_ref_pixel_multi(g, half_patch_size, phase_data, thresh, params)) refxy = find_min_mean(mean_sds, grid) - if isinstance(refxy, ValueError): - raise ValueError('Refpixel calculation not possible!') + if isinstance(refxy, RefPixelError): + raise RefPixelError('Refpixel calculation not possible!') refy, refx = refxy @@ -92,7 +185,7 @@ def find_min_mean(mean_sds, grid): try: refp_index = np.nanargmin(mean_sds) return grid[refp_index] - except ValueError as v: + except RefPixelError as v: log.error(v) return v @@ -250,7 +343,7 @@ def _validate_chipsize(chipsize, head): if chipsize < 3 or chipsize > head.ncols or (chipsize % 2 == 0): msg = "Chipsize setting must be >=3 and at least <= grid width" - raise ValueError(msg) + raise RefPixelError(msg) log.debug('Chipsize validation successful') @@ -262,7 +355,7 @@ def _validate_minimum_fraction(min_frac): raise cf.ConfigException('Minimum fraction is None') if min_frac < 0.0 or min_frac > 1.0: - raise ValueError("Minimum fraction setting must be >= 0.0 and <= 1.0 ") + raise RefPixelError("Minimum fraction setting must be >= 0.0 and <= 1.0 ") def _validate_search_win(refnx, refny, chipsize, head): @@ -275,7 +368,7 @@ def _validate_search_win(refnx, refny, chipsize, head): max_width = (head.ncols - (chipsize-1)) if refnx < 1 or refnx > max_width: msg = "Invalid refnx setting, must be > 0 and <= %s" - raise ValueError(msg % max_width) + raise RefPixelError(msg % max_width) if refny is None: raise cf.ConfigException('refny is None') @@ -283,10 +376,33 @@ def _validate_search_win(refnx, refny, chipsize, head): max_rows = (head.nrows - (chipsize-1)) if refny < 1 or refny > max_rows: msg = "Invalid refny setting, must be > 0 and <= %s" - raise ValueError(msg % max_rows) + raise RefPixelError(msg % max_rows) + + +def validate_supplied_lat_lon(params: dict) -> None: + """ + Function to validate that the user supplied lat/lon values sit within image bounds + """ + lon, lat = params[cf.REFX], params[cf.REFY] + if lon == -1 or lat == -1: + return + xmin, ymin, xmax, ymax = prepifg_helper.get_analysis_extent( + crop_opt=params[cf.IFG_CROP_OPT], + rasters=[prepifg_helper.dem_or_ifg(p.sampled_path) for p in params[cf.INTERFEROGRAM_FILES]], + xlooks=params[cf.IFG_LKSX], ylooks=params[cf.IFG_LKSY], + user_exts=(params[cf.IFG_XFIRST], params[cf.IFG_YFIRST], params[cf.IFG_XLAST], params[cf.IFG_YLAST]) + ) + msg = "Supplied {} value is outside the bounds of the interferogram data" + lat_lon_txt = '' + if (lon < xmin) or (lon > xmax): + lat_lon_txt += 'longitude' + if (lat < ymin) or (lat > ymax): + lat_lon_txt += ' and latitude' if lat_lon_txt else 'latitude' + if lat_lon_txt: + raise RefPixelError(msg.format(lat_lon_txt)) class RefPixelError(Exception): - ''' + """ Generic exception for reference pixel errors. - ''' + """ diff --git a/pyrate/core/roipac.py b/pyrate/core/roipac.py index a6dfcc193..55c4861c0 100644 --- a/pyrate/core/roipac.py +++ b/pyrate/core/roipac.py @@ -211,9 +211,9 @@ def roipac_header(file_path, params): projection = parse_header(rsc_file)[ifc.PYRATE_DATUM] else: raise RoipacException('No DEM resource/header file is provided') - if file_path.endswith('dem.tif'): + if file_path.endswith('_dem.tif'): header_file = os.path.join(params[cf.DEM_HEADER_FILE]) - elif file_path.endswith('unw.tif'): + elif file_path.endswith('unw_ifg.tif') or file_path.endswith('unw.tif'): # TODO: improve this interferogram_epoches = extract_epochs_from_filename(p.name) for header_path in params[cf.HEADER_FILE_PATHS]: diff --git a/pyrate/core/shared.py b/pyrate/core/shared.py index be78fd870..642b21165 100644 --- a/pyrate/core/shared.py +++ b/pyrate/core/shared.py @@ -27,9 +27,11 @@ from math import floor import os from os.path import basename, join +from pathlib import Path import struct from datetime import date from itertools import product +from enum import Enum import numpy as np from numpy import where, nan, isnan, sum as nsum, isclose import pyproj @@ -62,6 +64,13 @@ GDAL_Y_FIRST = 3 +class InputTypes(Enum): + IFG = '_ifg' + COH = '_coh' + DEM = '_dem' + HEADER = '_header' + + def joblib_log_level(level: str) -> int: """ Convert python log level to joblib int verbosity. @@ -266,12 +275,13 @@ class Ifg(RasterBase): interferometric phase raster band data and related data. """ # pylint: disable=too-many-instance-attributes - def __init__(self, path): + def __init__(self, path: Union[str, Path]): """ Interferogram constructor, for 2-band Ifg raster datasets. :param str path: Path to interferogram file """ + path = path.as_posix() if isinstance(path, Path) else path RasterBase.__init__(self, path) self._phase_band = None self._phase_data = None @@ -469,6 +479,14 @@ def write_modified_phase(self, data=None): self.dataset.SetMetadataItem(k, v) self.dataset.FlushCache() + def add_metadata(self, **kwargs): + if (not self.is_open) or self.is_read_only: + raise IOError("Ifg not open or readonly. Cannot write!") + + for k, v in kwargs.items(): + self.dataset.SetMetadataItem(k, v) + self.dataset.FlushCache() # write to disc + class IfgPart(object): """ @@ -682,7 +700,16 @@ def _is_interferogram(hdr): """ Convenience function to determine if file is interferogram """ - return ifc.PYRATE_WAVELENGTH_METRES in hdr + return (ifc.PYRATE_WAVELENGTH_METRES in hdr) and \ + (hdr[ifc.INPUT_TYPE] == InputTypes.IFG if ifc.INPUT_TYPE in hdr else True) + + +def _is_coherence(hdr): + """ + Convenience function to determine if file is interferogram + """ + return (ifc.PYRATE_WAVELENGTH_METRES in hdr) and \ + (hdr[ifc.INPUT_TYPE] == InputTypes.COH if ifc.INPUT_TYPE in hdr else False) def _is_incidence(hdr): @@ -728,7 +755,7 @@ def write_fullres_geotiff(header, data_path, dest, nodata): raise GeotiffException(msg) wkt = srs.ExportToWkt() - dtype = 'float32' if (_is_interferogram(header) or _is_incidence(header)) else 'int16' + dtype = 'float32' if (_is_interferogram(header) or _is_incidence(header) or _is_coherence(header)) else 'int16' # get subset of metadata relevant to PyRate md = collate_metadata(header) @@ -797,19 +824,27 @@ def collate_metadata(header): :return: dict of relevant metadata for PyRate """ md = dict() - if _is_interferogram(header): + + def __common_ifg_coh_update(header, md): for k in [ifc.PYRATE_WAVELENGTH_METRES, ifc.PYRATE_TIME_SPAN, ifc.PYRATE_INSAR_PROCESSOR, ifc.MASTER_DATE, ifc.SLAVE_DATE, - ifc.DATA_UNITS, ifc.DATA_TYPE]: + ifc.DATA_UNITS]: md.update({k: str(header[k])}) if header[ifc.PYRATE_INSAR_PROCESSOR] == GAMMA: for k in [ifc.MASTER_TIME, ifc.SLAVE_TIME, ifc.PYRATE_INCIDENCE_DEGREES]: md.update({k: str(header[k])}) + + if _is_coherence(header): + __common_ifg_coh_update(header, md) + md.update({ifc.DATA_TYPE: ifc.COH}) + elif _is_interferogram(header): + __common_ifg_coh_update(header, md) + md.update({ifc.DATA_TYPE: ifc.ORIG}) elif _is_incidence(header): - md.update({ifc.DATA_TYPE:ifc.INCIDENCE}) - else: # must be dem - md.update({ifc.DATA_TYPE:ifc.DEM}) + md.update({ifc.DATA_TYPE: ifc.INCIDENCE}) + else: # must be dem + md.update({ifc.DATA_TYPE: ifc.DEM}) return md @@ -915,21 +950,13 @@ def write_output_geotiff(md, gt, wkt, data, dest, nodata): # set other metadata ds.SetMetadataItem('DATA_TYPE', str(md['DATA_TYPE'])) + # sequence position for time series products - if "SEQUENCE_POSITION" in md: - ds.SetMetadataItem("SEQUENCE_POSITION", str(md["SEQUENCE_POSITION"])) - if "PYRATE_REFPIX_LAT" in md: - ds.SetMetadataItem("PYRATE_REFPIX_LAT", str(md["PYRATE_REFPIX_LAT"])) - if "PYRATE_REFPIX_LON" in md: - ds.SetMetadataItem("PYRATE_REFPIX_LON", str(md["PYRATE_REFPIX_LON"])) - if "PYRATE_REFPIX_X" in md: - ds.SetMetadataItem("PYRATE_REFPIX_X", str(md["PYRATE_REFPIX_X"])) - if "PYRATE_REFPIX_Y" in md: - ds.SetMetadataItem("PYRATE_REFPIX_Y", str(md["PYRATE_REFPIX_Y"])) - if "PYRATE_MEAN_REF_AREA" in md: - ds.SetMetadataItem("PYRATE_MEAN_REF_AREA", str(md["PYRATE_MEAN_REF_AREA"])) - if "STANDARD_DEVIATION_REF_AREA" in md: - ds.SetMetadataItem("PYRATE_STANDARD_DEVIATION_REF_AREA", str(md["PYRATE_STANDARD_DEVIATION_REF_AREA"])) + + for k in [ifc.SEQUENCE_POSITION, ifc.PYRATE_REFPIX_X, ifc.PYRATE_REFPIX_Y, ifc.PYRATE_REFPIX_LAT, + ifc.PYRATE_REFPIX_LON, ifc.PYRATE_MEAN_REF_AREA, ifc.PYRATE_STDDEV_REF_AREA]: + if k in md: + ds.SetMetadataItem(k, str(md[k])) # write data to geotiff band = ds.GetRasterBand(1) @@ -1285,13 +1312,13 @@ def extract_epochs_from_filename(filename_with_epochs: str) -> List[str]: def mpi_vs_multiprocess_logging(step, params): if mpiops.size > 1: # Over-ride input options if this is an MPI job - log.info(f"Running {step} step using mpi processing. Disabling parallel processing.") + log.info(f"Running {step} step using MPI processing. Disabling parallel processing.") params[cf.PARALLEL] = 0 else: if params[cf.PARALLEL] == 1: - log.info(f"Running {step} using {params[cf.PROCESSES]} processes") + log.info(f"Running {step} step in parallel using {params[cf.PROCESSES]} processes") else: - log.info(f"Running {step} serially") + log.info(f"Running {step} step in serial") def dem_or_ifg(data_path): @@ -1308,4 +1335,4 @@ def dem_or_ifg(data_path): if ifc.MASTER_DATE in md: # ifg return Ifg(data_path) else: - return DEM(data_path) \ No newline at end of file + return DEM(data_path) diff --git a/pyrate/core/stack.py b/pyrate/core/stack.py index 5b5e33709..3037ecc49 100644 --- a/pyrate/core/stack.py +++ b/pyrate/core/stack.py @@ -21,7 +21,7 @@ import itertools from scipy.linalg import solve, cholesky, qr, inv -from numpy import nan, isnan, sqrt, diag, delete, array, float32 +from numpy import nan, isnan, sqrt, diag, delete, array, float32, size import numpy as np from joblib import Parallel, delayed from pyrate.core import config as cf @@ -29,12 +29,13 @@ from pyrate.core.logger import pyratelogger as log -def stack_rate(ifgs, params, vcmt, mst=None): +def stack_rate_array(ifgs, params, vcmt, mst=None): """ - Pixel-by-pixel rate (velocity) estimation using iterative - weighted least-squares stacking method. + This function loops over all interferogram pixels in a 3-dimensional array and estimates + the pixel rate (velocity) by applying the iterative weighted least-squares stacking + algorithm 'pyrate.core.stack.stack_rate_pixel'. - :param Ifg.object ifgs: Sequence of interferogram objects from which to extract observations + :param Ifg.object ifgs: Sequence of objects containing the interferometric observations :param dict params: Configuration parameters :param ndarray vcmt: Derived positive definite temporal variance covariance matrix :param ndarray mst: Pixel-wise matrix describing the minimum spanning tree network @@ -43,17 +44,17 @@ def stack_rate(ifgs, params, vcmt, mst=None): :rtype: ndarray :return: error: Standard deviation of the rate map :rtype: ndarray - :return: samples: Number of observations used in rate calculation per pixel + :return: samples: Number of observations used in rate calculation for each pixel :rtype: ndarray """ - maxsig, nsig, pthresh, cols, error, mst, obs, parallel, _, rate, rows, samples, span = _stack_setup(ifgs, mst, params) + nsig, pthresh, cols, error, mst, obs, parallel, _, rate, rows, samples, span = _stack_setup(ifgs, mst, params) # pixel-by-pixel calculation. # nested loops to loop over the 2 image dimensions if parallel: log.info('Calculating stack rate in parallel') res = Parallel(n_jobs=params[cf.PROCESSES], verbose=joblib_log_level(cf.LOG_LEVEL))( - delayed(_stack_rate_by_pixel)(r, c, mst, nsig, obs, pthresh, span, vcmt) for r, c in itertools.product(range(rows), range(cols)) + delayed(stack_rate_pixel)(obs[:, r, c], mst[:, r, c], vcmt, span, nsig, pthresh) for r, c in itertools.product(range(rows), range(cols)) ) res = np.array(res) @@ -64,61 +65,70 @@ def stack_rate(ifgs, params, vcmt, mst=None): log.info('Calculating stack rate in serial') for i in range(rows): for j in range(cols): - rate[i, j], error[i, j], samples[i, j] = _stack_rate_by_pixel(i, j, mst, nsig, obs, pthresh, span, vcmt) + rate[i, j], error[i, j], samples[i, j] = stack_rate_pixel(obs[:, i, j], mst[:, i, j], vcmt, span, nsig, pthresh) + + return rate, error, samples + - # overwrite the data whose error is larger than the - # maximum sigma user threshold +def mask_rate(rate, error, maxsig): + """ + Function to mask pixels in the rate and error arrays when the error + is greater than the error threshold 'maxsig'. + + :param ndarray rate: array of pixel rates derived by stacking + :param ndarray error: array of errors for the pixel rates + :param int maxsig: error threshold for masking (in millimetres). + + :return: rate: Masked rate (velocity) map + :rtype: ndarray + :return: error: Masked error (standard deviation) map + :rtype: ndarray + """ + log.info('Masking stack rate and error maps where sigma is greater than {} millimetres'.format(maxsig)) + # initialise mask array with existing NaNs mask = ~isnan(error) + # original Nan count + orig = np.count_nonzero(mask) + # determine where error is larger than the maximum sigma threshold mask[mask] &= error[mask] > maxsig + # replace values with NaNs rate[mask] = nan error[mask] = nan - # samples[mask] = nan # should we also mask the samples? + # calculate percentage of masked pixels + nummasked = int(np.count_nonzero(mask)/orig*100) + log.info('Percentage of pixels masked = {}%'.format(nummasked)) - return rate, error, samples + return rate, error -def _stack_setup(ifgs, mst, params): +def stack_rate_pixel(obs, mst, vcmt, span, nsig, pthresh): """ - Convenience function for stack rate setup - """ - # MULTIPROCESSING parameters - parallel = params[cf.PARALLEL] - processes = params[cf.PROCESSES] - # stack rate parameters from config file - # n-sigma ratio used to threshold 'model minus observation' residuals - nsig = params[cf.LR_NSIG] - # Threshold for maximum allowable standard error - maxsig = params[cf.LR_MAXSIG] - # Pixel threshold; minimum number of coherent observations for a pixel - pthresh = params[cf.LR_PTHRESH] - rows, cols = ifgs[0].phase_data.shape - # make 3D block of observations - obs = array([np.where(isnan(x.phase_data), 0, x.phase_data) for x in ifgs]) - span = array([[x.time_span for x in ifgs]]) - # Update MST in case additional NaNs generated by APS filtering - if mst is None: # dummy mst if none is passed in - mst = ~isnan(obs) - else: - mst[isnan(obs)] = 0 - - # preallocate empty arrays. No need to preallocation NaNs with new code - error = np.empty([rows, cols], dtype=float32) - rate = np.empty([rows, cols], dtype=float32) - samples = np.empty([rows, cols], dtype=np.float32) - return maxsig, nsig, pthresh, cols, error, mst, obs, parallel, processes, rate, rows, samples, span - + Algorithm to estimate the rate (velocity) for a single pixel using iterative + weighted least-squares stacking method. -def _stack_rate_by_pixel(row, col, mst, nsig, obs, pthresh, span, vcmt): - """helper function for computing stack rate for one pixel""" + :param ndarray obs: Vector of interferometric phase observations for the pixel + :param ndarray mst: Vector describing the minimum spanning tree network for the pixel + :param ndarray vcmt: Derived positive definite temporal variance covariance matrix + :param ndarray span: Vector of interferometric time spans + :param int nsig: Threshold for iterative removal of interferometric observations + :param int pthresh: Threshold for minimum number of observations for the pixel + + :return: rate: Estimated rate (velocity) for the pixel + :rtype: float64 + :return: error: Standard deviation of the observations with respect to the estimated rate + :rtype: float64 + :return: samples: Number of observations used in the rate estimation for the pixel + :rtype: int + """ - # find the indices of independent ifgs for given pixel from MST - ind = np.nonzero(mst[:, row, col])[0] # only True's in mst are chosen + # find the indices of independent ifgs from MST + ind = np.nonzero(mst)[0] # only True's in mst are chosen # iterative loop to calculate 'robust' velocity for pixel default_no_samples = len(ind) while len(ind) >= pthresh: - # make vector of selected ifg observations - ifgv = obs[ind, row, col] + # select ifg observations + ifgv = obs[ind] # form design matrix from appropriate ifg time spans B = span[:, ind] @@ -165,3 +175,32 @@ def _stack_rate_by_pixel(row, col, mst, nsig, obs, pthresh, span, vcmt): return v[0], err[0], ifgv.shape[0] # dummy return for no change return np.nan, np.nan, default_no_samples + + +def _stack_setup(ifgs, mst, params): + """ + Convenience function for stack rate setup + """ + # MULTIPROCESSING parameters + parallel = params[cf.PARALLEL] + processes = params[cf.PROCESSES] + # stack rate parameters from config file + # n-sigma ratio used to threshold 'model minus observation' residuals + nsig = params[cf.LR_NSIG] + # Pixel threshold; minimum number of observations for a pixel + pthresh = params[cf.LR_PTHRESH] + rows, cols = ifgs[0].phase_data.shape + # make 3D block of observations + obs = array([np.where(isnan(x.phase_data), 0, x.phase_data) for x in ifgs]) + span = array([[x.time_span for x in ifgs]]) + # Update MST in case additional NaNs generated by APS filtering + if mst is None: # dummy mst if none is passed in + mst = ~isnan(obs) + else: + mst[isnan(obs)] = 0 + + # preallocate empty arrays (no need to preallocate NaNs) + error = np.empty([rows, cols], dtype=float32) + rate = np.empty([rows, cols], dtype=float32) + samples = np.empty([rows, cols], dtype=np.float32) + return nsig, pthresh, cols, error, mst, obs, parallel, processes, rate, rows, samples, span diff --git a/pyrate/default_parameters.py b/pyrate/default_parameters.py index 5a7b3afb8..5f8220dd7 100644 --- a/pyrate/default_parameters.py +++ b/pyrate/default_parameters.py @@ -428,5 +428,14 @@ "MaxValue": 1000, "PossibleValues": None, "Required": False - } + }, + "savenpy": { + "DataType": int, + "DefaultValue": 0, + "MinValue": 0, + "MaxValue": 1, + "PossibleValues": [1, 0], + "Required": False + }, + } diff --git a/pyrate/merge.py b/pyrate/merge.py index 1a8068648..6f969aa0f 100644 --- a/pyrate/merge.py +++ b/pyrate/merge.py @@ -15,29 +15,22 @@ # limitations under the License. """ This Python module does post-processing steps to assemble the -rate and time series outputs and save as geotiff files +stack rate and time series outputs and save as geotiff files """ -import os -from os.path import join -import pickle as cp +from os.path import join, isfile +import pickle import numpy as np from osgeo import gdal import subprocess from pathlib import Path +from math import floor -from pyrate.core import shared, ifgconstants as ifc, mpiops, config as cf -from pyrate.core.shared import PrereadIfg -from pyrate.constants import REF_COLOR_MAP_PATH -from pyrate.core.config import OBS_DIR, OUT_DIR, ConfigException +from pyrate.core import shared, stack, ifgconstants as ifc, mpiops, config as cf from pyrate.core.logger import pyratelogger as log gdal.SetCacheMax(64) -# Constants -MASTER_PROCESS = 0 - - def main(params): """ PyRate merge main function. Assembles product tiles in to @@ -45,244 +38,265 @@ def main(params): """ # setup paths rows, cols = params["rows"], params["cols"] - _merge_stack(rows, cols, params) + mpiops.run_once(_merge_stack, rows, cols, params) + mpiops.run_once(_create_png_from_tif, params[cf.OUT_DIR]) if params[cf.TIME_SERIES_CAL]: _merge_timeseries(rows, cols, params) - #mpiops.run_once(_delete_tsincr_files, params) + # mpiops.run_once(_delete_tsincr_files, params) - log.info('Creating quicklook images.') - mpiops.run_once(create_png_from_tif, params[cf.OUT_DIR]) - log.debug('Finished creating quicklook images.') +def _merge_stack(rows, cols, params): + """ + Merge stacking outputs + """ + shape, tiles, ifgs_dict = _merge_setup(rows, cols, params) -def create_png_from_tif(output_folder_path): + log.info('Merging and writing Stack Rate product geotiffs') - # open raster and choose band to find min, max - raster_path = os.path.join(output_folder_path, "stack_rate.tif") + # read and assemble tile outputs + rate = assemble_tiles(shape, params[cf.TMPDIR], tiles, out_type='stack_rate') + error = assemble_tiles(shape, params[cf.TMPDIR], tiles, out_type='stack_error') + samples = assemble_tiles(shape, params[cf.TMPDIR], tiles, out_type='stack_samples') - if not os.path.isfile(raster_path): - raise Exception("stack_rate.tif file not found at: "+raster_path) - gtif = gdal.Open(raster_path) - srcband = gtif.GetRasterBand(1) + # mask pixels according to threshold + if params[cf.LR_MAXSIG] > 0: + rate, error = stack.mask_rate(rate, error, params[cf.LR_MAXSIG]) + else: + log.info('Skipping stack product masking (maxsig = 0)') + + # save geotiff and numpy array files + for out, ot in zip([rate, error, samples], ['stack_rate', 'stack_error', 'stack_samples']): + _save_merged_files(ifgs_dict, params[cf.OUT_DIR], out, ot, savenpy=params["savenpy"]) + + +def _merge_timeseries(rows, cols, params): + """ + Merge time series output + """ + log.info("Merging timeseries output") + shape, tiles, ifgs_dict = _merge_setup(rows, cols, params) + + # load the first tsincr file to determine the number of time series tifs + tsincr_file = join(params[cf.TMPDIR], 'tsincr_0.npy') + tsincr = np.load(file=tsincr_file) + # pylint: disable=no-member + no_ts_tifs = tsincr.shape[2] + + # create 2 x no_ts_tifs as we are splitting tsincr and tscuml to all processes. + process_tifs = mpiops.array_split(range(2 * no_ts_tifs)) + + # depending on nvelpar, this will not fit in memory + # e.g. nvelpar=100, nrows=10000, ncols=10000, 32bit floats need 40GB memory + # 32 * 100 * 10000 * 10000 / 8 bytes = 4e10 bytes = 40 GB + # the double for loop helps us overcome the memory limit + log.info('Process {} writing {} timeseries tifs of ' + 'total {}'.format(mpiops.rank, len(process_tifs), no_ts_tifs * 2)) + for i in process_tifs: + if i < no_ts_tifs: + tscum_g = assemble_tiles(shape, params[cf.TMPDIR], tiles, out_type='tscuml', index=i) + _save_merged_files(ifgs_dict, params[cf.OUT_DIR], tscum_g, out_type='tscuml', index=i, + savenpy=params["savenpy"]) + else: + i %= no_ts_tifs + tsincr_g = assemble_tiles(shape, params[cf.TMPDIR], tiles, out_type='tsincr', index=i) + _save_merged_files(ifgs_dict, params[cf.OUT_DIR], tsincr_g, out_type='tsincr', index=i, + savenpy=params["savenpy"]) + mpiops.comm.barrier() + log.debug('Process {} finished writing {} timeseries tifs of ' + 'total {}'.format(mpiops.rank, len(process_tifs), no_ts_tifs * 2)) + + +def _create_png_from_tif(output_folder_path): + """ + Wrapper for rate and error png/kml generation + """ + create_png_and_kml_from_tif(output_folder_path, output_type='rate') + create_png_and_kml_from_tif(output_folder_path, output_type='error') + + +def create_png_and_kml_from_tif(output_folder_path: str, output_type: str) -> None: + """ + Function to create a preview PNG format image from a geotiff, and a KML file + """ + log.info(f'Creating quicklook image for stack_{output_type}') + # open raster and choose band to find min, max + raster_path = join(output_folder_path, f"stack_{output_type}.tif") + if not isfile(raster_path): + raise Exception(f"stack_{output_type}.tif file not found at: " + raster_path) + gtif = gdal.Open(raster_path) + # find bounds of image west, north, east, south = "", "", "", "" for line in gdal.Info(gtif).split('\n'): if "Upper Left" in line: west, north = line.split(")")[0].split("(")[1].split(",") if "Lower Right" in line: east, south = line.split(")")[0].split("(")[1].split(",") - - kml_file_path = os.path.join(output_folder_path, "stack_rate.kml") - kml_file_content = """ + # write KML file + kml_file_path = join(output_folder_path, f"stack_{output_type}.kml") + kml_file_content = f""" - stack_rate.kml + stack_{output_type}.kml - stack_rate.png + stack_{output_type}.png - stack_rate.png + stack_{output_type}.png - """+north+""" - """+south+""" - """+east+""" - """+west+""" + """ + north + """ + """ + south + """ + """ + east + """ + """ + west + """ """ - with open(kml_file_path, "w") as f: f.write(kml_file_content) - + # Get raster statistics - minimum, maximum, mean, stddev = srcband.GetStatistics(True, True) - maximum = max(abs(minimum), abs(maximum)) - minimum = -1 * maximum - step = (maximum - minimum) / 256.0 - - del gtif # manually close raster - - # read color map from utilities and write it to the output folder - - with open(REF_COLOR_MAP_PATH, "r") as f: - color_map_list = [] - for line in f.readlines(): - color_map_list.append(line.strip().split(" ")) - - no_of_data_value = len(np.arange(minimum, maximum, step)) - for i, no in enumerate(np.arange(minimum, maximum, step)): - color_map_list[i+1][0] = str(no) - - color_map_path = os.path.join(output_folder_path, "colourmap.txt") + srcband = gtif.GetRasterBand(1) + minimum, maximum, _, _ = srcband.GetStatistics(True, True) + del gtif # close geotiff (used to calculate statistics) + # steps used for the colourmap, must be even (currently hard-coded to 254 resulting in 255 values) + no_of_steps = 254 + # slightly different code required for rate map and rate error map + if output_type == 'rate': + # minimum value might be negative + maximum = max(abs(minimum), abs(maximum)) + minimum = -1 * maximum + # colours: blue -> white -> red (white==0) + # note that an extra value will be added for zero (i.e. white: 255 255 255) + # generate a colourmap for odd number of values (currently hard-coded to 255) + mid = int(no_of_steps * 0.5) + # allocate RGB values to three numpy arrays r, g, b + r = np.arange(0, mid) / mid + g = r + r = np.concatenate((r, np.ones(mid + 1))) + g = np.concatenate((g, np.array([1]), np.flipud(g))) + b = np.flipud(r) + # change direction of colours (blue: positve, red: negative) + r = np.flipud(r) * 255 + g = np.flipud(g) * 255 + b = np.flipud(b) * 255 + if output_type == 'error': + # colours: white -> red (minimum error -> maximum error + # allocate RGB values to three numpy arrays r, g, b + r = np.ones(no_of_steps+1)*255 + g = np.arange(0,no_of_steps+1)/(no_of_steps) + g = np.flipud(g)*255 + b = g + # generate the colourmap file in the output folder + color_map_path = join(output_folder_path, f"colourmap_{output_type}.txt") + log.info( + 'Saving colour map to file {}; min/max values: {:.2f}/{:.2f}'.format(color_map_path, minimum, + maximum)) with open(color_map_path, "w") as f: - for i in range(no_of_data_value): - f.write(' '.join(color_map_list[i]) + "\n") - - input_tif_path = os.path.join(output_folder_path, "stack_rate.tif") - output_png_path = os.path.join(output_folder_path, "stack_rate.png") + f.write("nan 0 0 0 0\n") + for i, value in enumerate(np.linspace(minimum, maximum, no_of_steps+1)): + f.write("%f %f %f %f 255\n" % (value, r[i], g[i], b[i])) + input_tif_path = join(output_folder_path, f"stack_{output_type}.tif") + output_png_path = join(output_folder_path, f"stack_{output_type}.png") subprocess.check_call(["gdaldem", "color-relief", "-of", "PNG", input_tif_path, "-alpha", color_map_path, output_png_path, "-nearest_color_entry"]) + log.debug(f'Finished creating quicklook image for stack_{output_type}') -def _merge_stack(rows, cols, params): - """ - Merge stacking outputs +def assemble_tiles(s, dir, tiles, out_type, index=None): """ - # setup paths - xlks, _, crop = cf.transform_params(params) - base_unw_paths = [] + Function to reassemble tiles from numpy files in to a merged array - for p in Path(params[OUT_DIR]).rglob("*rlks_*cr.tif"): - if "dem" not in str(p): - base_unw_paths.append(str(p)) + :param tuple s: shape for merged array. + :param str dir: path to directory containing numpy tile files. + :param str out_type: product type string, used to construct numpy tile file name. + :param int index: array third dimension index to extract from 3D time series array tiles. - if not base_unw_paths: - raise ConfigException(f"No rlooked files available in {params[OBS_DIR]} or {params[OBS_DIR]}") - - if "tif" in base_unw_paths[0].split(".")[1]: - dest_tifs = base_unw_paths # cf.get_dest_paths(base_unw_paths, crop, params, xlks) - for i, dest_tif in enumerate(dest_tifs): - dest_tifs[i] = dest_tif.replace("_tif","") - else: - dest_tifs = base_unw_paths # cf.get_dest_paths(base_unw_paths, crop, params, xlks) + :return: merged_array: array assembled from all tiles. + :rtype: ndarray + """ + log.debug('Re-assembling tiles for {}'.format(out_type)) + # pre-allocate dest array + merged_array = np.empty(shape=s, dtype=np.float32) - # load previously saved prepread_ifgs dict - preread_ifgs_file = join(params[cf.TMPDIR], 'preread_ifgs.pk') - ifgs = cp.load(open(preread_ifgs_file, 'rb')) - tiles = shared.get_tiles(dest_tifs[0], rows, cols) + # loop over each tile, load and slot in to correct spot + for t in tiles: + tile_file = Path(join(dir, out_type + '_'+str(t.index)+'.npy')) + tile = np.load(file=tile_file) + if index is None: #2D array + merged_array[t.top_left_y:t.bottom_right_y, t.top_left_x:t.bottom_right_x] = tile + else: #3D array + merged_array[t.top_left_y:t.bottom_right_y, t.top_left_x:t.bottom_right_x] = tile[:, :, index] - # stacking aggregation - if mpiops.size >= 3: - [_save_stack(ifgs, params, tiles, out_type=t) - for i, t in enumerate(['stack_rate', 'stack_error', 'stack_samples']) - if i == mpiops.rank] - else: - if mpiops.rank == MASTER_PROCESS: - [_save_stack(ifgs, params, tiles, out_type=t) - for t in ['stack_rate', 'stack_error', 'stack_samples']] - mpiops.comm.barrier() + log.debug('Finished assembling tiles for {}'.format(out_type)) + return merged_array -def _save_stack(ifgs_dict, params, tiles, out_type): +def _save_merged_files(ifgs_dict, outdir, array, out_type, index=None, savenpy=None): """ - Save stacking outputs + Convenience function to save PyRate geotiff and numpy array files """ - log.info('Merging PyRate outputs {}'.format(out_type)) + log.debug('Saving PyRate outputs {}'.format(out_type)) gt, md, wkt = ifgs_dict['gt'], ifgs_dict['md'], ifgs_dict['wkt'] epochlist = ifgs_dict['epochlist'] - ifgs = [v for v in ifgs_dict.values() if isinstance(v, PrereadIfg)] - dest = os.path.join(params[cf.OUT_DIR], out_type + ".tif") - md[ifc.EPOCH_DATE] = epochlist.dates + + if out_type in ('tsincr', 'tscuml'): + epoch = epochlist.dates[index + 1] + dest = join(outdir, out_type + "_" + str(epoch) + ".tif") + npy_file = join(outdir, out_type + "_" + str(epoch) + ".npy") + # sequence position; first time slice is #0 + md['SEQUENCE_POSITION'] = index+1 + md[ifc.EPOCH_DATE] = epoch + else: + dest = join(outdir, out_type + ".tif") + npy_file = join(outdir, out_type + '.npy') + md[ifc.EPOCH_DATE] = epochlist.dates + if out_type == 'stack_rate': md[ifc.DATA_TYPE] = ifc.STACKRATE elif out_type == 'stack_error': md[ifc.DATA_TYPE] = ifc.STACKERROR - else: + elif out_type == 'stack_samples': md[ifc.DATA_TYPE] = ifc.STACKSAMP + elif out_type == 'tsincr': + md[ifc.DATA_TYPE] = ifc.INCR + else: #tscuml + md[ifc.DATA_TYPE] = ifc.CUML - rate = np.zeros(shape=ifgs[0].shape, dtype=np.float32) + shared.write_output_geotiff(md, gt, wkt, array, dest, np.nan) + if savenpy: + np.save(file=npy_file, arr=array) - for t in tiles: - rate_file = os.path.join(params[cf.TMPDIR], out_type + '_'+str(t.index)+'.npy') - rate_file = Path(rate_file) - rate_tile = np.load(file=rate_file) - rate[t.top_left_y:t.bottom_right_y, t.top_left_x:t.bottom_right_x] = rate_tile - shared.write_output_geotiff(md, gt, wkt, rate, dest, np.nan) - npy_rate_file = os.path.join(params[cf.OUT_DIR], out_type + '.npy') - np.save(file=npy_rate_file, arr=rate) + log.debug('Finished saving {}'.format(out_type)) - log.debug('Finished PyRate merging {}'.format(out_type)) - -def _merge_timeseries(rows, cols, params): +def _merge_setup(rows, cols, params): """ - Merge time series output + Convenience function for Merge set up steps """ + # setup paths xlks, _, crop = cf.transform_params(params) - base_unw_paths = [] - for p in Path(params[OUT_DIR]).rglob("*rlks_*cr.tif"): + for p in Path(params[cf.OUT_DIR]).rglob("*rlks_*cr.tif"): if "dem" not in str(p): base_unw_paths.append(str(p)) if "tif" in base_unw_paths[0].split(".")[1]: - dest_tifs = base_unw_paths # cf.get_dest_paths(base_unw_paths, crop, params, xlks) + dest_tifs = base_unw_paths # cf.get_dest_paths(base_unw_paths, crop, params, xlks) for i, dest_tif in enumerate(dest_tifs): dest_tifs[i] = dest_tif.replace("_tif", "") else: - dest_tifs = base_unw_paths # cf.get_dest_paths(base_unw_paths, crop, params, xlks) - - output_dir = params[cf.TMPDIR] - # load previously saved prepread_ifgs dict - preread_ifgs_file = join(output_dir, 'preread_ifgs.pk') - ifgs = cp.load(open(preread_ifgs_file, 'rb')) - - # metadata and projections - gt, md, wkt = ifgs['gt'], ifgs['md'], ifgs['wkt'] - epochlist = ifgs['epochlist'] - ifgs = [v for v in ifgs.values() if isinstance(v, PrereadIfg)] + dest_tifs = base_unw_paths # cf.get_dest_paths(base_unw_paths, crop, params, xlks) + # load previously saved preread_ifgs dict + preread_ifgs_file = join(params[cf.TMPDIR], 'preread_ifgs.pk') + ifgs_dict = pickle.load(open(preread_ifgs_file, 'rb')) + ifgs = [v for v in ifgs_dict.values() if isinstance(v, shared.PrereadIfg)] + shape = ifgs[0].shape tiles = shared.get_tiles(dest_tifs[0], rows, cols) - - # load the first tsincr file to determine the number of time series tifs - tsincr_file = os.path.join(output_dir, 'tsincr_0.npy') - - tsincr = np.load(file=tsincr_file) - - # pylint: disable=no-member - no_ts_tifs = tsincr.shape[2] - # we create 2 x no_ts_tifs as we are splitting tsincr and tscuml - # to all processes. - process_tifs = mpiops.array_split(range(2 * no_ts_tifs)) - - # depending on nvelpar, this will not fit in memory - # e.g. nvelpar=100, nrows=10000, ncols=10000, 32bit floats need 40GB memory - # 32 * 100 * 10000 * 10000 / 8 bytes = 4e10 bytes = 40 GB - # the double for loop helps us overcome the memory limit - log.info('Process {} writing {} timeseries tifs of ' - 'total {}'.format(mpiops.rank, len(process_tifs), no_ts_tifs * 2)) - for i in process_tifs: - tscum_g = np.empty(shape=ifgs[0].shape, dtype=np.float32) - if i < no_ts_tifs: - for n, t in enumerate(tiles): - _assemble_tiles(i, n, t, tscum_g, output_dir, 'tscuml') - md[ifc.EPOCH_DATE] = epochlist.dates[i + 1] - # sequence position; first time slice is #0 - md['SEQUENCE_POSITION'] = i+1 - dest = os.path.join(params[cf.OUT_DIR], - 'tscuml' + "_" + - str(epochlist.dates[i + 1]) + ".tif") - md[ifc.DATA_TYPE] = ifc.CUML - shared.write_output_geotiff(md, gt, wkt, tscum_g, dest, np.nan) - else: - tsincr_g = np.empty(shape=ifgs[0].shape, dtype=np.float32) - i %= no_ts_tifs - for n, t in enumerate(tiles): - _assemble_tiles(i, n, t, tsincr_g, output_dir, 'tsincr') - md[ifc.EPOCH_DATE] = epochlist.dates[i + 1] - # sequence position; first time slice is #0 - md['SEQUENCE_POSITION'] = i+1 - dest = os.path.join(params[cf.OUT_DIR], - 'tsincr' + "_" + str( - epochlist.dates[i + 1]) + ".tif") - md[ifc.DATA_TYPE] = ifc.INCR - shared.write_output_geotiff(md, gt, wkt, tsincr_g, dest, np.nan) - mpiops.comm.barrier() - log.debug('Process {} finished writing {} timeseries tifs of ' - 'total {}'.format(mpiops.rank, len(process_tifs), no_ts_tifs * 2)) - - -def _assemble_tiles(i, n, tile, tsincr_g, output_dir, outtype): - # pylint: disable=too-many-arguments - """ - A reusable time series tile assembly function - """ - tsincr_file = os.path.join(output_dir, '{}_{}.npy'.format(outtype, n)) - tsincr = np.load(file=tsincr_file) - tsincr_g[tile.top_left_y:tile.bottom_right_y, tile.top_left_x:tile.bottom_right_x] = tsincr[:, :, i] + return shape, tiles, ifgs_dict def _delete_tsincr_files(params): diff --git a/pyrate/prepifg.py b/pyrate/prepifg.py index 03fbdd084..a9803df24 100644 --- a/pyrate/prepifg.py +++ b/pyrate/prepifg.py @@ -20,7 +20,7 @@ # -*- coding: utf-8 -*- import os from subprocess import check_call -from typing import List +from typing import List, Tuple from pathlib import Path from joblib import Parallel, delayed import numpy as np @@ -29,6 +29,7 @@ from pyrate.core.prepifg_helper import PreprocessError from pyrate.core.logger import pyratelogger as log from pyrate.core.shared import output_tiff_filename +from pyrate.configuration import MultiplePaths GAMMA = 1 ROIPAC = 0 @@ -53,60 +54,58 @@ def main(params): if params[cf.DEM_FILE] is not None: # optional DEM conversion ifg_paths.append(params[cf.DEM_FILE_PATH]) + if params[cf.COH_MASK]: + ifg_paths.extend(params[cf.COHERENCE_FILE_PATHS]) + shared.mkdir_p(params[cf.OUT_DIR]) # create output dir - process_ifgs_paths = np.array_split(ifg_paths, mpiops.size)[mpiops.rank] + user_exts = (params[cf.IFG_XFIRST], params[cf.IFG_YFIRST], params[cf.IFG_XLAST], params[cf.IFG_YLAST]) + xlooks, ylooks, crop = cf.transform_params(params) + ifgs = [prepifg_helper.dem_or_ifg(p.converted_path) for p in ifg_paths] + exts = prepifg_helper.get_analysis_extent(crop, ifgs, xlooks, ylooks, user_exts=user_exts) - gtiff_paths = [p.converted_path for p in process_ifgs_paths] - do_prepifg(gtiff_paths, params) + process_ifgs_paths = np.array_split(ifg_paths, mpiops.size)[mpiops.rank] + do_prepifg(process_ifgs_paths, exts, params) mpiops.comm.barrier() log.info("Finished prepifg") -def do_prepifg(gtiff_paths: List[str], params: dict) -> None: +def do_prepifg(multi_paths: List[MultiplePaths], exts: Tuple[float, float, float, float], params: dict) -> None: """ Prepare interferograms by applying multilooking/cropping operations. - :param list gtiff_paths: List of full-res geotiffs - :param dict params: Parameters dictionary corresponding to config file """ # pylint: disable=expression-not-assigned parallel = params[cf.PARALLEL] if mpiops.size > 1: parallel = False - for f in gtiff_paths: - if not os.path.isfile(f): + for f in multi_paths: + if not os.path.isfile(f.converted_path): raise FileNotFoundError("Can not find geotiff: " + str(f) + ". Ensure you have converted your " "interferograms to geotiffs.") - ifgs = [prepifg_helper.dem_or_ifg(p) for p in gtiff_paths] - xlooks, ylooks, crop = cf.transform_params(params) - user_exts = (params[cf.IFG_XFIRST], params[cf.IFG_YFIRST], params[cf.IFG_XLAST], params[cf.IFG_YLAST]) - exts = prepifg_helper.get_analysis_extent(crop, ifgs, xlooks, ylooks, user_exts=user_exts) - thresh = params[cf.NO_DATA_AVERAGING_THRESHOLD] - if params[cf.LARGE_TIFS]: log.info("Using gdal system calls to process prepifg") - ifg = ifgs[0] + ifg = prepifg_helper.dem_or_ifg(multi_paths[0].converted_path) + ifg.open() + xlooks, ylooks = params[cf.IFG_LKSX], params[cf.IFG_LKSY] res_str = [xlooks * ifg.x_step, ylooks * ifg.y_step] res_str = ' '.join([str(e) for e in res_str]) if parallel: Parallel(n_jobs=params[cf.PROCESSES], verbose=50)( - delayed(__prepifg_system)( - crop, exts, gtiff_path, params, res_str, thresh, xlooks, ylooks) for gtiff_path in gtiff_paths - ) + delayed(__prepifg_system)(exts, gtiff_path, params, res_str) for gtiff_path in multi_paths) else: - for gtiff_path in gtiff_paths: - __prepifg_system(crop, exts, gtiff_path, params, res_str, thresh, xlooks, ylooks) + for m_path in multi_paths: + __prepifg_system(exts, m_path, params, res_str) else: if parallel: Parallel(n_jobs=params[cf.PROCESSES], verbose=50)( - delayed(_prepifg_multiprocessing)(p, xlooks, ylooks, exts, thresh, crop, params) for p in gtiff_paths + delayed(_prepifg_multiprocessing)(p, exts, params) for p in multi_paths ) else: - for gtiff_path in gtiff_paths: - _prepifg_multiprocessing(gtiff_path, xlooks, ylooks, exts, thresh, crop, params) + for m_path in multi_paths: + _prepifg_multiprocessing(m_path, exts, params) COMMON_OPTIONS = "-co BLOCKXSIZE=256 -co BLOCKYSIZE=256 -co TILED=YES --config GDAL_CACHEMAX=64 -q" @@ -114,8 +113,9 @@ def do_prepifg(gtiff_paths: List[str], params: dict) -> None: GDAL_CALC = 'gdal_calc_local.py' -def __prepifg_system(crop, exts, gtiff, params, res, thresh, xlooks, ylooks): - p, c, l = _prepifg_multiprocessing(gtiff, xlooks, ylooks, exts, thresh, crop, params) +def __prepifg_system(exts, gtiff, params, res): + thresh = params[cf.NO_DATA_AVERAGING_THRESHOLD] + p, c, l = _prepifg_multiprocessing(gtiff, exts, params) log.info("Multilooking {p} into {l}".format(p=p, l=l)) extents = ' '.join([str(e) for e in exts]) @@ -145,8 +145,8 @@ def __prepifg_system(crop, exts, gtiff, params, res, thresh, xlooks, ylooks): # coh masking check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} --overwrite -A {p_unset} -B {c}\t' - f'--calc=\"A*(B>={params[cf.COH_THRESH]})' - f'-99999*logical_or((B<{params[cf.COH_THRESH]}),isclose(A,0,atol=0.000001))\"\t' + f'--calc=\"A*(B>={params[cf.COH_THRESH]}) - ' + f'99999*logical_or((B<{params[cf.COH_THRESH]}), isclose(A,0,atol=0.000001))\"\t' f'--outfile={corrected_p}\t' f'--NoDataValue=nan', shell=True) else: @@ -156,7 +156,7 @@ def __prepifg_system(crop, exts, gtiff, params, res, thresh, xlooks, ylooks): f'--outfile={nan_frac}\t' f'--NoDataValue=nan', shell=True) check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} --overwrite -A {p_unset}\t' - f'--calc=\"A - 99999 *isclose(A, 0, atol=0.000001)\"\t' + f'--calc=\"A - 99999*isclose(A, 0, atol=0.000001)\"\t' f'--outfile={corrected_p}\t' f'--NoDataValue=nan', shell=True) @@ -194,6 +194,8 @@ def __update_meta_data(p_unset, c, l): else: if v == ifc.DEM: # it's a dem md_str = '-mo {k}={v}'.format(k=ifc.DATA_TYPE, v=ifc.MLOOKED_DEM) + elif v == ifc.COH: + md_str = '-mo {k}={v}'.format(k=ifc.DATA_TYPE, v=ifc.MULTILOOKED_COH) else: # it's an ifg md_str = '-mo {k}={v}'.format(k=ifc.DATA_TYPE, v=ifc.MULTILOOKED) for k, v in md.items(): @@ -202,31 +204,42 @@ def __update_meta_data(p_unset, c, l): ds = None -def _prepifg_multiprocessing(path, xlooks, ylooks, exts, thresh, crop, params): +def _prepifg_multiprocessing(m_path: MultiplePaths, exts: Tuple[float, float, float, float], params: dict): """ Multiprocessing wrapper for prepifg """ - processor = params[cf.PROCESSOR] # roipac, gamma or geotif - if (processor == GAMMA) or (processor == GEOTIF): - header = gamma.gamma_header(path, params) - elif processor == ROIPAC: - log.info("Warning: ROI_PAC support will be deprecated in a future PyRate release") - header = roipac.roipac_header(path, params) - else: - raise PreprocessError('Processor must be ROI_PAC (0) or GAMMA (1)') + xlooks, ylooks, crop = cf.transform_params(params) + thresh = params[cf.NO_DATA_AVERAGING_THRESHOLD] + header = find_header(m_path, params) + header[ifc.INPUT_TYPE] = m_path.input_type # If we're performing coherence masking, find the coherence file for this IFG. if params[cf.COH_MASK] and shared._is_interferogram(header): - coherence_path = cf.coherence_paths_for(path, params, tif=True) + coherence_path = cf.coherence_paths_for(m_path.converted_path, params, tif=True) coherence_thresh = params[cf.COH_THRESH] else: coherence_path = None coherence_thresh = None if params[cf.LARGE_TIFS]: - op = output_tiff_filename(path, params[cf.OUT_DIR]) + op = output_tiff_filename(m_path.converted_path, params[cf.OUT_DIR]) looks_path = cf.mlooked_path(op, ylooks, crop) - return path, coherence_path, looks_path + return m_path.converted_path, coherence_path, looks_path + else: + prepifg_helper.prepare_ifg(m_path.converted_path, xlooks, ylooks, exts, thresh, crop, + out_path=params[cf.OUT_DIR], header=header, coherence_path=coherence_path, + coherence_thresh=coherence_thresh) + + +def find_header(path: MultiplePaths, params: dict): + processor = params[cf.PROCESSOR] # roipac, gamma or geotif + tif_path = path.converted_path + if (processor == GAMMA) or (processor == GEOTIF): + header = gamma.gamma_header(tif_path, params) + elif processor == ROIPAC: + log.info("Warning: ROI_PAC support will be deprecated in a future PyRate release") + header = roipac.roipac_header(tif_path, params) else: - prepifg_helper.prepare_ifg(path, xlooks, ylooks, exts, thresh, crop, out_path=params[cf.OUT_DIR], - header=header, coherence_path=coherence_path, coherence_thresh=coherence_thresh) + raise PreprocessError('Processor must be ROI_PAC (0) or GAMMA (1)') + header[ifc.INPUT_TYPE] = path.input_type + return header diff --git a/pyrate/process.py b/pyrate/process.py index f9be74ee2..c9ff7700b 100644 --- a/pyrate/process.py +++ b/pyrate/process.py @@ -17,11 +17,11 @@ """ This Python module runs the main PyRate processing workflow """ -import logging import os from os.path import join import pickle as cp from collections import OrderedDict +from typing import List, Tuple import numpy as np from pyrate.core import (shared, algorithm, orbital, ref_phs_est as rpe, @@ -31,6 +31,8 @@ from pyrate.core.aps import wrap_spatio_temporal_filter from pyrate.core.shared import Ifg, PrereadIfg, get_tiles, mpi_vs_multiprocess_logging from pyrate.core.logger import pyratelogger as log +from pyrate.prepifg import find_header +from pyrate.configuration import MultiplePaths MASTER_PROCESS = 0 @@ -45,7 +47,7 @@ def _join_dicts(dicts): return assembled_dict -def _create_ifg_dict(dest_tifs, params, tiles): +def _create_ifg_dict(dest_tifs, params): """ 1. Convert ifg phase data into numpy binary files. 2. Save the preread_ifgs dict with information about the ifgs that are @@ -62,7 +64,6 @@ def _create_ifg_dict(dest_tifs, params, tiles): ifgs_dict = {} nifgs = len(dest_tifs) process_tifs = mpiops.array_split(dest_tifs) - shared.save_numpy_phase(dest_tifs, tiles, params) for d in process_tifs: ifg = shared._prep_ifg(d, params) ifgs_dict[d] = PrereadIfg(path=d, @@ -119,20 +120,22 @@ def _save_mst_tile(tile, i, preread_ifgs): mpiops.comm.barrier() -def _ref_pixel_calc(ifg_paths, params): +def _ref_pixel_calc(ifg_paths: List[str], params: dict) -> Tuple[int, int]: """ Wrapper for reference pixel calculation """ - refx = params[cf.REFX] - refy = params[cf.REFY] + + lon = params[cf.REFX] + lat = params[cf.REFY] ifg = Ifg(ifg_paths[0]) ifg.open(readonly=True) + # assume all interferograms have same projection and will share the same transform + transform = ifg.dataset.GetGeoTransform() - if refx == -1 or refy == -1: + if lon == -1 or lat == -1: log.info('Searching for best reference pixel location') - half_patch_size, thresh, grid = refpixel.ref_pixel_setup(ifg_paths, params) process_grid = mpiops.array_split(grid) refpixel.save_ref_pixel_blocks(process_grid, half_patch_size, ifg_paths, params) @@ -149,17 +152,25 @@ def _ref_pixel_calc(ifg_paths, params): "Reference pixel calculation returned an all nan slice!\n" "Cannot continue downstream computation. Please change reference pixel algorithm used before " "continuing.") + refy, refx = refpixel_returned # row first means first value is latitude + log.info('Selected reference pixel coordinate (x, y): ({}, {})'.format(refx, refy)) + lon, lat = refpixel.convert_pixel_value_to_geographic_coordinate(refx, refy, transform) + log.info('Selected reference pixel coordinate (lon, lat): ({}, {})'.format(lon, lat)) - refy, refx = refpixel_returned - - log.info('Selected reference pixel coordinate: ({}, {})'.format(refx, refy)) else: - log.info('Reusing reference pixel from config file: ({}, {})'.format(refx, refy)) + log.info('Using reference pixel from config file (lon, lat): ({}, {})'.format(lon, lat)) + log.warning("Ensure user supplied reference pixel values are in lon/lat") + refx, refy = refpixel.convert_geographic_coordinate_to_pixel_value(lon, lat, transform) + log.info('Converted reference pixel coordinate (x, y): ({}, {})'.format(refx, refy)) + + refpixel.update_refpix_metadata(ifg_paths, refx, refy, transform, params) + + log.debug("refpx, refpy: "+str(refx) + " " + str(refy)) ifg.close() - return refx, refy + return int(refx), int(refy) -def _orb_fit_calc(ifg_paths, params, preread_ifgs=None): +def _orb_fit_calc(multi_paths: List[MultiplePaths], params, preread_ifgs=None) -> None: """ MPI wrapper for orbital fit correction """ @@ -169,6 +180,7 @@ def _orb_fit_calc(ifg_paths, params, preread_ifgs=None): return log.info('Calculating orbital correction') + ifg_paths = [p.sampled_path for p in multi_paths] if preread_ifgs: # don't check except for mpi tests # perform some general error/sanity checks log.debug('Checking Orbital error correction status') @@ -186,7 +198,8 @@ def _orb_fit_calc(ifg_paths, params, preread_ifgs=None): # 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) + headers = [find_header(p, params) for p in multi_paths] + orbital.remove_orbital_error(ifg_paths, params, headers, preread_ifgs=preread_ifgs) mpiops.comm.barrier() log.debug('Finished Orbital error correction') @@ -195,7 +208,7 @@ def _ref_phase_estimation(ifg_paths, params, refpx, refpy): """ Wrapper for reference phase estimation. """ - log.info("Calculating reference phase") + log.info("Calculating reference phase and correcting each interferogram") if len(ifg_paths) < 2: raise rpe.ReferencePhaseError( "At least two interferograms required for reference phase correction ({len_ifg_paths} " @@ -203,7 +216,7 @@ def _ref_phase_estimation(ifg_paths, params, refpx, refpy): ) if mpiops.run_once(shared.check_correction_status, ifg_paths, ifc.PYRATE_REF_PHASE): - log.debug('Finished reference phase estimation') + log.debug('Finished reference phase correction') return if params[cf.REF_EST_METHOD] == 1: @@ -228,7 +241,7 @@ def _ref_phase_estimation(ifg_paths, params, refpx, refpy): np.save(file=ref_phs_file, arr=collected_ref_phs) else: mpiops.comm.Send(ref_phs, dest=MASTER_PROCESS, tag=mpiops.rank) - log.debug('Finished reference phase estimation') + log.debug('Finished reference phase correction') # Preserve old return value so tests don't break. if isinstance(ifg_paths[0], Ifg): @@ -280,23 +293,27 @@ def process_ifgs(ifg_paths, params, rows, cols): if mpiops.size > 1: # turn of multiprocessing during mpi jobs params[cf.PARALLEL] = False + outdir = params[cf.TMPDIR] + if not os.path.exists(outdir): + shared.mkdir_p(outdir) tiles = mpiops.run_once(get_tiles, ifg_paths[0], rows, cols) - preread_ifgs = _create_ifg_dict(ifg_paths, params=params, tiles=tiles) - # _mst_calc(ifg_paths, params, tiles, preread_ifgs) + preread_ifgs = _create_ifg_dict(ifg_paths, params=params) + # validate user supplied ref pixel + refpixel.validate_supplied_lat_lon(params) refpx, refpy = _ref_pixel_calc(ifg_paths, params) - log.debug("refpx, refpy: "+str(refpx) + " " + str(refpy)) - # remove non ifg keys _ = [preread_ifgs.pop(k) for k in ['gt', 'epochlist', 'md', 'wkt']] - _orb_fit_calc(ifg_paths, params, preread_ifgs) + multi_paths = params[cf.INTERFEROGRAM_FILES] + _orb_fit_calc(multi_paths, params, preread_ifgs) _ref_phase_estimation(ifg_paths, params, refpx, refpy) + shared.save_numpy_phase(ifg_paths, tiles, params) _mst_calc(ifg_paths, params, tiles, preread_ifgs) # spatio-temporal aps filter @@ -326,7 +343,7 @@ def _stack_calc(ifg_paths, params, vcmt, tiles, preread_ifgs): log.info('Stacking of tile {}'.format(t.index)) ifg_parts = [shared.IfgPart(p, t, preread_ifgs, params) for p in ifg_paths] mst_grid_n = np.load(os.path.join(output_dir, 'mst_mat_{}.npy'.format(t.index))) - rate, error, samples = stack.stack_rate(ifg_parts, params, vcmt, mst_grid_n) + rate, error, samples = stack.stack_rate_array(ifg_parts, params, vcmt, mst_grid_n) # declare file names np.save(file=os.path.join(output_dir, 'stack_rate_{}.npy'.format(t.index)), arr=rate) np.save(file=os.path.join(output_dir, 'stack_error_{}.npy'.format(t.index)), arr=error) diff --git a/requirements-dev.txt b/requirements-dev.txt index a3b342687..97cb4ee5e 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,10 +1,10 @@ ghp-import==0.5.5 -sphinx==2.1.2 +sphinx==3.0.3 sphinx_rtd_theme==0.4.3 -recommonmark==0.5.0 -matplotlib==3.1.1 -IPython==7.6.1 -sphinxcontrib-programoutput==0.14 +sphinxcontrib-programoutput==0.16 +recommonmark==0.6.0 +matplotlib==3.2.1 +ipython==7.14.0 setuptools==41.0.1 -twine==1.13.0 -wheel==0.33.4 +twine==3.1.1 +wheel==0.34.2 diff --git a/requirements-test.txt b/requirements-test.txt index 58f008b4e..7c1e08958 100644 --- a/requirements-test.txt +++ b/requirements-test.txt @@ -1,5 +1,5 @@ -pytest-cov==2.5.1 -coverage==4.5.3 -codecov==2.0.15 -tox==3.13.2 -pytest==5.0.1 +pytest-cov==2.8.1 +coverage==5.1 +codecov==2.1.0 +tox==3.15.0 +pytest==5.4.2 diff --git a/scripts/gdal_calc_local.py b/scripts/gdal_calc_local.py index a525ecbfd..b5093b1a0 100755 --- a/scripts/gdal_calc_local.py +++ b/scripts/gdal_calc_local.py @@ -286,7 +286,7 @@ def doit(opts, args): raise if myNDVs is not None: - myResult[numpy.isclose(myResult, -99999, atol=1e-6)] = myOutNDV + myResult[numpy.isclose(myResult, -99999, atol=1e-5)] = myOutNDV elif not isinstance(myResult, numpy.ndarray): myResult = numpy.ones( (nYValid,nXValid) ) * myResult diff --git a/setup.cfg b/setup.cfg index b93be300f..a1adf6209 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,10 @@ # Inside of setup.cfg [metadata] description-file = README.rst + +[aliases] +test = pytest + +[tool:pytest] +markers: + slow: Skipped unless --runslow passed diff --git a/tests/common.py b/tests/common.py index eabe0311a..f9d7a4c82 100644 --- a/tests/common.py +++ b/tests/common.py @@ -24,6 +24,7 @@ import stat import tempfile from os.path import join +from subprocess import check_output from pathlib import Path import numpy as np @@ -34,6 +35,14 @@ from pyrate.core.shared import (Ifg, nan_and_mm_convert, get_geotiff_header_info, write_output_geotiff, dem_or_ifg) from pyrate.constants import PYRATEPATH +from pyrate.configuration import Configuration + +TRAVIS = True if 'TRAVIS' in os.environ else False +PYTHON3P6 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.6') else False +PYTHON3P7 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.7') else False +PYTHON3P8 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.8') else False +GDAL_VERSION = check_output(["gdal-config", "--version"]).decode(encoding="utf-8").split('\n')[0] + TEMPDIR = tempfile.gettempdir() TESTDIR = join(PYRATEPATH, 'tests') @@ -159,6 +168,26 @@ def assert_tifs_equal(tif1, tif2): sds = None +def copy_small_ifg_file_list(): + temp_dir = tempfile.mkdtemp() + move_files(SML_TEST_TIF, temp_dir, file_type='*.tif', copy=True) + datafiles = glob.glob(join(temp_dir, "*.tif")) + for d in datafiles: + Path(d).chmod(0o664) # assign write permission as conv2tif output is readonly + return temp_dir, datafiles + + +def copy_and_setup_small_data(): + temp_dir, datafiles = copy_small_ifg_file_list() + datafiles.sort() + ifgs = [dem_or_ifg(i) for i in datafiles] + + for i in ifgs: + i.open() + i.nodata_value = 0 + return temp_dir, ifgs + + def small_ifg_file_list(datafiles=None): """Returns the file list of all the .tif files after prepifg conversion input phase data is in radians; these ifgs are in radians - not converted to mm""" @@ -274,9 +303,12 @@ def reconstruct_mst(shape, tiles, output_dir): return mst -def move_files(source_dir, dest_dir, file_type='*.tif'): +def move_files(source_dir, dest_dir, file_type='*.tif', copy=False): for filename in glob.glob(os.path.join(source_dir, file_type)): - shutil.move(filename, dest_dir) + if copy: + shutil.copy(filename, dest_dir) + else: + shutil.move(filename, dest_dir) def assert_ifg_phase_equal(ifg_path1, ifg_path2): @@ -382,12 +414,13 @@ def write_timeseries_geotiff(ifgs, params, tsincr, pr_type): def calculate_stack_rate(ifgs, params, vcmt, mst_mat=None): # log.info('Calculating stacked rate') - res = stack.stack_rate(ifgs, params, vcmt, mst_mat) + res = stack.stack_rate_array(ifgs, params, vcmt, mst_mat) for r in res: if r is None: raise ValueError('TODO: bad value') - rate, error, samples = res + r, e, samples = res + rate, error = stack.mask_rate(r, e, params['maxsig']) write_stackrate_tifs(ifgs, params, res) # log.info('Stacked rate calculated') return rate, error, samples @@ -459,7 +492,7 @@ def copytree(src, dst, symlinks=False, ignore=None): def pre_prepare_ifgs(ifg_paths, params): """ - Open ifg for reading + nan and mm convert ifgs """ ifgs = [Ifg(p) for p in ifg_paths] for i in ifgs: @@ -480,6 +513,8 @@ def assert_two_dirs_equal(dir1, dir2, ext, num_files=None): if num_files is not None: assert len(dir1_files) == num_files assert len(dir2_files) == num_files + else: + assert len(dir1_files) == len(dir2_files) if dir1_files[0].suffix == '.tif': for m_f, s_f in zip(dir1_files, dir2_files): assert m_f.name == s_f.name @@ -489,6 +524,8 @@ def assert_two_dirs_equal(dir1, dir2, ext, num_files=None): for m_f, s_f in zip(dir1_files, dir2_files): assert m_f.name == s_f.name np.testing.assert_array_almost_equal(np.load(m_f), np.load(s_f)) + elif dir1_files[0].suffix in {'.kml', '.png'}: + return else: raise diff --git a/tests/conftest.py b/tests/conftest.py index 3a2d39d58..19240a3ba 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -3,7 +3,7 @@ import string import tempfile import pytest -from pyrate.core import mpiops, config as cf, prepifg_helper +from pyrate.core import mpiops, config as cf, prepifg_helper, shared from pyrate.configuration import Configuration from tests.common import TEST_CONF_ROIPAC, TEST_CONF_GAMMA from tests.common import ROIPAC_SYSTEM_CONF, GAMMA_SYSTEM_CONF, GEOTIF_SYSTEM_CONF, SML_TEST_COH_LIST @@ -46,11 +46,21 @@ def fin(): return mpiops.comm +@pytest.fixture(params=[0, 1]) +def coh_mask(request): + return request.param + + @pytest.fixture(params=[1, 2]) def ref_est_method(request): return request.param +@pytest.fixture(params=[(-1, -1), (150.941666654, -34.218333314)]) +def ref_pixel(request): + return request.param + + @pytest.fixture(params=[1, 3]) def orbfit_lks(request): return request.param @@ -87,6 +97,8 @@ def params(conf_file): @pytest.fixture def gamma_params(): params = Configuration(TEST_CONF_GAMMA).__dict__ + shutil.rmtree(params[cf.OUT_DIR]) + shared.mkdir_p(params[cf.OUT_DIR]) yield params shutil.rmtree(params[cf.OUT_DIR]) @@ -107,7 +119,7 @@ def roipac_or_gamma_conf(request): def gamma_conf(request): params = Configuration(TEST_CONF_GAMMA).__dict__ yield request.param - shutil.rmtree(params[cf.OUT_DIR]) + shutil.rmtree(params[cf.OUT_DIR], ignore_errors=True) @pytest.fixture diff --git a/tests/test_coherence.py b/tests/test_coherence.py index 7e1f92459..b93ed4f08 100644 --- a/tests/test_coherence.py +++ b/tests/test_coherence.py @@ -1,4 +1,5 @@ import os +import stat import tempfile import numpy as np from osgeo import osr @@ -24,14 +25,14 @@ def test_small_data_coherence(gamma_params): conv2tif.main(gamma_params) for i in ifg_multilist: - p = i.converted_path - ifg = pyrate.core.shared.dem_or_ifg(data_path=p) + p = Path(i.converted_path) + p.chmod(0o664) # assign write permission as conv2tif output is readonly + ifg = pyrate.core.shared.dem_or_ifg(data_path=p.as_posix()) if not isinstance(ifg, Ifg): continue ifg.open() - # now do coherence masking and compare - ifg = pyrate.core.shared.dem_or_ifg(data_path=p) + ifg = pyrate.core.shared.dem_or_ifg(data_path=p.as_posix()) ifg.open() converted_coh_file_path = cf.coherence_paths_for(p, gamma_params, tif=True) gdal_python.coherence_masking(ifg.dataset, diff --git a/tests/test_config.py b/tests/test_config.py index 3669d8fc3..51468e881 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -29,19 +29,10 @@ from tests.common import SML_TEST_CONF, SML_TEST_TIF from tests.common import TEST_CONF_ROIPAC, TEST_CONF_GAMMA from pyrate.core import config -from pyrate.core.config import ( - validate_parameters, validate_optional_parameters, validate_epochs, - validate_obs_thresholds, validate_ifgs, validate_gamma_headers, - validate_coherence_files, validate_tifs_exist, validate_minimum_epochs, - validate_epoch_thresholds, validate_epoch_cutoff, - validate_crop_parameters, validate_slpf_cutoff, - validate_reference_pixel_search_windows, - validate_multilook_parameters, - validate_prepifg_tifs_exist, - _get_temporal_info, _get_prepifg_info, write_config_file) - -from pyrate.constants import SIXTEEN_DIGIT_EPOCH_PAIR, TWELVE_DIGIT_EPOCH_PAIR, \ - EIGHT_DIGIT_EPOCH +from pyrate.core.config import write_config_file + +from pyrate.constants import SIXTEEN_DIGIT_EPOCH_PAIR, TWELVE_DIGIT_EPOCH_PAIR, EIGHT_DIGIT_EPOCH +from pyrate.core import config as cf from pyrate.core.config import ( _COHERENCE_VALIDATION, @@ -321,238 +312,6 @@ def validate(key, value): self.assertFalse(validate(TIME_SERIES_METHOD, 0)) self.assertFalse(validate(TIME_SERIES_METHOD, 3)) - def test_gamma_validation_only_performed_when_on(self): - self.params[PROCESSOR] = 1 - self.params[HDR_FILE_LIST] = None - with pytest.raises(ConfigException): - validate_optional_parameters(self.params) - - self.params[PROCESSOR] = 0 - validate_parameters(self.params) - - def test_coherence_validation_only_performed_when_on(self): - self.params[COH_MASK] = 1 - self.params[COH_FILE_LIST] = None - with pytest.raises(ConfigException): - validate_optional_parameters(self.params) - - self.params[COH_MASK] = 0 - validate_optional_parameters(self.params) - - def test_orbital_validation_only_performed_when_on(self): - self.params[ORBITAL_FIT] = 1 - self.params[ORBITAL_FIT_METHOD] = 0 - with pytest.raises(ConfigException): - validate_optional_parameters(self.params) - - self.params[ORBITAL_FIT] = 0 - validate_optional_parameters(self.params) - - def test_apsest_validation_only_performed_when_on(self): - self.params[APSEST] = 1 - self.params[TLPF_METHOD] = 0 - with pytest.raises(ConfigException): - validate_optional_parameters(self.params) - - self.params[APSEST] = 0 - validate_optional_parameters(self.params) - - def test_time_series_validation_only_performed_when_on(self): - self.params[TIME_SERIES_CAL] = 1 - self.params[TIME_SERIES_PTHRESH] = 0 - with pytest.raises(ConfigException): - validate_optional_parameters(self.params) - - self.params[TIME_SERIES_CAL] = 0 - validate_optional_parameters(self.params) - - def test_epochs_in_gamma_obs(self): - validate_epochs(self.params[IFG_FILE_LIST], SIXTEEN_DIGIT_EPOCH_PAIR) - self.params[IFG_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'bad_epochs_ifms_17') - with pytest.raises(ConfigException): - validate_epochs(self.params[IFG_FILE_LIST], SIXTEEN_DIGIT_EPOCH_PAIR) - - def test_epochs_in_roipac_obs(self): - validate_epochs(self.roipac_params[IFG_FILE_LIST], TWELVE_DIGIT_EPOCH_PAIR) - self.roipac_params[IFG_FILE_LIST] = os.path.join(common.SML_TEST_OBS, 'bad_epochs_ifms_17') - with pytest.raises(ConfigException): - validate_epochs(self.roipac_params[IFG_FILE_LIST], TWELVE_DIGIT_EPOCH_PAIR) - - def test_epochs_in_gamma_headers(self): - validate_epochs(self.params[HDR_FILE_LIST], EIGHT_DIGIT_EPOCH) - self.params[HDR_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'bad_epochs_headers') - with pytest.raises(ConfigException): - validate_epochs(self.params[HDR_FILE_LIST], EIGHT_DIGIT_EPOCH) - - def test_epochs_in_coherence_files(self): - self.params[COH_MASK] = 1 - self.params[COH_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'ifms_17') - validate_epochs(self.params[COH_FILE_LIST], SIXTEEN_DIGIT_EPOCH_PAIR) - self.params[COH_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'bad_epochs_ifms_17') - with pytest.raises(ConfigException): - validate_epochs(self.params[COH_FILE_LIST], SIXTEEN_DIGIT_EPOCH_PAIR) - - def test_ifgs_exist(self): - validate_ifgs(self.params[IFG_FILE_LIST]) - self.params[IFG_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'bad_epochs_ifms_17') - with pytest.raises(ConfigException): - validate_ifgs(self.params[IFG_FILE_LIST]) - - def test_gamma_headers_exist(self): - validate_gamma_headers(self.params[IFG_FILE_LIST], self.params[HDR_FILE_LIST]) - self.params[HDR_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'bad_epochs_headers') - with pytest.raises(ConfigException): - validate_gamma_headers(self.params[IFG_FILE_LIST], self.params[HDR_FILE_LIST]) - - def test_coherence_files_exist(self): - self.params[COH_MASK] = 1 - self.params[COH_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'ifms_17') - validate_coherence_files(self.params[IFG_FILE_LIST], self.params) - self.params[COH_FILE_LIST] = os.path.join(common.SML_TEST_GAMMA, 'bad_epochs_ifms_17') - tmp_config = tempfile.mktemp('.conf') - write_config_file(self.params, output_conf_file=tmp_config) - with pytest.raises(ConfigException): - Configuration(tmp_config).__dict__ - - def test_obs_thresholds(self): - self.params[TIME_SERIES_PTHRESH] = 16 - self.params[TLPF_PTHR] = 16 - validate_obs_thresholds(self.params[IFG_FILE_LIST], self.params) - - self.params[TIME_SERIES_PTHRESH] = 18 - self.params[TLPF_PTHR] = 18 - with pytest.raises(ConfigException): - validate_obs_thresholds(self.params[IFG_FILE_LIST], self.params) - - -class TestConfigValidationWithFullResGeotiffs(unittest.TestCase): - @classmethod - def setUpClass(cls): - from pyrate import conv2tif - cls.params = Configuration(TEST_CONF_GAMMA).__dict__ - conv2tif.main(cls.params) - - @classmethod - def tearDownClass(cls): - common.remove_tifs(cls.params[OBS_DIR]) - - def setUp(self): - self.params = config.get_config_params(TEST_CONF_GAMMA, step='prepifg') - self.dummy_dir = '/i/should/not/exist' - crop_opts = config._crop_opts(self.params) - self.min_extents, self.n_cols, self.n_rows = \ - config._get_fullres_info(self.params[IFG_FILE_LIST], self.params[OUT_DIR], crop_opts) - if os.path.exists(self.dummy_dir): - raise IOError("{dummy_dir} needs to not exist for test purposes.") - - def test_validate_tifs_exist(self): - validate_tifs_exist(self.params[IFG_FILE_LIST], self.params[OUT_DIR]) - with pytest.raises(ConfigException): - validate_tifs_exist(self.params[IFG_FILE_LIST], self.dummy_dir) - - def test_validate_ifg_crop_coordinates(self): - self.params[IFG_CROP_OPT] = 3 - self.params[IFG_XFIRST] = 150.0 - with pytest.raises(ConfigException): - validate_crop_parameters(self.min_extents, self.params) - self.params[IFG_XFIRST] = 150.92 - - self.params[IFG_XLAST] = 150.95 - with pytest.raises(ConfigException): - validate_crop_parameters(self.min_extents, self.params) - self.params[IFG_XLAST] = 150.94 - - self.params[IFG_YFIRST] = -34.16 - with pytest.raises(ConfigException): - validate_crop_parameters(self.min_extents, self.params) - self.params[IFG_YFIRST] = -34.18 - - self.params[IFG_YLAST] = -34.24 - with pytest.raises(ConfigException): - validate_crop_parameters(self.min_extents, self.params) - - def test_validate_ifg_multilook(self): - self.params[IFG_LKSX] = 48 - self.params[IFG_LKSY] = 73 - with pytest.raises(ConfigException): - validate_multilook_parameters(self.n_cols, self.n_rows, - IFG_LKSX, IFG_LKSY, self.params) - - -class TestConfigValidationWithPrepifgGeotiffs(unittest.TestCase): - @classmethod - def setUpClass(cls): - from pyrate import conv2tif, prepifg - cls.params = Configuration(TEST_CONF_GAMMA).__dict__ - conv2tif.main(cls.params) - prepifg.main(cls.params) - - @classmethod - def tearDownClass(cls): - common.remove_tifs(cls.params[OBS_DIR]) - common.remove_tifs(cls.params[OUT_DIR]) - - def setUp(self): - self.params = config.get_config_params(TEST_CONF_GAMMA, - 'process') - self.dummy_dir = '/i/should/not/exist' - crop_opts = config._crop_opts(self.params) - self.extents, self.n_cols, self.n_rows, _ = \ - _get_prepifg_info(self.params[IFG_FILE_LIST], self.params[OUT_DIR], self.params) - self.n_epochs, self.max_span = \ - _get_temporal_info(self.params[IFG_FILE_LIST], self.params[OUT_DIR]) - if os.path.exists(self.dummy_dir): - raise IOError("{dummy_dir} needs to not exist for test purposes.") - - def test_validate_prepifg_tifs_exist(self): - validate_prepifg_tifs_exist(self.params[IFG_FILE_LIST], self.params[OUT_DIR], self.params) - with pytest.raises(ConfigException): - self.params[IFG_LKSX] = 100 - validate_prepifg_tifs_exist(self.params[IFG_FILE_LIST], self.params[OUT_DIR], self.params) - - def test_validate_epoch_thresholds(self): - with pytest.raises(ConfigException): - validate_minimum_epochs(self.n_epochs, 50) - - validate_epoch_thresholds(self.n_epochs, self.params) - self.params[LR_PTHRESH] = 20 - with pytest.raises(ConfigException): - validate_epoch_thresholds(self.n_epochs, self.params) - self.params[LR_PTHRESH] = 5 - - self.params[SLPF_CUTOFF] = 1000 - with pytest.raises(ConfigException): - validate_epoch_cutoff(self.max_span, SLPF_CUTOFF, self.params) - - - - def test_validate_search_windows(self): - self.params[REF_CHIP_SIZE] = 21 - self.params[REFNX] = 2 - self.params[REFNY] = 3 - validate_reference_pixel_search_windows(self.n_cols, self.n_rows, self.params) - - self.params[REFNX] = 3 - with pytest.raises(ConfigException): - validate_reference_pixel_search_windows(self.n_cols, self.n_rows, self.params) - - self.params[REFNX] = 2 - self.params[REFNY] = 4 - with pytest.raises(ConfigException): - validate_reference_pixel_search_windows(self.n_cols, self.n_rows, self.params) - - def test_validate_multilook_parameters(self): - self.params[ORBITAL_FIT_LOOKS_X] = 48 - self.params[ORBITAL_FIT_LOOKS_Y] = 73 - with pytest.raises(ConfigException): - validate_multilook_parameters(self.n_cols, self.n_rows, - ORBITAL_FIT_LOOKS_X, ORBITAL_FIT_LOOKS_Y, self.params) - - def test_validate_slpf_cutoff(self): - self.params[SLPF_CUTOFF] = 9999 - with pytest.raises(ConfigException): - validate_slpf_cutoff(self.extents, self.params) - class ConfigTest(unittest.TestCase): diff --git a/tests/test_conv2tif.py b/tests/test_conv2tif.py index 5c67f0978..6954180a8 100644 --- a/tests/test_conv2tif.py +++ b/tests/test_conv2tif.py @@ -1,10 +1,16 @@ import os +import shutil import pytest import glob import copy +import itertools +from pathlib import Path import pyrate.core.config as cf -from pyrate import conv2tif, prepifg +from pyrate.core.shared import Ifg, DEM +from pyrate.core import ifgconstants as ifc +from pyrate import conv2tif, prepifg, configuration +from tests.common import manipulate_test_conf def test_dem_and_incidence_not_converted(gamma_params): @@ -18,6 +24,38 @@ def test_dem_and_incidence_not_converted(gamma_params): assert len(dem_tif) == 0 +def test_conv2tif_file_types(tempdir, gamma_conf): + tdir = Path(tempdir()) + params = manipulate_test_conf(gamma_conf, tdir) + params[cf.COH_MASK] = 1 + output_conf_file = 'conf.conf' + output_conf = tdir.joinpath(output_conf_file) + cf.write_config_file(params=params, output_conf_file=output_conf) + params_s = configuration.Configuration(output_conf).__dict__ + conv2tif.main(params_s) + ifg_files = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_ifg.tif')) + coh_files = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_coh.tif')) + dem_file = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_dem.tif'))[0] + # assert coherence and ifgs have correct metadata + for i in itertools.chain(*[ifg_files, coh_files]): + ifg = Ifg(i) + ifg.open() + md = ifg.meta_data + if i.name.endswith('_ifg.tif'): + assert md[ifc.DATA_TYPE] == ifc.ORIG + continue + if i.name.endswith('_coh.tif'): + assert md[ifc.DATA_TYPE] == ifc.COH + continue + + # assert dem has correct metadata + dem = DEM(dem_file.as_posix()) + dem.open() + md = dem.dataset.GetMetadata() + assert md[ifc.DATA_TYPE] == ifc.DEM + shutil.rmtree(tdir) + + def test_tifs_placed_in_out_dir(gamma_params): # Test no tifs in obs dir tifs = glob.glob(os.path.join(gamma_params[cf.OUT_DIR], '*.tif')) @@ -33,6 +71,8 @@ def test_num_gamma_tifs_equals_num_unws(gamma_params): # 17 unws + dem assert len(gtifs) == 18 + for g, _ in gtifs: # assert all output from conv2tfi are readonly + assert Path(g).stat().st_mode == 33060 def test_num_roipac_tifs_equals_num_unws(roipac_params): gtifs = conv2tif.main(roipac_params) diff --git a/tests/test_covariance.py b/tests/test_covariance.py index 63a8486da..76d51b754 100644 --- a/tests/test_covariance.py +++ b/tests/test_covariance.py @@ -17,6 +17,7 @@ This Python module contains tests for the covariance.py PyRate module. """ import os +from pathlib import Path import shutil import sys import tempfile @@ -30,6 +31,8 @@ from pyrate.core.covariance import cvd, get_vcmt, RDist from pyrate.configuration import Configuration import pyrate.core.orbital +from pyrate.core import roipac +from pyrate.core.config import parse_namelist from tests import common from tests.common import (small5_mock_ifgs, small5_ifgs, TEST_CONF_ROIPAC, small_data_setup, prepare_ifgs_without_phase) @@ -179,7 +182,6 @@ class LegacyEqualityTest(unittest.TestCase): @classmethod def setUpClass(cls): - params = Configuration(TEST_CONF_ROIPAC).__dict__ cls.temp_out_dir = tempfile.mkdtemp() sys.argv = ['prepifg.py', TEST_CONF_ROIPAC] @@ -190,13 +192,15 @@ def setUpClass(cls): conv2tif.main(params) prepifg.main(params) cls.params = params - xlks, ylks, crop = cf.transform_params(params) - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], - params[cf.OBS_DIR]) - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, params, xlks) + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + dest_paths = [c.converted_path for c in params[cf.INTERFEROGRAM_FILES]] + dest_paths = dest_paths[:-2] + for i in dest_paths: + Path(i).chmod(0o664) # assign write permission as conv2tif output is readonly ifgs = common.pre_prepare_ifgs(dest_paths, params) refx, refy = process._ref_pixel_calc(dest_paths, params) - pyrate.core.orbital.remove_orbital_error(ifgs, params) + headers = [roipac.roipac_header(i, cls.params) for i in base_ifg_paths] + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) ifgs = prepare_ifgs_without_phase(dest_paths, params) for ifg in ifgs: ifg.close() @@ -205,8 +209,7 @@ def setUpClass(cls): r_dist = RDist(ifgs[0])() ifgs[0].close() # Calculate interferogram noise - cls.maxvar = [cvd(i, params, r_dist, calc_alpha=True, - save_acg=True, write_vals=True)[0] for i in dest_paths] + cls.maxvar = [cvd(i, params, r_dist, calc_alpha=True, save_acg=True, write_vals=True)[0] for i in dest_paths] cls.vcmt = get_vcmt(ifgs, cls.maxvar) for ifg in ifgs: ifg.close() @@ -226,8 +229,7 @@ def test_legacy_maxvar_equality_small_test_files(self): def test_legacy_vcmt_equality_small_test_files(self): from tests.common import SML_TEST_DIR LEGACY_VCM_DIR = os.path.join(SML_TEST_DIR, 'vcm') - legacy_vcm = np.genfromtxt(os.path.join(LEGACY_VCM_DIR, - 'vcmt.csv'), delimiter=',') + legacy_vcm = np.genfromtxt(os.path.join(LEGACY_VCM_DIR, 'vcmt.csv'), delimiter=',') np.testing.assert_array_almost_equal(legacy_vcm, self.vcmt, decimal=3) def test_metadata(self): @@ -243,8 +245,7 @@ def test_save_cvd_data(self): if not ifg.is_open: ifg.open() data_file = join(self.params[cf.TMPDIR], - 'cvd_data_{b}.npy'.format( - b=basename(ifg.data_path).split('.')[0])) + 'cvd_data_{b}.npy'.format(b=basename(ifg.data_path).split('.')[0])) assert isfile(data_file) if __name__ == "__main__": diff --git a/tests/test_data/merge/colourmap.txt b/tests/test_data/merge/colourmap.txt deleted file mode 100644 index eb826c22e..000000000 --- a/tests/test_data/merge/colourmap.txt +++ /dev/null @@ -1,256 +0,0 @@ -nan 0 0 0 0 --16.042793273926 103 0 31 255 --15.917458951473455 162 25 43 255 --15.792124629020908 201 72 67 255 --15.66679030656836 228 130 104 255 --15.541455984115814 246 182 154 255 --15.416121661663267 251 221 205 255 --15.29078733921072 242 239 238 255 --15.165453016758173 213 231 240 255 --15.040118694305626 166 207 227 255 --14.914784371853079 107 172 208 255 --14.789450049400532 58 132 187 255 --14.664115726947985 29 92 154 255 --14.538781404495438 103 0 31 255 --14.413447082042891 106 1 31 255 --14.288112759590344 109 2 32 255 --14.162778437137797 112 3 32 255 --14.03744411468525 115 4 33 255 --13.912109792232703 118 5 33 255 --13.786775469780157 120 6 34 255 --13.66144114732761 123 7 34 255 --13.536106824875063 126 8 35 255 --13.410772502422516 129 9 35 255 --13.285438179969969 132 10 36 255 --13.160103857517422 135 11 36 255 --13.034769535064875 137 12 37 255 --12.909435212612328 140 13 38 255 --12.784100890159781 143 14 38 255 --12.658766567707234 146 16 39 255 --12.533432245254687 148 17 39 255 --12.40809792280214 151 18 40 255 --12.282763600349593 153 20 41 255 --12.157429277897046 156 21 41 255 --12.0320949554445 158 23 42 255 --11.906760632991952 161 24 43 255 --11.781426310539405 163 26 44 255 --11.656091988086859 166 27 44 255 --11.530757665634312 168 29 45 255 --11.405423343181765 170 31 46 255 --11.280089020729218 172 33 47 255 --11.15475469827667 174 35 48 255 --11.029420375824124 177 37 49 255 --10.904086053371577 179 39 50 255 --10.77875173091903 180 41 51 255 --10.653417408466483 182 43 52 255 --10.528083086013936 184 45 53 255 --10.402748763561389 186 48 54 255 --10.277414441108842 188 50 55 255 --10.152080118656295 190 52 57 255 --10.026745796203748 191 55 58 255 --9.901411473751201 193 57 59 255 --9.776077151298654 194 60 60 255 --9.650742828846107 196 62 62 255 --9.52540850639356 198 65 63 255 --9.400074183941014 199 68 64 255 --9.274739861488467 200 70 66 255 --9.14940553903592 202 73 67 255 --9.024071216583373 203 76 69 255 --8.898736894130826 205 78 70 255 --8.773402571678279 206 81 72 255 --8.648068249225732 207 84 73 255 --8.522733926773185 209 87 75 255 --8.397399604320638 210 89 76 255 --8.272065281868091 211 92 78 255 --8.146730959415544 213 95 80 255 --8.021396636962997 214 98 82 255 --7.89606231451045 215 100 83 255 --7.770727992057903 217 103 85 255 --7.645393669605356 218 106 87 255 --7.520059347152809 219 109 89 255 --7.3947250247002625 220 111 91 255 --7.2693907022477156 222 114 92 255 --7.144056379795169 223 117 94 255 --7.018722057342622 224 119 96 255 --6.893387734890075 225 122 98 255 --6.768053412437528 226 125 100 255 --6.642719089984981 227 127 102 255 --6.517384767532434 228 130 104 255 --6.392050445079887 230 133 106 255 --6.26671612262734 231 135 108 255 --6.141381800174793 232 138 111 255 --6.016047477722246 233 141 113 255 --5.890713155269699 234 143 115 255 --5.765378832817152 235 146 117 255 --5.640044510364605 236 148 119 255 --5.514710187912058 236 151 122 255 --5.389375865459511 237 153 124 255 --5.2640415430069645 238 156 126 255 --5.1387072205544175 239 158 128 255 --5.013372898101871 240 161 131 255 --4.888038575649324 241 163 133 255 --4.762704253196777 241 165 136 255 --4.63736993074423 242 168 138 255 --4.512035608291683 243 170 140 255 --4.386701285839136 243 172 143 255 --4.261366963386589 244 175 145 255 --4.136032640934042 245 177 148 255 --4.010698318481495 245 179 150 255 --3.885363996028948 246 181 153 255 --3.760029673576401 246 184 155 255 --3.6346953511238542 247 186 158 255 --3.5093610286713073 247 188 160 255 --3.3840267062187603 248 190 163 255 --3.2586923837662134 248 192 165 255 --3.1333580613136665 248 194 168 255 --3.0080237388611195 249 196 170 255 --2.8826894164085726 249 198 173 255 --2.7573550939560256 249 200 175 255 --2.6320207715034787 249 202 178 255 --2.5066864490509317 250 204 180 255 --2.381352126598385 250 205 183 255 --2.256017804145838 250 207 185 255 --2.130683481693291 250 209 188 255 --2.005349159240744 250 211 190 255 --1.880014836788197 250 212 192 255 --1.75468051433565 250 214 195 255 --1.6293461918831031 251 216 197 255 --1.5040118694305562 251 217 199 255 --1.3786775469780093 251 219 201 255 --1.2533432245254623 251 220 204 255 --1.1280089020729154 251 222 206 255 --1.0026745796203684 250 223 208 255 --0.8773402571678215 250 224 210 255 --0.7520059347152745 250 226 212 255 --0.6266716122627276 250 227 214 255 --0.5013372898101807 250 228 216 255 --0.3760029673576337 250 229 218 255 --0.2506686449050868 249 230 219 255 --0.12533432245253984 249 231 221 255 -7.105427357601002e-15 249 232 223 255 -0.12533432245255227 248 233 224 255 -0.250668644905101 248 234 226 255 -0.3760029673576497 248 235 227 255 -0.5013372898101949 247 236 229 255 -0.62667161226274 247 236 230 255 -0.7520059347152888 246 237 232 255 -0.8773402571678375 245 238 233 255 -1.0026745796203826 245 238 234 255 -1.1280089020729278 244 238 235 255 -1.2533432245254765 243 239 236 255 -1.3786775469780252 243 239 237 255 -1.5040118694305704 242 239 238 255 -1.6293461918831156 241 239 239 255 -1.7546805143356643 240 240 239 255 -1.880014836788213 239 240 240 255 -2.005349159240758 238 240 240 255 -2.1306834816933033 237 239 241 255 -2.256017804145852 236 239 241 255 -2.381352126598401 234 239 242 255 -2.506686449050946 233 239 242 255 -2.632020771503491 232 238 242 255 -2.75735509395604 231 238 242 255 -2.8826894164085886 229 238 242 255 -3.0080237388611337 228 237 242 255 -3.133358061313679 227 237 242 255 -3.2586923837662276 225 236 242 255 -3.3840267062187763 224 235 242 255 -3.5093610286713215 222 235 242 255 -3.6346953511238667 220 234 241 255 -3.7600296735764154 219 233 241 255 -3.885363996028964 217 233 241 255 -4.010698318481509 215 232 241 255 -4.136032640934054 213 231 240 255 -4.261366963386603 212 230 240 255 -4.386701285839152 210 229 239 255 -4.512035608291697 208 228 239 255 -4.637369930744242 206 227 238 255 -4.762704253196791 204 226 238 255 -4.88803857564934 202 225 237 255 -5.013372898101885 200 224 237 255 -5.13870722055443 198 223 236 255 -5.264041543006979 196 222 236 255 -5.389375865459527 194 221 235 255 -5.514710187912073 191 220 235 255 -5.640044510364618 189 219 234 255 -5.7653788328171665 187 218 233 255 -5.890713155269715 185 217 233 255 -6.01604747772226 182 215 232 255 -6.1413818001748055 180 214 232 255 -6.266716122627354 178 213 231 255 -6.392050445079903 175 212 230 255 -6.517384767532448 173 210 229 255 -6.642719089984993 170 209 229 255 -6.768053412437542 168 208 228 255 -6.893387734890091 165 206 227 255 -7.018722057342636 163 205 226 255 -7.144056379795181 160 203 226 255 -7.26939070224773 157 202 225 255 -7.3947250247002785 155 200 224 255 -7.520059347152824 152 199 223 255 -7.645393669605369 149 197 222 255 -7.7707279920579175 147 196 222 255 -7.896062314510466 144 194 221 255 -8.021396636963011 141 193 220 255 -8.146730959415557 138 191 219 255 -8.272065281868105 135 189 218 255 -8.397399604320654 133 188 217 255 -8.5227339267732 130 186 216 255 -8.648068249225744 127 184 215 255 -8.773402571678293 124 183 214 255 -8.898736894130842 121 181 213 255 -9.024071216583387 118 179 212 255 -9.149405539035932 116 177 211 255 -9.27473986148848 113 175 210 255 -9.40007418394103 110 174 209 255 -9.525408506393575 107 172 208 255 -9.65074282884612 104 170 207 255 -9.776077151298669 102 168 206 255 -9.901411473751217 99 166 205 255 -10.026745796203762 96 164 204 255 -10.152080118656308 94 162 203 255 -10.277414441108856 91 161 202 255 -10.402748763561405 88 159 201 255 -10.52808308601395 86 157 200 255 -10.653417408466495 83 155 199 255 -10.778751730919044 81 153 198 255 -10.904086053371593 79 151 197 255 -11.029420375824138 76 149 196 255 -11.154754698276683 74 147 195 255 -11.280089020729232 72 146 194 255 -11.40542334318178 70 144 194 255 -11.530757665634326 68 142 193 255 -11.656091988086871 66 140 192 255 -11.78142631053942 64 138 191 255 -11.906760632991968 62 136 190 255 -12.032094955444514 60 135 189 255 -12.157429277897059 58 133 188 255 -12.282763600349607 57 131 187 255 -12.408097922802156 55 129 185 255 -12.533432245254701 54 127 184 255 -12.658766567707247 52 126 183 255 -12.784100890159795 50 124 182 255 -12.909435212612344 49 122 181 255 -13.03476953506489 48 120 180 255 -13.160103857517434 46 118 178 255 -13.285438179969983 45 116 177 255 -13.410772502422532 43 114 176 255 -13.536106824875077 42 113 174 255 -13.661441147327622 41 111 173 255 -13.78677546978017 39 109 171 255 -13.91210979223272 38 107 169 255 -14.037444114685265 37 105 168 255 -14.16277843713781 36 103 166 255 -14.288112759590359 34 101 164 255 -14.413447082042907 33 99 162 255 -14.538781404495452 32 97 160 255 -14.664115726947998 31 95 158 255 -14.789450049400546 30 93 156 255 -14.914784371853095 29 91 154 255 -15.04011869430564 27 89 151 255 -15.165453016758185 26 87 149 255 -15.290787339210734 25 85 147 255 -15.416121661663283 24 83 144 255 -15.541455984115828 23 81 142 255 -15.666790306568373 22 79 139 255 -15.792124629020922 20 77 137 255 diff --git a/tests/test_data/merge/stack_rate.kml b/tests/test_data/merge/stack_rate.kml deleted file mode 100644 index bbe75dfbc..000000000 --- a/tests/test_data/merge/stack_rate.kml +++ /dev/null @@ -1,18 +0,0 @@ - - - - stack_rate.kml - - stack_rate.png - - stack_rate.png - - - -34.1700000 - -34.2300000 - 150.9491667 - 150.9100000 - - - - \ No newline at end of file diff --git a/tests/test_data/merge/stack_rate.png b/tests/test_data/merge/stack_rate.png deleted file mode 100644 index eb662ffdb..000000000 Binary files a/tests/test_data/merge/stack_rate.png and /dev/null differ diff --git a/tests/test_data/merge/stack_rate.png.aux.xml b/tests/test_data/merge/stack_rate.png.aux.xml deleted file mode 100644 index 42bcad758..000000000 --- a/tests/test_data/merge/stack_rate.png.aux.xml +++ /dev/null @@ -1,7 +0,0 @@ - - GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]] - 1.5091000000000000e+02, 8.3333300000000001e-04, 0.0000000000000000e+00, -3.4170000000000002e+01, 0.0000000000000000e+00, -8.3333300000000001e-04 - - PIXEL - - diff --git a/tests/test_data/merge/stack_rate.tif b/tests/test_data/merge/stack_rate.tif deleted file mode 100644 index 7706fe119..000000000 Binary files a/tests/test_data/merge/stack_rate.tif and /dev/null differ diff --git a/tests/test_data/merge/stack_rate.tif.aux.xml b/tests/test_data/merge/stack_rate.tif.aux.xml deleted file mode 100644 index cfae54584..000000000 --- a/tests/test_data/merge/stack_rate.tif.aux.xml +++ /dev/null @@ -1,21 +0,0 @@ - - - STACKED_RATE_MAP - [datetime.date(2006, 6, 19) datetime.date(2006, 8, 28) - datetime.date(2006, 10, 2) datetime.date(2006, 11, 6) - datetime.date(2006, 12, 11) datetime.date(2007, 1, 15) - datetime.date(2007, 2, 19) datetime.date(2007, 3, 26) - datetime.date(2007, 4, 30) datetime.date(2007, 6, 4) - datetime.date(2007, 7, 9) datetime.date(2007, 8, 13) - datetime.date(2007, 9, 17)] - - - - -0 - -8.4672326960268 - -16.042793273926 - 2.3288197678111 - 91.58 - - - diff --git a/tests/test_data/merge/stackrate.kml b/tests/test_data/merge/stackrate.kml deleted file mode 100644 index caf8f5737..000000000 --- a/tests/test_data/merge/stackrate.kml +++ /dev/null @@ -1,18 +0,0 @@ - - - - stackrate.kml - - stackrate.png - - stackrate.png - - - -34.1700000 - -34.2300000 - 150.9491667 - 150.9100000 - - - - diff --git a/tests/test_data/prepifg/tif/geo_060619-061002.tif b/tests/test_data/prepifg/tif/geo_060619-061002_unw.tif similarity index 100% rename from tests/test_data/prepifg/tif/geo_060619-061002.tif rename to tests/test_data/prepifg/tif/geo_060619-061002_unw.tif diff --git a/tests/test_data/prepifg/tif/geo_070326-070917.tif b/tests/test_data/prepifg/tif/geo_070326-070917_unw.tif similarity index 100% rename from tests/test_data/prepifg/tif/geo_070326-070917.tif rename to tests/test_data/prepifg/tif/geo_070326-070917_unw.tif diff --git a/tests/test_data/small_test/conf/pyrate_gamma_test.conf b/tests/test_data/small_test/conf/pyrate_gamma_test.conf index 1d28949a0..e5f431490 100644 --- a/tests/test_data/small_test/conf/pyrate_gamma_test.conf +++ b/tests/test_data/small_test/conf/pyrate_gamma_test.conf @@ -39,7 +39,7 @@ nan_conversion: 1 # gamma prepifg runs in parallel in single machine if parallel != 0 # parallel = 1, stacking/timeseries computation is done in parallel # parallel = 0, stacking/timeseries computation is done in serial -parallel: 1 +parallel: 0 processes: 4 #------------------------------------ @@ -62,8 +62,8 @@ ifgylast: -34.22 # refnx/y: number of search grid points in x/y direction # refchipsize: chip 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 @@ -158,5 +158,9 @@ nsig: 3 pthr: 5 maxsig: 2 +# optionally supply rows and cols for tiles used during process and merge step rows: 3 cols: 2 + +# optionally save stack numpy files during merge step +savenpy: 1 diff --git a/tests/test_data/small_test/dem/roipac_test_trimmed_4rlks_3cr.tif b/tests/test_data/small_test/dem/roipac_test_trimmed_4rlks_3cr.tif index feceffb49..fc1eb0053 100644 Binary files a/tests/test_data/small_test/dem/roipac_test_trimmed_4rlks_3cr.tif and b/tests/test_data/small_test/dem/roipac_test_trimmed_4rlks_3cr.tif differ diff --git a/tests/test_gamma.py b/tests/test_gamma.py index c4e5ede83..907383267 100644 --- a/tests/test_gamma.py +++ b/tests/test_gamma.py @@ -290,7 +290,7 @@ def test_fail_bad_date_order(self): self.assertRaises(self.err, gamma.combine_headers, H1, H0, self.dh) -glob_prefix = "*utm_unw_1rlks_1cr.tif" +glob_prefix = "*utm_unw_ifg_1rlks_1cr.tif" @pytest.fixture(scope='module') diff --git a/tests/test_gamma_vs_roipac.py b/tests/test_gamma_vs_roipac.py index e83f64577..a0eb6295e 100644 --- a/tests/test_gamma_vs_roipac.py +++ b/tests/test_gamma_vs_roipac.py @@ -42,7 +42,7 @@ def test_files_are_same(tempdir, get_config): gamma_params = __workflow(gamma_params, gamma_tdir) # conv2tif output equal - __assert_same_files_produced(roipac_params[cf.OUT_DIR], gamma_params[cf.OUT_DIR], "*_unw.tif", 17) + __assert_same_files_produced(roipac_params[cf.OUT_DIR], gamma_params[cf.OUT_DIR], "*_unw_ifg.tif", 17) # prepifg output equal __assert_same_files_produced(roipac_params[cf.OUT_DIR], gamma_params[cf.OUT_DIR], diff --git a/tests/test_merge.py b/tests/test_merge.py index 1c30b8334..ca975b97d 100644 --- a/tests/test_merge.py +++ b/tests/test_merge.py @@ -18,37 +18,49 @@ This Python module contains tests for the Merge step of PyRate. """ import os -import unittest +from subprocess import check_call +import itertools +import pytest from pathlib import Path -from pyrate.merge import create_png_from_tif -from tests.common import TESTDIR +from pyrate.merge import _create_png_from_tif +from pyrate.core import config as cf +from pyrate.merge import _merge_stack +from pyrate.configuration import Configuration +from tests.common import manipulate_test_conf +@pytest.fixture +def create_stack_output(tempdir, gamma_conf): + tdir = Path(tempdir()) + params = manipulate_test_conf(gamma_conf, tdir) + output_conf_file = tdir.joinpath('conf.cfg') + output_conf = tdir.joinpath(output_conf_file) + cf.write_config_file(params=params, output_conf_file=output_conf) + check_call(f"pyrate conv2tif -f {output_conf}", shell=True) + check_call(f"pyrate prepifg -f {output_conf}", shell=True) + check_call(f"pyrate process -f {output_conf}", shell=True) -class MergingTest(unittest.TestCase): + params = Configuration(output_conf).__dict__ + return params, tdir - def test_png_creation(self): - output_folder_path = Path(TESTDIR).joinpath("test_data", "merge") - create_png_from_tif(output_folder_path) +@pytest.mark.slow +def test_png_creation(create_stack_output): + params, tdir = create_stack_output - # check if color map is created - output_color_map_path = os.path.join(output_folder_path, "colourmap.txt") - if not os.path.isfile(output_color_map_path): - self.assertTrue(False, "Output color map file not found at: " + output_color_map_path) + output_folder_path = params[cf.OUT_DIR] - # check if png is created - output_image_path = os.path.join(output_folder_path, "stack_rate.png") - if not os.path.isfile(output_image_path): - self.assertTrue(False, "Output png file not found at: " + output_image_path) + rows, cols = params["rows"], params["cols"] + _merge_stack(rows, cols, params) + _create_png_from_tif(output_folder_path) - # check if kml is created - output_kml_path = os.path.join(output_folder_path, "stack_rate.kml") - if not os.path.isfile(output_kml_path): - self.assertTrue(False, "Output kml file not found at: " + output_kml_path) + # check if color map is created + for out_type in ['rate', 'error']: + output_color_map_path = os.path.join(output_folder_path, f"colourmap_{out_type}.txt") + assert Path(output_color_map_path).exists(), "Output color map file not found at: " + output_color_map_path - self.assertTrue(True) - - -if __name__ == '__main__': - unittest.main() + # check if merged files are created + for _type, output_type in itertools.product(["stack_rate", "stack_error"], ['.tif', '.png', '.kml']): + output_image_path = os.path.join(output_folder_path, _type + output_type) + print(f"checking {output_image_path}") + assert Path(output_image_path).exists(), f"Output {output_type} file not found at {output_image_path}" diff --git a/tests/test_mpi.py b/tests/test_mpi.py index c0ccc2830..d61381186 100644 --- a/tests/test_mpi.py +++ b/tests/test_mpi.py @@ -22,9 +22,6 @@ import numpy as np import os from pathlib import Path - -import pyrate.core.orbital -import pyrate.core.shared from pyrate import process, prepifg, conv2tif, configuration from pyrate.core import mpiops, config as cf from tests import common @@ -33,29 +30,33 @@ def test_vcm_legacy_vs_mpi(mpisync, tempdir, roipac_or_gamma_conf): + params = configuration.Configuration(roipac_or_gamma_conf).__dict__ LEGACY_VCM_DIR = os.path.join(SML_TEST_DIR, 'vcm') legacy_vcm = np.genfromtxt(os.path.join(LEGACY_VCM_DIR, 'vcmt.csv'), delimiter=',') tmpdir = Path(mpiops.run_once(tempdir)) mpiops.run_once(common.copytree, params[cf.OBS_DIR], tmpdir) params[cf.OUT_DIR] = tmpdir.joinpath('out') - params[cf.PARALLEL] = False - xlks, ylks, crop = cf.transform_params(params) - base_unw_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], params[cf.OBS_DIR]) - # dest_paths are tifs that have been geotif converted and multilooked - dest_paths = cf.get_dest_paths(base_unw_paths, crop, params, xlks) + params[cf.PARALLEL] = 0 + output_conf = Path(tmpdir).joinpath('conf.cfg') + cf.write_config_file(params=params, output_conf_file=output_conf) + params = configuration.Configuration(output_conf).__dict__ + dest_paths = [p.sampled_path for p in params[cf.INTERFEROGRAM_FILES]] # run conv2tif and prepifg, create the dest_paths files conv2tif.main(params) + params[cf.INTERFEROGRAM_FILES].pop() prepifg.main(params) - - tiles = pyrate.core.shared.get_tiles(dest_paths[0], rows=1, cols=1) - preread_ifgs = process._create_ifg_dict(dest_paths, params=params, tiles=tiles) + params[cf.INTERFEROGRAM_FILES].pop() + preread_ifgs = process._create_ifg_dict(dest_paths, params=params) refpx, refpy = process._ref_pixel_calc(dest_paths, params) - process._orb_fit_calc(dest_paths, params) + process._orb_fit_calc(params[cf.INTERFEROGRAM_FILES], params) process._ref_phase_estimation(dest_paths, params, refpx, refpy) maxvar, vcmt = process._maxvar_vcm_calc(dest_paths, params, preread_ifgs) + + # phase data after ref pixel has changed due to commit bf2f7ebd + # Legacy tests won't match anymore np.testing.assert_array_almost_equal(maxvar, legacy_maxvar, decimal=4) np.testing.assert_array_almost_equal(legacy_vcm, vcmt, decimal=3) mpiops.run_once(shutil.rmtree, tmpdir) diff --git a/tests/test_mpi_vs_multiprocess_vs_single_process.py b/tests/test_mpi_vs_multiprocess_vs_single_process.py index 5ad20cbaf..ff782cb12 100644 --- a/tests/test_mpi_vs_multiprocess_vs_single_process.py +++ b/tests/test_mpi_vs_multiprocess_vs_single_process.py @@ -5,13 +5,17 @@ from subprocess import check_call, check_output, CalledProcessError import numpy as np from pyrate.core import config as cf -from tests.common import assert_same_files_produced, assert_two_dirs_equal, manipulate_test_conf +from tests.common import ( + assert_same_files_produced, + assert_two_dirs_equal, + manipulate_test_conf, + TRAVIS, + PYTHON3P6, + PYTHON3P7, + PYTHON3P8, + GDAL_VERSION +) -TRAVIS = True if 'TRAVIS' in os.environ else False -PYTHON3P6 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.6') else False -PYTHON3P7 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.7') else False -PYTHON3P8 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.8') else False -GDAL_VERSION = check_output(["gdal-config", "--version"]).decode(encoding="utf-8").split('\n')[0] # python3.7 and gdal3.0.4 REGRESSION = PYTHON3P7 and (GDAL_VERSION == '3.0.4') # python3.7 and gdal3.0.2 @@ -23,7 +27,7 @@ def parallel(request): return request.param -@pytest.fixture(params=[1, 2, 3, 4]) +@pytest.fixture(params=[1, 2, 4]) def local_crop(request): return request.param @@ -50,6 +54,7 @@ def modify_params(conf_file, parallel_vs_serial, output_conf_file): params[cf.ORBITAL_FIT_DEGREE] = orbfit_degrees params[cf.REF_EST_METHOD] = ref_est_method params["rows"], params["cols"] = 3, 2 + params["savenpy"] = 1 params["tiles"] = params["rows"] * params["cols"] print(params) @@ -65,7 +70,7 @@ def modify_params(conf_file, parallel_vs_serial, output_conf_file): @pytest.mark.skipif(REGRESSION or PYTHON3P8 or PYTHON3P6, reason="Skip if not python3.7 and gdal=3.0.4") def test_pipeline_parallel_vs_mpi(modified_config, gamma_conf): - if TRAVIS and np.random.randint(0, 1000) > 149: # skip 85% of tests randomly + if np.random.randint(0, 1000) > 149: # skip 85% of tests randomly pytest.skip("Randomly skipping as part of 85 percent") print("\n\n") @@ -86,32 +91,28 @@ def test_pipeline_parallel_vs_mpi(modified_config, gamma_conf): mr_conf, params_m = modified_config(gamma_conf, 1, 'multiprocess_conf.conf') - check_call(f"pyrate conv2tif -f {mr_conf}", shell=True) - check_call(f"pyrate prepifg -f {mr_conf}", shell=True) - check_call(f"pyrate process -f {mr_conf}", shell=True) - check_call(f"pyrate merge -f {mr_conf}", shell=True) + check_call(f"pyrate workflow -f {mr_conf}", shell=True) sr_conf, params_s = modified_config(gamma_conf, 0, 'singleprocess_conf.conf') - check_call(f"pyrate conv2tif -f {sr_conf}", shell=True) - check_call(f"pyrate prepifg -f {sr_conf}", shell=True) - check_call(f"pyrate process -f {sr_conf}", shell=True) - check_call(f"pyrate merge -f {sr_conf}", shell=True) + check_call(f"pyrate workflow -f {sr_conf}", shell=True) # convert2tif tests, 17 interferograms - assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "*_unw.tif", 17) + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "*_unw_ifg.tif", 17) # if coherence masking, comprare coh files were converted if params[cf.COH_MASK]: - assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "*_utm.tif", 17) + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "*_coh.tif", 17) print("coherence files compared") - - # prepifg + process steps that overwrite tifs test - - # 17 ifgs + 1 dem - assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], + # 17 ifgs + 1 dem + 17 mlooked coh files + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], + f"*{params[cf.IFG_CROP_OPT]}cr.tif", 35) + else: + # 17 ifgs + 1 dem + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 18) + # prepifg + process steps that overwrite tifs test # ifg phase checking in the previous step checks the process pipeline upto APS correction # 2 x because of aps files @@ -130,8 +131,11 @@ def test_pipeline_parallel_vs_mpi(modified_config, gamma_conf): # compare merge step assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "stack*.tif", 3) + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "stack*.kml", 2) + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "stack*.png", 2) assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "stack*.npy", 3) - assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "tscuml*.tif") + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "tscuml*.tif", 12) + assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "tsincr*.tif", 12) print("==========================xxx===========================") @@ -146,7 +150,7 @@ def coh_mask(request): @pytest.fixture() -def modified_config_short(tempdir, local_crop, get_lks, coh_mask): +def modified_config_short(tempdir, local_crop, get_lks, coh_mask, ref_pixel): orbfit_lks = 1 orbfit_method = 1 orbfit_degrees = 1 @@ -161,6 +165,7 @@ def modify_params(conf_file, parallel, output_conf_file, largetifs): params[cf.APSEST] = 1 params[cf.LARGE_TIFS] = largetifs params[cf.IFG_LKSX], params[cf.IFG_LKSY] = get_lks, get_lks + params[cf.REFX], params[cf.REFY] = ref_pixel params[cf.REFNX], params[cf.REFNY] = 4, 4 params[cf.IFG_CROP_OPT] = local_crop @@ -170,6 +175,7 @@ def modify_params(conf_file, parallel, output_conf_file, largetifs): params[cf.ORBITAL_FIT_DEGREE] = orbfit_degrees params[cf.REF_EST_METHOD] = ref_est_method params["rows"], params["cols"] = 3, 2 + params["savenpy"] = 1 params["tiles"] = params["rows"] * params["cols"] print(params) @@ -214,7 +220,7 @@ def test_stack_and_ts_mpi_vs_parallel_vs_serial(modified_config_short, gamma_con 3. Doing 1 and 2 means we have checked single vs parallel python multiprocess pipelines 4. This also checks the entire pipeline using largetifs (new prepifg) vs old perpifg (python based) """ - if TRAVIS and np.random.randint(0, 1000) > 399: # skip 60% of tests randomly + if np.random.randint(0, 1000) > 399: # skip 60% of tests randomly pytest.skip("Randomly skipping as part of 60 percent") print("\n\n") @@ -225,24 +231,22 @@ def test_stack_and_ts_mpi_vs_parallel_vs_serial(modified_config_short, gamma_con sr_conf, params_p = modified_config_short(gamma_conf, parallel, 'parallel_conf.conf', 0) - check_call(f"pyrate conv2tif -f {sr_conf}", shell=True) - check_call(f"pyrate prepifg -f {sr_conf}", shell=True) - check_call(f"pyrate process -f {sr_conf}", shell=True) - check_call(f"pyrate merge -f {sr_conf}", shell=True) - + check_call(f"pyrate workflow -f {sr_conf}", shell=True) # convert2tif tests, 17 interferograms - assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_unw.tif", 17) + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_unw_ifg.tif", 17) # if coherence masking, compare coh files were converted if params[cf.COH_MASK]: - assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_utm.tif", 17) + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_coh.tif", 17) print("coherence files compared") # prepifg + process steps that overwrite tifs test - - # 17 ifgs + 1 dem - assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 18) + # 17 mlooked ifgs + 1 dem + 17 mlooked coherence files + if params[cf.COH_MASK]: + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 35) + else: + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 18) # ifg phase checking in the previous step checks the process pipeline upto APS correction assert_two_dirs_equal(params[cf.TMPDIR], params_p[cf.TMPDIR], "tsincr_*.npy", params['tiles'] * 2) @@ -253,6 +257,8 @@ def test_stack_and_ts_mpi_vs_parallel_vs_serial(modified_config_short, gamma_con # compare merge step assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "stack*.tif", 3) + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "stack*.kml", 2) + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "stack*.png", 2) assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "stack*.npy", 3) assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "tscuml*.tif") diff --git a/tests/test_orbital.py b/tests/test_orbital.py index 427360168..5676a7647 100644 --- a/tests/test_orbital.py +++ b/tests/test_orbital.py @@ -31,7 +31,7 @@ from numpy.testing import assert_array_equal, assert_array_almost_equal from scipy.linalg import lstsq -from .common import small5_mock_ifgs, MockIfg +from tests.common import small5_mock_ifgs, MockIfg from pyrate.core import algorithm, config as cf from pyrate.core.orbital import INDEPENDENT_METHOD, NETWORK_METHOD, PLANAR, \ QUADRATIC, PART_CUBIC @@ -40,6 +40,9 @@ from pyrate.core.orbital import _get_num_params, remove_orbital_error from pyrate.core.shared import Ifg from pyrate.core.shared import nanmedian +from pyrate.core import roipac +from pyrate.configuration import Configuration +from tests import common from tests.common import TEST_CONF_ROIPAC, IFMS16 from tests.common import SML_TEST_LEGACY_ORBITAL_DIR from tests.common import SML_TEST_TIF, small_data_setup @@ -697,8 +700,7 @@ def setUp(self): self.params[cf.PARALLEL] = False data_paths = [os.path.join(SML_TEST_TIF, p) for p in IFMS16] - self.ifg_paths = [os.path.join(self.BASE_DIR, os.path.basename(d)) - for d in data_paths] + self.ifg_paths = [os.path.join(self.BASE_DIR, os.path.basename(d)) for d in data_paths] for d in data_paths: shutil.copy(d, os.path.join(self.BASE_DIR, os.path.basename(d))) @@ -708,8 +710,12 @@ def tearDown(self): def test_orbital_correction_legacy_equality(self): from pyrate import process + from pyrate.configuration import MultiplePaths - process._orb_fit_calc(self.ifg_paths, self.params) + multi_paths = [MultiplePaths(self.BASE_DIR, p) for p in self.ifg_paths] + for m in multi_paths: # cheat + m.sampled_path = m.converted_path + process._orb_fit_calc(multi_paths, self.params) onlyfiles = [f for f in os.listdir(SML_TEST_LEGACY_ORBITAL_DIR) if os.path.isfile(os.path.join(SML_TEST_LEGACY_ORBITAL_DIR, f)) @@ -754,12 +760,11 @@ class LegacyComparisonTestsOrbfitMethod2(unittest.TestCase): """ def setUp(self): self.BASE_DIR = tempfile.mkdtemp() - self.params = cf.get_config_params(TEST_CONF_ROIPAC) # change to orbital error correction method 2 + self.params = Configuration(common.TEST_CONF_ROIPAC).__dict__ self.params[cf.ORBITAL_FIT_METHOD] = NETWORK_METHOD self.params[cf.ORBITAL_FIT_LOOKS_X] = 1 self.params[cf.ORBITAL_FIT_LOOKS_Y] = 1 - data_paths = [os.path.join(SML_TEST_TIF, p) for p in small_ifg_file_list()] self.new_data_paths = [os.path.join(self.BASE_DIR, os.path.basename(d)) for d in data_paths] for d in data_paths: @@ -768,6 +773,7 @@ def setUp(self): os.chmod(d_copy, 0o660) self.ifgs = small_data_setup(datafiles=self.new_data_paths) + self.headers = [roipac.roipac_header(i.data_path, self.params) for i in self.ifgs] for i in self.ifgs: if not i.is_open: @@ -785,7 +791,7 @@ def tearDown(self): shutil.rmtree(self.BASE_DIR) def test_orbital_correction_legacy_equality_orbfit_method_2(self): - remove_orbital_error(self.ifgs, self.params) + remove_orbital_error(self.ifgs, self.params, headers=self.headers) onlyfiles = [f for f in os.listdir(SML_TEST_LEGACY_ORBITAL_DIR) if os.path.isfile(os.path.join(SML_TEST_LEGACY_ORBITAL_DIR, f)) @@ -818,7 +824,7 @@ def test_orbital_error_method2_dummy(self): self.params[cf.ORBITAL_FIT_LOOKS_X] = 2 self.params[cf.ORBITAL_FIT_LOOKS_Y] = 2 - remove_orbital_error(self.ifgs, self.params) + remove_orbital_error(self.ifgs, self.params, self.headers) onlyfiles = [f for f in os.listdir(SML_TEST_LEGACY_ORBITAL_DIR) if os.path.isfile(os.path.join(SML_TEST_LEGACY_ORBITAL_DIR, f)) diff --git a/tests/test_prepifg.py b/tests/test_prepifg.py index d79782a82..b666dec3c 100644 --- a/tests/test_prepifg.py +++ b/tests/test_prepifg.py @@ -23,6 +23,8 @@ import unittest from math import floor from os.path import exists, join +from pathlib import Path +from subprocess import check_call import numpy as np from numpy import isnan, nanmax, nanmin, nanmean, ones, nan, reshape, sum as npsum @@ -35,12 +37,16 @@ from pyrate.core.shared import Ifg, DEM from pyrate.core.prepifg_helper import CUSTOM_CROP, MAXIMUM_CROP, MINIMUM_CROP, \ ALREADY_SAME_SIZE +from pyrate.core import roipac from pyrate.core.prepifg_helper import prepare_ifgs, _resample, PreprocessError, CustomExts +from pyrate.core import ifgconstants as ifc +from pyrate.configuration import Configuration +from pyrate import conv2tif, prepifg + from tests import common from tests.common import SML_TEST_LEGACY_PREPIFG_DIR -from tests.common import PREP_TEST_TIF, SML_TEST_DEM_DIR -from tests.common import SML_TEST_DEM_TIF -from pyrate import conv2tif, prepifg +from tests.common import PREP_TEST_TIF, SML_TEST_DEM_DIR, PREP_TEST_OBS +from tests.common import SML_TEST_DEM_TIF, SML_TEST_DEM_HDR, manipulate_test_conf gdal.UseExceptions() DUMMY_SECTION_NAME = 'pyrate' @@ -49,14 +55,99 @@ sys.exit("ERROR: Missing 'prepifg' dir for unittests\n") +def test_prepifg_treat_inputs_read_only(gamma_conf, tempdir, coh_mask): + tdir = Path(tempdir()) + params = common.manipulate_test_conf(gamma_conf, tdir) + params[cf.COH_MASK] = coh_mask + output_conf = tdir.joinpath('conf.cfg') + cf.write_config_file(params=params, output_conf_file=output_conf) + check_call(f"mpirun -n 3 pyrate conv2tif -f {output_conf}", shell=True) + tifs = list(Path(params[cf.OUT_DIR]).glob('*_unw_ifg.tif')) + assert len(tifs) == 17 + + check_call(f"mpirun -n 3 pyrate prepifg -f {output_conf}", shell=True) + cropped = list(Path(params[cf.OUT_DIR]).glob('*cr.tif')) + + if coh_mask: # 17 + 1 dem + 17 coh files + assert len(cropped) == 35 + else: # 17 + 1 dem + assert len(cropped) == 18 + # check all tifs from conv2tif are still readonly + for t in tifs: + assert t.stat().st_mode == 33060 + + +def test_prepifg_file_types(tempdir, gamma_conf, coh_mask): + tdir = Path(tempdir()) + params = manipulate_test_conf(gamma_conf, tdir) + params[cf.COH_MASK] = coh_mask + params[cf.PARALLEL] = 0 + output_conf_file = 'conf.conf' + output_conf = tdir.joinpath(output_conf_file) + cf.write_config_file(params=params, output_conf_file=output_conf) + params_s = Configuration(output_conf).__dict__ + conv2tif.main(params_s) + # reread params from config + params_s = Configuration(output_conf).__dict__ + prepifg.main(params_s) + ifg_files = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_ifg.tif')) + assert len(ifg_files) == 17 + mlooked_files = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_ifg_1rlks_1cr.tif')) + assert len(mlooked_files) == 17 + coh_files = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_cc_coh.tif')) + mlooked_coh_files = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_cc_coh_1rlks_1cr.tif')) + if coh_mask: + assert len(coh_files) == 17 + assert len(mlooked_coh_files) == 17 + dem_file = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_dem.tif'))[0] + mlooked_dem_file = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_dem_1rlks_1cr.tif'))[0] + import itertools + + # assert coherence and ifgs have correct metadata + for i in itertools.chain(*[ifg_files, mlooked_files, coh_files, mlooked_coh_files]): + ifg = Ifg(i) + ifg.open() + md = ifg.meta_data + if i.name.endswith('_ifg.tif'): + assert md[ifc.DATA_TYPE] == ifc.ORIG + continue + if i.name.endswith('_coh.tif'): + assert md[ifc.DATA_TYPE] == ifc.COH + continue + if i.name.endswith('_cc_coh_1rlks_1cr.tif'): + assert md[ifc.DATA_TYPE] == ifc.MULTILOOKED_COH + continue + if i.name.endswith('_ifg_1rlks_1cr.tif'): + if coh_mask: + assert md[ifc.DATA_TYPE] == ifc.COHERENCE + else: + assert md[ifc.DATA_TYPE] == ifc.MULTILOOKED + continue + + # assert dem has correct metadata + dem = DEM(dem_file.as_posix()) + dem.open() + md = dem.dataset.GetMetadata() + assert md[ifc.DATA_TYPE] == ifc.DEM + + dem = DEM(mlooked_dem_file.as_posix()) + dem.open() + md = dem.dataset.GetMetadata() + assert md[ifc.DATA_TYPE] == ifc.MLOOKED_DEM + shutil.rmtree(tdir) + + # convenience ifg creation funcs def diff_exts_ifgs(): """Returns pair of test Ifgs with different extents""" - bases = ['geo_060619-061002.tif', 'geo_070326-070917.tif'] + bases = ['geo_060619-061002_unw.tif', 'geo_070326-070917_unw.tif'] + headers = ['geo_060619-061002.unw.rsc', 'geo_070326-070917.unw.rsc'] random_dir = tempfile.mkdtemp() - for p in bases: + for p, h in zip(bases, headers): shutil.copy(src=os.path.join(PREP_TEST_TIF, p), dst=os.path.join(random_dir, p)) + shutil.copy(src=os.path.join(PREP_TEST_OBS, h), + dst=os.path.join(random_dir, h)) return [Ifg(join(random_dir, p)) for p in bases], random_dir @@ -110,14 +201,17 @@ def setUp(self): self.ys = -self.xs self.ifgs, self.random_dir = diff_exts_ifgs() self.ifg_paths = [i.data_path for i in self.ifgs] - paths = ["geo_060619-061002_1rlks_1cr.tif", - "geo_060619-061002_1rlks_2cr.tif", - "geo_060619-061002_1rlks_3cr.tif", - "geo_060619-061002_4rlks_3cr.tif", - "geo_070326-070917_1rlks_1cr.tif", - "geo_070326-070917_1rlks_2cr.tif", - "geo_070326-070917_1rlks_3cr.tif", - "geo_070326-070917_4rlks_3cr.tif"] + + params = Configuration(common.TEST_CONF_ROIPAC).__dict__ + self.headers = [roipac.roipac_header(i.data_path, params) for i in self.ifgs] + paths = ["geo_060619-061002_unw_1rlks_1cr.tif", + "geo_060619-061002_unw_1rlks_2cr.tif", + "geo_060619-061002_unw_1rlks_3cr.tif", + "geo_060619-061002_unw_4rlks_3cr.tif", + "geo_070326-070917_unw_1rlks_1cr.tif", + "geo_070326-070917_unw_1rlks_2cr.tif", + "geo_070326-070917_unw_1rlks_3cr.tif", + "geo_070326-070917_unw_4rlks_3cr.tif"] self.exp_files = [join(self.random_dir, p) for p in paths] @staticmethod @@ -164,7 +258,7 @@ def assert_projection_equal(self, files): def test_multilooked_projection_same_as_geotiff(self): xlooks = ylooks = 1 - prepare_ifgs(self.ifg_paths, MAXIMUM_CROP, xlooks, ylooks) + prepare_ifgs(self.ifg_paths, MAXIMUM_CROP, xlooks, ylooks, headers=self.headers) mlooked_paths = [mlooked_path(f, crop_out=MAXIMUM_CROP, looks=xlooks) for f in self.ifg_paths] self.assert_projection_equal(self.ifg_paths + mlooked_paths) @@ -172,7 +266,7 @@ def test_multilooked_projection_same_as_geotiff(self): def test_default_max_extents(self): """Test ifgcropopt=2 crops datasets to max bounding box extents.""" xlooks = ylooks = 1 - prepare_ifgs(self.ifg_paths, MAXIMUM_CROP, xlooks, ylooks) + prepare_ifgs(self.ifg_paths, MAXIMUM_CROP, xlooks, ylooks, self.headers) for f in [self.exp_files[1], self.exp_files[5]]: self.assertTrue(exists(f), msg="Output files not created") @@ -196,7 +290,7 @@ def test_default_max_extents(self): def test_min_extents(self): """Test ifgcropopt=1 crops datasets to min extents.""" xlooks = ylooks = 1 - prepare_ifgs(self.ifg_paths, MINIMUM_CROP, xlooks, ylooks) + prepare_ifgs(self.ifg_paths, MINIMUM_CROP, xlooks, ylooks, headers=self.headers) ifg = Ifg(self.exp_files[0]) ifg.open() @@ -204,8 +298,7 @@ def test_min_extents(self): # NB: also verifies gdalwarp correctly copies geotransform across # NB: expected data copied from gdalinfo output gt = ifg.dataset.GetGeoTransform() - exp_gt = (150.911666666, 0.000833333, 0, - -34.172499999, 0, -0.000833333) + exp_gt = (150.911666666, 0.000833333, 0, -34.172499999, 0, -0.000833333) for i, j in zip(gt, exp_gt): self.assertAlmostEqual(i, j) self.assert_geotransform_equal([self.exp_files[0], self.exp_files[4]]) @@ -217,8 +310,7 @@ def test_min_extents(self): def test_custom_extents(self): xlooks = ylooks = 1 cext = self._custom_extents_tuple() - prepare_ifgs(self.ifg_paths, CUSTOM_CROP, xlooks, ylooks, - user_exts=cext) + prepare_ifgs(self.ifg_paths, CUSTOM_CROP, xlooks, ylooks, headers=self.headers, user_exts=cext) ifg = Ifg(self.exp_files[2]) ifg.open() @@ -242,10 +334,10 @@ def test_exception_without_all_4_crop_parameters(self): for i in [None, '']: cext = (150.92, -34.18, 150.94, i) self.assertRaises(PreprocessError, prepare_ifgs, self.ifg_paths, - CUSTOM_CROP, xlooks, ylooks, user_exts=cext) + CUSTOM_CROP, xlooks, ylooks, self.headers, user_exts=cext) # three parameters provided self.assertRaises(PreprocessError, prepare_ifgs, self.ifg_paths, - CUSTOM_CROP, xlooks, ylooks, + CUSTOM_CROP, xlooks, ylooks, self.headers, user_exts=(150.92, -34.18, 150.94)) # close ifgs for i in self.ifgs: @@ -264,7 +356,7 @@ def test_custom_extents_misalignment(self): self.assertRaises(PreprocessError, prepare_ifgs, self.ifg_paths, CUSTOM_CROP, - xlooks, ylooks, user_exts=cext) + xlooks, ylooks, user_exts=cext, headers=self.headers) # close ifgs for i in self.ifgs: i.close() @@ -272,7 +364,7 @@ def test_custom_extents_misalignment(self): def test_nodata(self): """Verify NODATA value copied correctly (amplitude band not copied)""" xlooks = ylooks = 1 - prepare_ifgs(self.ifg_paths, MINIMUM_CROP, xlooks, ylooks) + prepare_ifgs(self.ifg_paths, MINIMUM_CROP, xlooks, ylooks, self.headers) for ex in [self.exp_files[0], self.exp_files[4]]: ifg = Ifg(ex) @@ -286,7 +378,7 @@ def test_nodata(self): def test_nans(self): """Verify NaNs replace 0 in the multilooked phase band""" xlooks = ylooks = 1 - prepare_ifgs(self.ifg_paths, MINIMUM_CROP, xlooks, ylooks) + prepare_ifgs(self.ifg_paths, MINIMUM_CROP, xlooks, ylooks, self.headers) for ex in [self.exp_files[0], self.exp_files[4]]: ifg = Ifg(ex) @@ -307,9 +399,11 @@ def test_multilook(self): scale = 4 # assumes square cells self.ifgs.append(DEM(SML_TEST_DEM_TIF)) self.ifg_paths = [i.data_path for i in self.ifgs] + # append the dem header + self.headers.append(SML_TEST_DEM_HDR) cext = self._custom_extents_tuple() xlooks = ylooks = scale - prepare_ifgs(self.ifg_paths, CUSTOM_CROP, xlooks, ylooks, thresh=1.0, user_exts=cext) + prepare_ifgs(self.ifg_paths, CUSTOM_CROP, xlooks, ylooks, thresh=1.0, user_exts=cext, headers=self.headers) for n, ipath in enumerate([self.exp_files[3], self.exp_files[7]]): i = Ifg(ipath) @@ -357,10 +451,12 @@ def test_output_datatype(self): scale = 4 # assumes square cells self.ifgs.append(DEM(SML_TEST_DEM_TIF)) self.ifg_paths = [i.data_path for i in self.ifgs] + [SML_TEST_DEM_TIF] + self.headers.append(SML_TEST_DEM_HDR) + cext = self._custom_extents_tuple() xlooks = ylooks = scale prepare_ifgs(self.ifg_paths, CUSTOM_CROP, xlooks, ylooks, - thresh=1.0, user_exts=cext) + thresh=1.0, user_exts=cext, headers=self.headers) for i in self.ifg_paths: mlooked_ifg = mlooked_path(i, xlooks, CUSTOM_CROP) @@ -377,10 +473,10 @@ def test_invalid_looks(self): values = [0, -1, -10, -100000.6, ""] for v in values: self.assertRaises(PreprocessError, prepare_ifgs, self.ifg_paths, - CUSTOM_CROP, xlooks=v, ylooks=1) + CUSTOM_CROP, xlooks=v, ylooks=1, headers=self.headers) self.assertRaises(PreprocessError, prepare_ifgs, self.ifg_paths, - CUSTOM_CROP, xlooks=1, ylooks=v) + CUSTOM_CROP, xlooks=1, ylooks=v, headers=self.headers) class ThresholdTests(unittest.TestCase): @@ -426,9 +522,22 @@ class SameSizeTests(unittest.TestCase): """Tests aspects of the prepifg.py script, such as resampling.""" def __init__(self, *args, **kwargs): + import datetime super(SameSizeTests, self).__init__(*args, **kwargs) self.xs = 0.000833333 self.ys = -self.xs + self.headers = [ + {'NCOLS': 47, 'NROWS': 72, 'LAT': -34.17, 'LONG': 150.91, 'X_STEP': 0.000833333, 'Y_STEP': -0.000833333, + 'WAVELENGTH_METRES': 0.0562356424, 'MASTER_DATE': datetime.date(2007, 3, 26), + 'SLAVE_DATE': datetime.date(2007, 9, 17), 'TIME_SPAN_YEAR': 0.4791238877481177, + 'DATA_UNITS': 'RADIANS', 'INSAR_PROCESSOR': 'ROIPAC', 'X_LAST': 150.94916665099998, + 'Y_LAST': -34.229999976, 'DATUM': 'WGS84', 'DATA_TYPE': 'ORIGINAL_IFG'}, + {'NCOLS': 47, 'NROWS': 72, 'LAT': -34.17, 'LONG': 150.91, 'X_STEP': 0.000833333, 'Y_STEP': -0.000833333, + 'WAVELENGTH_METRES': 0.0562356424, 'MASTER_DATE': datetime.date(2007, 3, 26), + 'SLAVE_DATE': datetime.date(2007, 9, 17), 'TIME_SPAN_YEAR': 0.4791238877481177, + 'DATA_UNITS': 'RADIANS', 'INSAR_PROCESSOR': 'ROIPAC', 'X_LAST': 150.94916665099998, + 'Y_LAST': -34.229999976, 'DATUM': 'WGS84', 'DATA_TYPE': 'ORIGINAL_IFG'} + ] # TODO: check output files for same extents? # TODO: make prepifg dir readonly to test output to temp dir @@ -437,15 +546,14 @@ def test_already_same_size(self): # should do nothing as layers are same size & no multilooking required ifgs = same_exts_ifgs() ifg_data_paths = [d.data_path for d in ifgs] - res_tup = prepare_ifgs(ifg_data_paths, ALREADY_SAME_SIZE, 1, 1) + res_tup = prepare_ifgs(ifg_data_paths, ALREADY_SAME_SIZE, 1, 1, self.headers) res = [r[1] for r in res_tup] self.assertTrue(all(res)) def test_already_same_size_mismatch(self): ifgs, random_dir = diff_exts_ifgs() ifg_data_paths = [d.data_path for d in ifgs] - self.assertRaises(PreprocessError, prepare_ifgs, - ifg_data_paths, ALREADY_SAME_SIZE, 1, 1) + self.assertRaises(PreprocessError, prepare_ifgs, ifg_data_paths, ALREADY_SAME_SIZE, 1, 1, self.headers) for i in ifgs: i.close() shutil.rmtree(random_dir) @@ -455,10 +563,9 @@ def test_same_size_multilooking(self): ifgs = same_exts_ifgs() ifg_data_paths = [d.data_path for d in ifgs] xlooks = ylooks = 2 - prepare_ifgs(ifg_data_paths, ALREADY_SAME_SIZE, xlooks, ylooks) + prepare_ifgs(ifg_data_paths, ALREADY_SAME_SIZE, xlooks, ylooks, self.headers) - looks_paths = [mlooked_path(d, looks=xlooks, crop_out=ALREADY_SAME_SIZE) - for d in ifg_data_paths] + looks_paths = [mlooked_path(d, looks=xlooks, crop_out=ALREADY_SAME_SIZE) for d in ifg_data_paths] mlooked = [Ifg(i) for i in looks_paths] for m in mlooked: m.open() @@ -467,17 +574,17 @@ def test_same_size_multilooking(self): for ifg in mlooked: self.assertAlmostEqual(ifg.x_step, xlooks * self.xs) self.assertAlmostEqual(ifg.x_step, ylooks * self.xs) - # os.remove(ifg.data_path) + os.remove(ifg.data_path) def test_mlooked_path(): - path = 'geo_060619-061002.tif' + path = 'geo_060619-061002_unw.tif' assert mlooked_path(path, looks=2, crop_out=4) == \ - 'geo_060619-061002_2rlks_4cr.tif' + 'geo_060619-061002_unw_2rlks_4cr.tif' - path = 'some/dir/geo_060619-061002.tif' + path = 'some/dir/geo_060619-061002_unw.tif' assert mlooked_path(path, looks=4, crop_out=2) == \ - 'some/dir/geo_060619-061002_4rlks_2cr.tif' + 'some/dir/geo_060619-061002_unw_4rlks_2cr.tif' path = 'some/dir/geo_060619-061002_4rlks.tif' assert mlooked_path(path, looks=4, crop_out=8) == \ @@ -563,7 +670,9 @@ def setUp(self): from tests.common import small_data_setup self.ifgs = small_data_setup() self.ifg_paths = [i.data_path for i in self.ifgs] - prepare_ifgs(self.ifg_paths, crop_opt=1, xlooks=1, ylooks=1) + params = Configuration(common.TEST_CONF_ROIPAC).__dict__ + self.headers = [roipac.roipac_header(i.data_path, params) for i in self.ifgs] + prepare_ifgs(self.ifg_paths, crop_opt=1, xlooks=1, ylooks=1, headers=self.headers) looks_paths = [mlooked_path(d, looks=1, crop_out=1) for d in self.ifg_paths] self.ifgs_with_nan = [Ifg(i) for i in looks_paths] @@ -697,7 +806,7 @@ def common_check(self, ele, inc): sys.argv = ['dummy', self.conf_file] prepifg.main(params) # test 17 geotiffs created - geotifs = glob.glob(os.path.join(params[cf.OUT_DIR], '*_unw.tif')) + geotifs = glob.glob(os.path.join(params[cf.OUT_DIR], '*_unw_ifg.tif')) self.assertEqual(17, len(geotifs)) # test dem geotiff created demtif = glob.glob(os.path.join(params[cf.OUT_DIR], '*_dem.tif')) @@ -713,7 +822,3 @@ def common_check(self, ele, inc): self.assertEqual(18, len(mlooked_tifs)) inc = glob.glob(os.path.join(self.base_dir, '*utm_{inc}.tif'.format(inc=inc))) self.assertEqual(0, len(inc)) - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/test_prepifg_system_vs_python.py b/tests/test_prepifg_system_vs_python.py index 9ae342c6f..c6772cdc0 100644 --- a/tests/test_prepifg_system_vs_python.py +++ b/tests/test_prepifg_system_vs_python.py @@ -1,17 +1,17 @@ -import os import shutil import pytest -import numpy as np from pathlib import Path -from subprocess import check_call, check_output, CalledProcessError +from subprocess import check_call import numpy as np from pyrate.core import config as cf -from tests.common import assert_two_dirs_equal, manipulate_test_conf -TRAVIS = True if 'TRAVIS' in os.environ else False -PYTHON3P6 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.6') else False -PYTHON3P7 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.7') else False -PYTHON3P8 = True if ('TRAVIS_PYTHON_VERSION' in os.environ and os.environ['TRAVIS_PYTHON_VERSION'] == '3.8') else False -GDAL_VERSION = check_output(["gdal-config", "--version"]).decode(encoding="utf-8").split('\n')[0] +from tests.common import ( + assert_two_dirs_equal, + manipulate_test_conf, + TRAVIS, + PYTHON3P6, + PYTHON3P7, + GDAL_VERSION +) # python3.7 and gdal3.0.4 REGRESSION = PYTHON3P7 and (GDAL_VERSION == '3.0.4') # python3.7 and gdal3.0.2 @@ -23,11 +23,6 @@ def local_crop(request): return request.param -@pytest.fixture(params=[0, 1]) -def coh_mask(request): - return request.param - - @pytest.fixture() def modified_config_short(tempdir, local_crop, get_lks, coh_mask): orbfit_lks = 1 @@ -132,16 +127,19 @@ def test_prepifg_largetifs_vs_python(modified_config_largetifs, gamma_conf, crea check_call(f"pyrate prepifg -f {sr_conf}", shell=True) # convert2tif tests, 17 interferograms - assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_unw.tif", 17) + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_unw_ifg.tif", 17) # if coherence masking, compare coh files were converted if params[cf.COH_MASK]: - assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_utm.tif", 17) + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_coh.tif", 17) print("coherence files compared") + # 17 ifgs + 1 dem + 17 mlooked file + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 35) + else: + # prepifg + # 17 ifgs + 1 dem + assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 18) - # prepifg - # 17 ifgs + 1 dem - assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 18) print("==========================xxx===========================") diff --git a/tests/test_pyrate.py b/tests/test_pyrate.py index bcc5d733f..1ff65f1d3 100644 --- a/tests/test_pyrate.py +++ b/tests/test_pyrate.py @@ -131,13 +131,15 @@ def setUpClass(cls): # Turn off validation because we're in a different working dir # and relative paths in config won't be work. - params = config.get_config_params(common.TEST_CONF_ROIPAC, - validate=False) + params = config.get_config_params(common.TEST_CONF_ROIPAC) params[cf.OUT_DIR] = cls.BASE_OUT_DIR params[cf.PROCESSOR] = 0 # roipac params[cf.APS_CORRECTION] = 0 paths = glob.glob(join(cls.BASE_OUT_DIR, 'geo_*-*.tif')) params[cf.PARALLEL] = False + params[cf.INTERFEROGRAM_FILES] = [MultiplePaths(cls.BASE_OUT_DIR, p) for p in paths] + for p in params[cf.INTERFEROGRAM_FILES]: # cheat + p.sampled_path = p.converted_path process.process_ifgs(sorted(paths), params, 2, 2) if not hasattr(cls, 'ifgs'): @@ -214,18 +216,21 @@ def setUpClass(cls): # and multilooked by run_prepifg base_unw_paths = list(cf.parse_namelist(params[cf.IFG_FILE_LIST])) - base_unw_paths_ = [MultiplePaths(params[cf.OUT_DIR], b, ifglksx=params[cf.IFG_LKSX], - ifgcropopt=params[cf.IFG_CROP_OPT]) for b in base_unw_paths] + multi_paths = [MultiplePaths(params[cf.OUT_DIR], b, ifglksx=params[cf.IFG_LKSX], + ifgcropopt=params[cf.IFG_CROP_OPT]) for b in base_unw_paths] # dest_paths are tifs that have been geotif converted and multilooked - cls.dest_paths = [b.sampled_path for b in base_unw_paths_] - gtif_paths_ = conv2tif.do_geotiff(base_unw_paths_, params) - gtif_paths = [gt for gt, b in gtif_paths_] - prepifg.do_prepifg(gtif_paths, params) - tiles = pyrate.core.shared.get_tiles(cls.dest_paths[0], rows, cols) + cls.converted_paths = [b.converted_path for b in multi_paths] + cls.sampled_paths = [b.sampled_path for b in multi_paths] + from copy import copy + orig_params = copy(params) + conv2tif.main(params) + prepifg.main(orig_params) + tiles = pyrate.core.shared.get_tiles(cls.sampled_paths[0], rows, cols) ifgs = common.small_data_setup() + params[cf.INTERFEROGRAM_FILES] = multi_paths - cls.refpixel_p, cls.maxvar_p, cls.vcmt_p = process.process_ifgs(cls.dest_paths, params, rows, cols) + cls.refpixel_p, cls.maxvar_p, cls.vcmt_p = process.process_ifgs(cls.sampled_paths, params, rows, cols) cls.mst_p = common.reconstruct_mst(ifgs[0].shape, tiles, params[cf.TMPDIR]) cls.rate_p, cls.error_p, cls.samples_p = \ [common.reconstruct_stack_rate(ifgs[0].shape, tiles, params[cf.TMPDIR], t) for t in rate_types] @@ -238,20 +243,20 @@ def setUpClass(cls): params[cf.PROCESSES] = 1 params[cf.OUT_DIR] = cls.tif_dir_s params[cf.TMPDIR] = os.path.join(params[cf.OUT_DIR], cf.TMPDIR) - - base_unw_paths_ = [MultiplePaths(params[cf.OUT_DIR], b, ifglksx=params[cf.IFG_LKSX], ifgcropopt=params[ - cf.IFG_CROP_OPT]) for b in base_unw_paths] - - cls.dest_paths_s = [b.sampled_path for b in base_unw_paths_] - gtif_paths_ = conv2tif.do_geotiff(base_unw_paths_, params) - gtif_paths = [gt for gt, b in gtif_paths_] - prepifg.do_prepifg(gtif_paths, params) - cls.refpixel, cls.maxvar, cls.vcmt = process.process_ifgs(cls.dest_paths_s, params, rows, cols) + multi_paths = [MultiplePaths(params[cf.OUT_DIR], b, ifglksx=params[cf.IFG_LKSX], + ifgcropopt=params[cf.IFG_CROP_OPT]) for b in base_unw_paths] + + cls.converted_paths_s = [b.converted_path for b in multi_paths] + cls.sampled_paths_s = [b.sampled_path for b in multi_paths] + orig_params = copy(params) + conv2tif.main(params) + prepifg.main(orig_params) + params[cf.INTERFEROGRAM_FILES] = multi_paths + cls.refpixel, cls.maxvar, cls.vcmt = process.process_ifgs(cls.sampled_paths_s, params, rows, cols) cls.mst = common.reconstruct_mst(ifgs[0].shape, tiles, params[cf.TMPDIR]) cls.rate, cls.error, cls.samples = \ [common.reconstruct_stack_rate(ifgs[0].shape, tiles, params[cf.TMPDIR], t) for t in rate_types] - @classmethod def tearDownClass(cls): shutil.rmtree(cls.tif_dir, ignore_errors=True) @@ -262,7 +267,7 @@ def test_orbital_correction(self): key = 'ORBITAL_ERROR' value = 'REMOVED' - for i in common.small_data_setup(datafiles=self.dest_paths): + for i in common.small_data_setup(datafiles=self.sampled_paths): self.key_check(i, key, value) def key_check(self, ifg, key, value): @@ -276,7 +281,7 @@ def test_phase_conversion(self): key = 'DATA_UNITS' value = 'MILLIMETRES' - for i in common.small_data_setup(datafiles=self.dest_paths): + for i in common.small_data_setup(datafiles=self.sampled_paths): self.key_check(i, key, value) def test_mst_equal(self): diff --git a/tests/test_ref_phs_est.py b/tests/test_ref_phs_est.py index d216825b7..709c94371 100644 --- a/tests/test_ref_phs_est.py +++ b/tests/test_ref_phs_est.py @@ -19,6 +19,7 @@ """ import glob import os +from pathlib import Path import shutil import sys import tempfile @@ -29,6 +30,7 @@ from pyrate.core import ifgconstants as ifc, config as cf from pyrate.core.ref_phs_est import ReferencePhaseError from pyrate.core.shared import CorrectionStatusError +from pyrate.core import roipac from pyrate import prepifg, process, conv2tif from pyrate.configuration import Configuration from tests import common @@ -89,7 +91,6 @@ def setUp(self): for ifg in self.ifgs: ifg.close() - def tearDown(self): try: shutil.rmtree(self.tmp_dir) @@ -129,7 +130,10 @@ class RefPhsEstimationLegacyTestMethod1Serial(unittest.TestCase): @classmethod def setUpClass(cls): - + from pyrate.core import roipac + # start with a clean output dir + params = Configuration(common.TEST_CONF_ROIPAC).__dict__ + shutil.rmtree(params[cf.OUT_DIR]) params = Configuration(common.TEST_CONF_ROIPAC).__dict__ cls.temp_out_dir = tempfile.mkdtemp() sys.argv = ['prepifg.py', common.TEST_CONF_ROIPAC] @@ -141,13 +145,10 @@ def setUpClass(cls): params[cf.REF_EST_METHOD] = 1 params[cf.PARALLEL] = False - xlks, ylks, crop = cf.transform_params(params) - - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], - params[cf.OBS_DIR]) - - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, - params, xlks) + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + headers = [roipac.roipac_header(i, params) for i in base_ifg_paths] + dest_paths = [Path(cls.temp_out_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in params[cf.INTERFEROGRAM_FILES][:-2]] # start run_pyrate copy ifgs = common.pre_prepare_ifgs(dest_paths, params) @@ -156,7 +157,7 @@ def setUpClass(cls): refx, refy = process._ref_pixel_calc(dest_paths, params) # Estimate and remove orbit errors - pyrate.core.orbital.remove_orbital_error(ifgs, params) + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) for i in ifgs: i.close() @@ -206,7 +207,7 @@ def test_ifgs_after_ref_phs_est(self): LEGACY_REF_PHASE_DIR, f), delimiter=',') for k, j in enumerate(self.ifgs): if f.split('_corrected')[-1].split('.')[0] == \ - os.path.split(j.data_path)[-1].split('_unw_1rlks')[0]: + os.path.split(j.data_path)[-1].split('_unw_ifg_1rlks')[0]: count += 1 # all numbers equal np.testing.assert_array_almost_equal(ifg_data, @@ -228,7 +229,8 @@ class RefPhsEstimationLegacyTestMethod1Parallel(unittest.TestCase): """ @classmethod def setUpClass(cls): - + params = Configuration(common.TEST_CONF_ROIPAC).__dict__ + shutil.rmtree(params[cf.OUT_DIR]) params = Configuration(common.TEST_CONF_ROIPAC).__dict__ cls.temp_out_dir = tempfile.mkdtemp() sys.argv = ['prepifg.py', common.TEST_CONF_ROIPAC] @@ -242,11 +244,10 @@ def setUpClass(cls): xlks, ylks, crop = cf.transform_params(params) - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], - params[cf.OBS_DIR]) - - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, - params, xlks) + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + headers = [roipac.roipac_header(i, params) for i in base_ifg_paths] + dest_paths = [Path(cls.temp_out_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in params[cf.INTERFEROGRAM_FILES][:-2]] # start run_pyrate copy ifgs = common.pre_prepare_ifgs(dest_paths, params) @@ -255,7 +256,7 @@ def setUpClass(cls): refx, refy = process._ref_pixel_calc(dest_paths, params) # Estimate and remove orbit errors - pyrate.core.orbital.remove_orbital_error(ifgs, params) + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) for i in ifgs: i.close() @@ -299,7 +300,7 @@ def test_ifgs_after_ref_phs_est(self): LEGACY_REF_PHASE_DIR, f), delimiter=',') for k, j in enumerate(self.ifgs): if f.split('_corrected')[-1].split('.')[0] == \ - os.path.split(j.data_path)[-1].split('_unw_1rlks')[0]: + os.path.split(j.data_path)[-1].split('_unw_ifg_1rlks')[0]: count += 1 # all numbers equal np.testing.assert_array_almost_equal( @@ -325,7 +326,8 @@ class RefPhsEstimationLegacyTestMethod2Serial(unittest.TestCase): @classmethod def setUpClass(cls): - + params = Configuration(common.TEST_CONF_ROIPAC).__dict__ + shutil.rmtree(params[cf.OUT_DIR]) params = Configuration(common.TEST_CONF_ROIPAC).__dict__ cls.temp_out_dir = tempfile.mkdtemp() sys.argv = ['prepifg.py', common.TEST_CONF_ROIPAC] @@ -337,14 +339,10 @@ def setUpClass(cls): params[cf.REF_EST_METHOD] = 2 params[cf.PARALLEL] = False - xlks, ylks, crop = cf.transform_params(params) - - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], - params[cf.OBS_DIR]) - - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, - params, xlks) - + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + headers = [roipac.roipac_header(i, params) for i in base_ifg_paths] + dest_paths = [Path(cls.temp_out_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in params[cf.INTERFEROGRAM_FILES][:-2]] # start run_pyrate copy ifgs = common.pre_prepare_ifgs(dest_paths, params) mst_grid = common.mst_calculation(dest_paths, params) @@ -352,7 +350,7 @@ def setUpClass(cls): refx, refy = process._ref_pixel_calc(dest_paths, params) # Estimate and remove orbit errors - pyrate.core.orbital.remove_orbital_error(ifgs, params) + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) for i in ifgs: i.close() @@ -390,7 +388,7 @@ def test_ifgs_after_ref_phs_est(self): LEGACY_REF_PHASE_DIR, f), delimiter=',') for k, j in enumerate(self.ifgs): if f.split('_corrected_method2')[-1].split('.')[0] == \ - os.path.split(j.data_path)[-1].split('_unw_1rlks')[0]: + os.path.split(j.data_path)[-1].split('_unw_ifg_1rlks')[0]: count += 1 # all numbers equal np.testing.assert_array_almost_equal(ifg_data, @@ -420,10 +418,10 @@ class RefPhsEstimationLegacyTestMethod2Parallel(unittest.TestCase): # TODO: Improve the parallel tests to remove duplication from serial tests @classmethod def setUpClass(cls): - + params = Configuration(common.TEST_CONF_ROIPAC).__dict__ + shutil.rmtree(params[cf.OUT_DIR]) params = Configuration(common.TEST_CONF_ROIPAC).__dict__ cls.temp_out_dir = tempfile.mkdtemp() - sys.argv = ['prepifg.py', common.TEST_CONF_ROIPAC] params[cf.OUT_DIR] = cls.temp_out_dir params[cf.TMPDIR] = cls.temp_out_dir conv2tif.main(params) @@ -431,15 +429,14 @@ def setUpClass(cls): params[cf.OUT_DIR] = cls.temp_out_dir params[cf.REF_EST_METHOD] = 2 - params[cf.PARALLEL] = True - - xlks, ylks, crop = cf.transform_params(params) + params[cf.PARALLEL] = 1 - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], - params[cf.OBS_DIR]) + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + headers = [roipac.roipac_header(i, params) for i in base_ifg_paths] - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, - params, xlks) + # leave 2 out due to conv2tif and prepifg dems + dest_paths = [Path(cls.temp_out_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in params[cf.INTERFEROGRAM_FILES][:-2]] # start run_pyrate copy ifgs = common.pre_prepare_ifgs(dest_paths, params) @@ -447,7 +444,7 @@ def setUpClass(cls): refx, refy = process._ref_pixel_calc(dest_paths, params) # Estimate and remove orbit errors - pyrate.core.orbital.remove_orbital_error(ifgs, params) + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) for i in ifgs: i.close() @@ -484,7 +481,7 @@ def test_ifgs_after_ref_phs_est(self): LEGACY_REF_PHASE_DIR, f), delimiter=',') for k, j in enumerate(self.ifgs): if f.split('_corrected_method2')[-1].split('.')[0] == \ - os.path.split(j.data_path)[-1].split('_unw_1rlks')[0]: + os.path.split(j.data_path)[-1].split('_unw_ifg_1rlks')[0]: count += 1 # all numbers equal np.testing.assert_array_almost_equal( diff --git a/tests/test_refpixel.py b/tests/test_refpixel.py index 19aa0cf63..17839b420 100644 --- a/tests/test_refpixel.py +++ b/tests/test_refpixel.py @@ -18,15 +18,21 @@ """ import copy import unittest -import tempfile import shutil +from subprocess import check_call, run +from pathlib import Path +import pytest from numpy import nan, mean, std, isnan from pyrate.core import config as cf -from pyrate.core.refpixel import ref_pixel, _step +from pyrate.core.refpixel import ref_pixel, _step, RefPixelError +from pyrate.core import shared, ifgconstants as ifc + from pyrate import process +from pyrate.configuration import Configuration from tests.common import TEST_CONF_ROIPAC -from tests.common import small_data_setup, MockIfg, small_ifg_file_list +from tests.common import small_data_setup, MockIfg, copy_small_ifg_file_list, \ + copy_and_setup_small_data, manipulate_test_conf, assert_two_dirs_equal, PYTHON3P6 # TODO: figure out how editing resource.setrlimit fixes the error @@ -62,7 +68,7 @@ def test_missing_chipsize(self): def test_chipsize_valid(self): for illegal in [0, -1, -15, 1, 2, self.ifgs[0].ncols+1, 4, 6, 10, 20]: self.params[cf.REF_CHIP_SIZE] = illegal - self.assertRaises(ValueError, ref_pixel, self.ifgs, self.params) + self.assertRaises(RefPixelError, ref_pixel, self.ifgs, self.params) def test_minimum_fraction_missing(self): self.params[cf.REF_MIN_FRAC] = None @@ -71,18 +77,18 @@ def test_minimum_fraction_missing(self): def test_minimum_fraction_threshold(self): for illegal in [-0.1, 1.1, 1.000001, -0.0000001]: self.params[cf.REF_MIN_FRAC] = illegal - self.assertRaises(ValueError, ref_pixel, self.ifgs, self.params) + self.assertRaises(RefPixelError, ref_pixel, self.ifgs, self.params) def test_search_windows(self): # 45 is max # cells a width 3 sliding window can iterate over for illegal in [-5, -1, 0, 46, 50, 100]: self.params[cf.REFNX] = illegal - self.assertRaises(ValueError, ref_pixel, self.ifgs, self.params) + self.assertRaises(RefPixelError, ref_pixel, self.ifgs, self.params) # 40 is max # cells a width 3 sliding window can iterate over for illegal in [-5, -1, 0, 71, 85, 100]: self.params[cf.REFNY] = illegal - self.assertRaises(ValueError, ref_pixel, self.ifgs, self.params) + self.assertRaises(RefPixelError, ref_pixel, self.ifgs, self.params) def test_missing_search_windows(self): self.params[cf.REFNX] = None @@ -100,8 +106,8 @@ class ReferencePixelTests(unittest.TestCase): """ def setUp(self): - self.ifgs = small_data_setup() self.params = cf.get_config_params(TEST_CONF_ROIPAC) + self.params[cf.OUT_DIR], self.ifgs = copy_and_setup_small_data() self.params[cf.REFNX] = REFNX self.params[cf.REFNY] = REFNY self.params[cf.REF_CHIP_SIZE] = CHIPSIZE @@ -223,10 +229,12 @@ def _expected_ref_pixel(ifgs, cs): class LegacyEqualityTest(unittest.TestCase): def setUp(self): - self.ifg_paths = small_ifg_file_list() self.params = cf.get_config_params(TEST_CONF_ROIPAC) - self.params[cf.PARALLEL] = False - self.params[cf.OUT_DIR] = tempfile.mkdtemp() + self.params[cf.PARALLEL] = 0 + self.params[cf.OUT_DIR], self.ifg_paths = copy_small_ifg_file_list() + conf_file = Path(self.params[cf.OUT_DIR], 'conf_file.conf') + cf.write_config_file(params=self.params, output_conf_file=conf_file) + self.params = Configuration(conf_file).__dict__ self.params_alt_ref_frac = copy.copy(self.params) self.params_alt_ref_frac[cf.REF_MIN_FRAC] = 0.5 self.params_all_2s = copy.copy(self.params) @@ -242,6 +250,13 @@ def setUp(self): def tearDown(self): shutil.rmtree(self.params[cf.OUT_DIR]) + def test_small_test_data_ref_pixel_lat_lon_provided(self): + self.params[cf.REFX], self.params[cf.REFY] = 150.941666654, -34.218333314 + refx, refy = process._ref_pixel_calc(self.ifg_paths, self.params) + self.assertEqual(refx, 38) + self.assertEqual(refy, 58) + self.assertAlmostEqual(0.8, self.params[cf.REF_MIN_FRAC]) + def test_small_test_data_ref_pixel(self): refx, refy = process._ref_pixel_calc(self.ifg_paths, self.params) self.assertEqual(refx, 38) @@ -255,11 +270,21 @@ def test_small_test_data_ref_chipsize_15(self): self.assertEqual(refy, 7) self.assertAlmostEqual(0.5, self.params_alt_ref_frac[cf.REF_MIN_FRAC]) - def test_small_test_data_ref_all_1(self): - - refx, refy = process._ref_pixel_calc(self.ifg_paths, - self.params_all_1s) + def test_metadata(self): + refx, refy = process._ref_pixel_calc(self.ifg_paths, self.params_chipsize_15) + for i in self.ifg_paths: + ifg = shared.Ifg(i) + ifg.open(readonly=True) + md = ifg.meta_data + for k, v in zip([ifc.PYRATE_REFPIX_X, ifc.PYRATE_REFPIX_Y, ifc.PYRATE_REFPIX_LAT, + ifc.PYRATE_REFPIX_LON, ifc.PYRATE_MEAN_REF_AREA, ifc.PYRATE_STDDEV_REF_AREA], + [str(refx), str(refy), 0, 0, 0, 0]): + assert k in md # metadata present + # assert values + ifg.close() + def test_small_test_data_ref_all_1(self): + refx, refy = process._ref_pixel_calc(self.ifg_paths, self.params_all_1s) self.assertAlmostEqual(0.7, self.params_all_1s[cf.REF_MIN_FRAC]) self.assertEqual(1, self.params_all_1s[cf.REFNX]) self.assertEqual(1, self.params_all_1s[cf.REFNY]) @@ -270,11 +295,9 @@ def test_small_test_data_ref_all_1(self): class LegacyEqualityTestMultiprocessParallel(unittest.TestCase): def setUp(self): - self.ifg_paths = small_ifg_file_list() self.params = cf.get_config_params(TEST_CONF_ROIPAC) self.params[cf.PARALLEL] = True - self.params[cf.OUT_DIR] = tempfile.mkdtemp() - + self.params[cf.OUT_DIR], self.ifg_paths = copy_small_ifg_file_list() self.params_alt_ref_frac = copy.copy(self.params) self.params_alt_ref_frac[cf.REF_MIN_FRAC] = 0.5 self.params_all_2s = copy.copy(self.params) @@ -329,5 +352,34 @@ def test_small_test_data_ref_all_1(self): self.assertEqual(refy, 2) -if __name__ == "__main__": - unittest.main() +@pytest.mark.slow +@pytest.mark.skip(PYTHON3P6, reason='Skipped in python3p6') +def test_error_msg_refpixel_out_out_bounds(tempdir, gamma_conf): + "check correct latitude/longitude refpixel error is raised when specified refpixel is out of bounds" + for x, (refx, refy) in zip(['longitude', 'latitude', 'longitude and latitude'], + [(150., -34.218333314), (150.941666654, -34.), (150, -34)]): + _, out = _get_mlooked_files(gamma_conf, Path(tempdir()), refx=refx, refy=refy) + msg = "Supplied {} value is outside the bounds of the interferogram data" + assert msg.format(x) in out.stderr + + +@pytest.mark.slow +@pytest.mark.skip(PYTHON3P6, reason='Skipped in python3p6') +def test_gamma_ref_pixel_search_vs_lat_lon(tempdir, gamma_conf): + params_1, _ = _get_mlooked_files(gamma_conf, Path(tempdir()), refx=-1, refy=-1) + params_2, _ = _get_mlooked_files(gamma_conf, Path(tempdir()), refx=150.941666654, refy=-34.218333314) + assert_two_dirs_equal(params_1[cf.OUT_DIR], params_2[cf.OUT_DIR], f"*{params_1[cf.IFG_CROP_OPT]}cr.tif", 18) + + +def _get_mlooked_files(gamma_conf, tdir, refx, refy): + params = manipulate_test_conf(gamma_conf, tdir) + params[cf.REFX] = refx + params[cf.REFY] = refy + output_conf_file = 'config.conf' + output_conf = tdir.joinpath(output_conf_file) + cf.write_config_file(params=params, output_conf_file=output_conf) + check_call(f"pyrate conv2tif -f {output_conf}", shell=True) + check_call(f"pyrate prepifg -f {output_conf}", shell=True) + stdout = run(f"pyrate process -f {output_conf}", shell=True, capture_output=True, text=True) + print("============================================", stdout) + return params, stdout diff --git a/tests/test_roipac.py b/tests/test_roipac.py index 2dd432582..020f77b4f 100644 --- a/tests/test_roipac.py +++ b/tests/test_roipac.py @@ -130,7 +130,7 @@ def test_to_geotiff_ifg(self): band = ds.GetRasterBand(1) self.assertEqual(ds.RasterCount, 1) - exp_path = join(PREP_TEST_TIF, 'geo_060619-061002.tif') + exp_path = join(PREP_TEST_TIF, 'geo_060619-061002_unw.tif') exp_ds = gdal.Open(exp_path) exp_band = exp_ds.GetRasterBand(1) @@ -152,7 +152,7 @@ def test_to_geotiff_ifg(self): def test_to_geotiff_wrong_input_data(self): # ensure failure if TIF/other file used instead of binary UNW data self.dest = os.path.join(TEMPDIR, 'tmp_roipac_ifg.tif') - data_path = join(PREP_TEST_TIF, 'geo_060619-061002.tif') + data_path = join(PREP_TEST_TIF, 'geo_060619-061002_unw.tif') self.assertRaises( GeotiffException, write_fullres_geotiff, diff --git a/tests/test_shared.py b/tests/test_shared.py index 388ff58ef..c58b82595 100644 --- a/tests/test_shared.py +++ b/tests/test_shared.py @@ -22,6 +22,7 @@ import sys import tempfile import unittest +from pathlib import Path from itertools import product from numpy import isnan, where, nan from os.path import join, basename, exists @@ -148,6 +149,7 @@ class IfgIOTests(unittest.TestCase): def setUp(self): self.ifg = Ifg(join(SML_TEST_TIF, 'geo_070709-070813_unw.tif')) + self.header = join(common.SML_TEST_OBS, 'geo_070709-070813_unw.rsc') def test_open(self): self.assertTrue(self.ifg.dataset is None) @@ -166,11 +168,13 @@ def test_open_ifg_from_dataset(self): gdal.Dataset object as Dataset has already been read in """ paths = [self.ifg.data_path] + headers = [self.header] mlooked_phase_data = prepifg_helper.prepare_ifgs(paths, crop_opt=prepifg_helper.ALREADY_SAME_SIZE, xlooks=2, ylooks=2, - write_to_disc=False) + write_to_disc=False, + headers=headers) mlooked = [Ifg(m[1]) for m in mlooked_phase_data] self.assertRaises(RasterException, mlooked[0].open) @@ -198,7 +202,6 @@ def test_write(self): i.close() os.remove(dest) - def test_write_fails_on_readonly(self): # check readonly status is same before # and after open() for readonly file @@ -339,10 +342,8 @@ def setUpClass(cls): cls.params = Configuration(cls.test_conf).__dict__ cls.params[cf.OBS_DIR] = common.SML_TEST_GAMMA cls.params[cf.PROCESSOR] = 1 # gamma - file_list = list(cf.parse_namelist(os.path.join(common.SML_TEST_GAMMA, - 'ifms_17'))) - fd, cls.params[cf.IFG_FILE_LIST] = tempfile.mkstemp(suffix='.conf', - dir=cls.tif_dir) + file_list = list(cf.parse_namelist(os.path.join(common.SML_TEST_GAMMA, 'ifms_17'))) + fd, cls.params[cf.IFG_FILE_LIST] = tempfile.mkstemp(suffix='.conf', dir=cls.tif_dir) os.close(fd) # write a short filelist with only 3 gamma unws with open(cls.params[cf.IFG_FILE_LIST], 'w') as fp: @@ -357,15 +358,13 @@ def setUpClass(cls): cls.params[cf.IFG_FILE_LIST], cls.params[cf.OBS_DIR]) cls.base_unw_paths.append(common.SML_TEST_DEM_GAMMA) - xlks, ylks, crop = cf.transform_params(cls.params) # dest_paths are tifs that have been geotif converted and multilooked conv2tif.main(cls.params) prepifg.main(cls.params) - # run_prepifg.gamma_prepifg(cls.base_unw_paths, cls.params) - cls.base_unw_paths.pop() # removed dem as we don't want it in ifgs - cls.dest_paths = cf.get_dest_paths( - cls.base_unw_paths, crop, cls.params, xlks) + cls.dest_paths = [Path(cls.tif_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in cls.params[cf.INTERFEROGRAM_FILES][:-2]] + cls.ifgs = common.small_data_setup(datafiles=cls.dest_paths) @classmethod diff --git a/tests/test_stackrate.py b/tests/test_stackrate.py index afbeca847..08c9f5e0d 100644 --- a/tests/test_stackrate.py +++ b/tests/test_stackrate.py @@ -21,15 +21,16 @@ import shutil import tempfile import unittest +from pathlib import Path -from numpy import eye, array, ones +from numpy import eye, array, ones, nan import numpy as np -from numpy.testing import assert_array_almost_equal +from numpy.testing import assert_array_almost_equal, assert_array_equal import pyrate.core.orbital import tests.common -from pyrate.core import shared, ref_phs_est as rpe, config as cf, covariance as vcm_module -from pyrate.core.stack import stack_rate +from pyrate.core import shared, ref_phs_est as rpe, config as cf, covariance as vcm_module, roipac +from pyrate.core.stack import stack_rate_array, stack_rate_pixel, mask_rate from pyrate import process, prepifg, conv2tif from pyrate.configuration import Configuration from tests.common import (SML_TEST_DIR, prepare_ifgs_without_phase, @@ -47,32 +48,60 @@ def __init__(self, timespan, phase): self.phase_data = array([[phase]]) -class StackRateTests(unittest.TestCase): +class StackRatePixelTests(unittest.TestCase): """ - Tests the weighted least squares algorithm for determinining + Tests the weighted least squares algorithm for determining the best fitting velocity """ def setUp(self): - phase = [0.5, 3.5, 4, 2.5, 3.5, 1] - timespan = [0.1, 0.7, 0.8, 0.5, 0.7, 0.2] - self.ifgs = [SinglePixelIfg(s, p) for s, p in zip(timespan, phase)] - - def test_stack_rate(self): + self.phase = array([0.5, 3.5, 4, 2.5, 3.5, 1]) + self.timespan = array([[0.1, 0.7, 0.8, 0.5, 0.7, 0.2]]) + self.vcmt = eye(6, 6) + self.mst = ones((6, 1, 1)) + self.mst[4] = 0 + self.params = default_params() + + def test_stack_rate_pixel(self): # Simple test with one pixel and equal weighting exprate = array([[5.0]]) experr = array([[0.836242010007091]]) expsamp = array([[5]]) - vcmt = eye(6, 6) - mst = ones((6, 1, 1)) - mst[4] = 0 - params = default_params() - rate, error, samples = stack_rate(self.ifgs, params, vcmt, mst) + rate, error, samples = stack_rate_pixel(self.phase, self.mst, self.vcmt, + self.timespan, self.params['nsig'], self.params['pthr']) assert_array_almost_equal(rate, exprate) assert_array_almost_equal(error, experr) assert_array_almost_equal(samples, expsamp) +class MaskRateTests(unittest.TestCase): + """ + Test the maxsig threshold masking algorithm + """ + + def setUp(self): + self.r = array([5.0, 4.5]) # rates for 2 pixels + self.e = array([1.1, 2.1]) # errors for 2 pixels + + def test_mask_rate_maxsig1(self): + # both rate and error values masked + rate, error = mask_rate(self.r, self.e, 1) + assert_array_equal(rate, array([nan, nan])) + assert_array_equal(error, array([nan, nan])) + + def test_mask_rate_maxsig2(self): + # one rate and one error masked + rate, error = mask_rate(self.r, self.e, 2) + assert_array_equal(rate, array([5.0, nan])) + assert_array_equal(error, array([1.1, nan])) + + def test_mask_rate_maxsig3(self): + # No values masked in rate or error + rate, error = mask_rate(self.r, self.e, 3) + assert_array_equal(rate, self.r) + assert_array_equal(error, self.e) + + class LegacyEqualityTest(unittest.TestCase): """ Tests equality with legacy data @@ -93,9 +122,10 @@ def setUpClass(cls): xlks, _, crop = cf.transform_params(params) - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], params[cf.OBS_DIR]) - - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, params, xlks) + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + headers = [roipac.roipac_header(i, params) for i in base_ifg_paths] + dest_paths = [Path(cls.temp_out_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in params[cf.INTERFEROGRAM_FILES][:-2]] # start run_pyrate copy ifgs = pre_prepare_ifgs(dest_paths, params) mst_grid = tests.common.mst_calculation(dest_paths, params) @@ -103,7 +133,7 @@ def setUpClass(cls): refx, refy = process._ref_pixel_calc(dest_paths, params) # Estimate and remove orbit errors - pyrate.core.orbital.remove_orbital_error(ifgs, params) + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) ifgs = prepare_ifgs_without_phase(dest_paths, params) for ifg in ifgs: ifg.close() @@ -146,34 +176,34 @@ def test_stack_rate_full_parallel(self): """ python multiprocessing by rows vs serial """ - np.testing.assert_array_almost_equal(self.rate, self.rate_s, decimal=3) + assert_array_almost_equal(self.rate, self.rate_s, decimal=3) def test_stackrate_error_parallel(self): """ python multiprocessing by rows vs serial """ - np.testing.assert_array_almost_equal(self.error, self.error_s, decimal=3) + assert_array_almost_equal(self.error, self.error_s, decimal=3) def test_stackrate_samples_parallel(self): """ python multiprocessing by rows vs serial """ - np.testing.assert_array_almost_equal(self.samples, self.samples_s, decimal=3) + assert_array_almost_equal(self.samples, self.samples_s, decimal=3) def test_stack_rate(self): """ Compare with legacy data """ - np.testing.assert_array_almost_equal(self.rate_s, self.rate_container, decimal=3) + assert_array_almost_equal(self.rate_s, self.rate_container, decimal=3) def test_stackrate_error(self): """ Compare with legacy data """ - np.testing.assert_array_almost_equal(self.error_s, self.error_container, decimal=3) + assert_array_almost_equal(self.error_s, self.error_container, decimal=3) def test_stackrate_samples(self): """ Compare with legacy data """ - np.testing.assert_array_almost_equal(self.samples_s, self.samples_container, decimal=3) + assert_array_almost_equal(self.samples_s, self.samples_container, decimal=3) diff --git a/tests/test_timeseries.py b/tests/test_timeseries.py index 850788432..87049a3d7 100644 --- a/tests/test_timeseries.py +++ b/tests/test_timeseries.py @@ -22,6 +22,7 @@ import sys import tempfile import unittest +from pathlib import Path from datetime import date, timedelta from numpy import nan, asarray, where import numpy as np @@ -29,7 +30,7 @@ import pyrate.core.orbital import tests.common as common -from pyrate.core import ref_phs_est as rpe, config as cf, mst, covariance +from pyrate.core import ref_phs_est as rpe, config as cf, mst, covariance, roipac from pyrate import process, prepifg, conv2tif from pyrate.configuration import Configuration from pyrate.core.timeseries import time_series @@ -120,18 +121,16 @@ def setUpClass(cls): params[cf.REF_EST_METHOD] = 2 - xlks, ylks, crop = cf.transform_params(params) - - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], - params[cf.OBS_DIR]) - - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, params, xlks) + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + headers = [roipac.roipac_header(i, params) for i in base_ifg_paths] + dest_paths = [Path(cls.temp_out_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in params[cf.INTERFEROGRAM_FILES][:-2]] # start run_pyrate copy ifgs = common.pre_prepare_ifgs(dest_paths, params) mst_grid = common.mst_calculation(dest_paths, params) refx, refy = process._ref_pixel_calc(dest_paths, params) # Estimate and remove orbit errors - pyrate.core.orbital.remove_orbital_error(ifgs, params) + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) ifgs = common.prepare_ifgs_without_phase(dest_paths, params) for ifg in ifgs: ifg.close() @@ -214,12 +213,10 @@ def setUpClass(cls): params[cf.REF_EST_METHOD] = 2 - xlks, ylks, crop = cf.transform_params(params) - - base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST], - params[cf.OBS_DIR]) - - dest_paths = cf.get_dest_paths(base_ifg_paths, crop, params, xlks) + base_ifg_paths = [c.unwrapped_path for c in params[cf.INTERFEROGRAM_FILES]] + headers = [roipac.roipac_header(i, params) for i in base_ifg_paths] + dest_paths = [Path(cls.temp_out_dir).joinpath(Path(c.sampled_path).name).as_posix() + for c in params[cf.INTERFEROGRAM_FILES][:-2]] # start run_pyrate copy ifgs = common.pre_prepare_ifgs(dest_paths, params) mst_grid = common.mst_calculation(dest_paths, params) @@ -227,7 +224,7 @@ def setUpClass(cls): refx, refy = process._ref_pixel_calc(dest_paths, params) # Estimate and remove orbit errors - pyrate.core.orbital.remove_orbital_error(ifgs, params) + pyrate.core.orbital.remove_orbital_error(ifgs, params, headers) ifgs = common.prepare_ifgs_without_phase(dest_paths, params) for ifg in ifgs: ifg.close() diff --git a/utils/colourmap.txt b/utils/colourmap.txt deleted file mode 100644 index 9073cacec..000000000 --- a/utils/colourmap.txt +++ /dev/null @@ -1,269 +0,0 @@ -nan 0 0 0 0 -elevation_value 103 0 31 255 -elevation_value 162 25 43 255 -elevation_value 201 72 67 255 -elevation_value 228 130 104 255 -elevation_value 246 182 154 255 -elevation_value 251 221 205 255 -elevation_value 242 239 238 255 -elevation_value 213 231 240 255 -elevation_value 166 207 227 255 -elevation_value 107 172 208 255 -elevation_value 58 132 187 255 -elevation_value 29 92 154 255 -elevation_value 103 0 31 255 -elevation_value 106 1 31 255 -elevation_value 109 2 32 255 -elevation_value 112 3 32 255 -elevation_value 115 4 33 255 -elevation_value 118 5 33 255 -elevation_value 120 6 34 255 -elevation_value 123 7 34 255 -elevation_value 126 8 35 255 -elevation_value 129 9 35 255 -elevation_value 132 10 36 255 -elevation_value 135 11 36 255 -elevation_value 137 12 37 255 -elevation_value 140 13 38 255 -elevation_value 143 14 38 255 -elevation_value 146 16 39 255 -elevation_value 148 17 39 255 -elevation_value 151 18 40 255 -elevation_value 153 20 41 255 -elevation_value 156 21 41 255 -elevation_value 158 23 42 255 -elevation_value 161 24 43 255 -elevation_value 163 26 44 255 -elevation_value 166 27 44 255 -elevation_value 168 29 45 255 -elevation_value 170 31 46 255 -elevation_value 172 33 47 255 -elevation_value 174 35 48 255 -elevation_value 177 37 49 255 -elevation_value 179 39 50 255 -elevation_value 180 41 51 255 -elevation_value 182 43 52 255 -elevation_value 184 45 53 255 -elevation_value 186 48 54 255 -elevation_value 188 50 55 255 -elevation_value 190 52 57 255 -elevation_value 191 55 58 255 -elevation_value 193 57 59 255 -elevation_value 194 60 60 255 -elevation_value 196 62 62 255 -elevation_value 198 65 63 255 -elevation_value 199 68 64 255 -elevation_value 200 70 66 255 -elevation_value 202 73 67 255 -elevation_value 203 76 69 255 -elevation_value 205 78 70 255 -elevation_value 206 81 72 255 -elevation_value 207 84 73 255 -elevation_value 209 87 75 255 -elevation_value 210 89 76 255 -elevation_value 211 92 78 255 -elevation_value 213 95 80 255 -elevation_value 214 98 82 255 -elevation_value 215 100 83 255 -elevation_value 217 103 85 255 -elevation_value 218 106 87 255 -elevation_value 219 109 89 255 -elevation_value 220 111 91 255 -elevation_value 222 114 92 255 -elevation_value 223 117 94 255 -elevation_value 224 119 96 255 -elevation_value 225 122 98 255 -elevation_value 226 125 100 255 -elevation_value 227 127 102 255 -elevation_value 228 130 104 255 -elevation_value 230 133 106 255 -elevation_value 231 135 108 255 -elevation_value 232 138 111 255 -elevation_value 233 141 113 255 -elevation_value 234 143 115 255 -elevation_value 235 146 117 255 -elevation_value 236 148 119 255 -elevation_value 236 151 122 255 -elevation_value 237 153 124 255 -elevation_value 238 156 126 255 -elevation_value 239 158 128 255 -elevation_value 240 161 131 255 -elevation_value 241 163 133 255 -elevation_value 241 165 136 255 -elevation_value 242 168 138 255 -elevation_value 243 170 140 255 -elevation_value 243 172 143 255 -elevation_value 244 175 145 255 -elevation_value 245 177 148 255 -elevation_value 245 179 150 255 -elevation_value 246 181 153 255 -elevation_value 246 184 155 255 -elevation_value 247 186 158 255 -elevation_value 247 188 160 255 -elevation_value 248 190 163 255 -elevation_value 248 192 165 255 -elevation_value 248 194 168 255 -elevation_value 249 196 170 255 -elevation_value 249 198 173 255 -elevation_value 249 200 175 255 -elevation_value 249 202 178 255 -elevation_value 250 204 180 255 -elevation_value 250 205 183 255 -elevation_value 250 207 185 255 -elevation_value 250 209 188 255 -elevation_value 250 211 190 255 -elevation_value 250 212 192 255 -elevation_value 250 214 195 255 -elevation_value 251 216 197 255 -elevation_value 251 217 199 255 -elevation_value 251 219 201 255 -elevation_value 251 220 204 255 -elevation_value 251 222 206 255 -elevation_value 250 223 208 255 -elevation_value 250 224 210 255 -elevation_value 250 226 212 255 -elevation_value 250 227 214 255 -elevation_value 250 228 216 255 -elevation_value 250 229 218 255 -elevation_value 249 230 219 255 -elevation_value 249 231 221 255 -elevation_value 249 232 223 255 -elevation_value 248 233 224 255 -elevation_value 248 234 226 255 -elevation_value 248 235 227 255 -elevation_value 247 236 229 255 -elevation_value 247 236 230 255 -elevation_value 246 237 232 255 -elevation_value 245 238 233 255 -elevation_value 245 238 234 255 -elevation_value 244 238 235 255 -elevation_value 243 239 236 255 -elevation_value 243 239 237 255 -elevation_value 242 239 238 255 -elevation_value 241 239 239 255 -elevation_value 240 240 239 255 -elevation_value 239 240 240 255 -elevation_value 238 240 240 255 -elevation_value 237 239 241 255 -elevation_value 236 239 241 255 -elevation_value 234 239 242 255 -elevation_value 233 239 242 255 -elevation_value 232 238 242 255 -elevation_value 231 238 242 255 -elevation_value 229 238 242 255 -elevation_value 228 237 242 255 -elevation_value 227 237 242 255 -elevation_value 225 236 242 255 -elevation_value 224 235 242 255 -elevation_value 222 235 242 255 -elevation_value 220 234 241 255 -elevation_value 219 233 241 255 -elevation_value 217 233 241 255 -elevation_value 215 232 241 255 -elevation_value 213 231 240 255 -elevation_value 212 230 240 255 -elevation_value 210 229 239 255 -elevation_value 208 228 239 255 -elevation_value 206 227 238 255 -elevation_value 204 226 238 255 -elevation_value 202 225 237 255 -elevation_value 200 224 237 255 -elevation_value 198 223 236 255 -elevation_value 196 222 236 255 -elevation_value 194 221 235 255 -elevation_value 191 220 235 255 -elevation_value 189 219 234 255 -elevation_value 187 218 233 255 -elevation_value 185 217 233 255 -elevation_value 182 215 232 255 -elevation_value 180 214 232 255 -elevation_value 178 213 231 255 -elevation_value 175 212 230 255 -elevation_value 173 210 229 255 -elevation_value 170 209 229 255 -elevation_value 168 208 228 255 -elevation_value 165 206 227 255 -elevation_value 163 205 226 255 -elevation_value 160 203 226 255 -elevation_value 157 202 225 255 -elevation_value 155 200 224 255 -elevation_value 152 199 223 255 -elevation_value 149 197 222 255 -elevation_value 147 196 222 255 -elevation_value 144 194 221 255 -elevation_value 141 193 220 255 -elevation_value 138 191 219 255 -elevation_value 135 189 218 255 -elevation_value 133 188 217 255 -elevation_value 130 186 216 255 -elevation_value 127 184 215 255 -elevation_value 124 183 214 255 -elevation_value 121 181 213 255 -elevation_value 118 179 212 255 -elevation_value 116 177 211 255 -elevation_value 113 175 210 255 -elevation_value 110 174 209 255 -elevation_value 107 172 208 255 -elevation_value 104 170 207 255 -elevation_value 102 168 206 255 -elevation_value 99 166 205 255 -elevation_value 96 164 204 255 -elevation_value 94 162 203 255 -elevation_value 91 161 202 255 -elevation_value 88 159 201 255 -elevation_value 86 157 200 255 -elevation_value 83 155 199 255 -elevation_value 81 153 198 255 -elevation_value 79 151 197 255 -elevation_value 76 149 196 255 -elevation_value 74 147 195 255 -elevation_value 72 146 194 255 -elevation_value 70 144 194 255 -elevation_value 68 142 193 255 -elevation_value 66 140 192 255 -elevation_value 64 138 191 255 -elevation_value 62 136 190 255 -elevation_value 60 135 189 255 -elevation_value 58 133 188 255 -elevation_value 57 131 187 255 -elevation_value 55 129 185 255 -elevation_value 54 127 184 255 -elevation_value 52 126 183 255 -elevation_value 50 124 182 255 -elevation_value 49 122 181 255 -elevation_value 48 120 180 255 -elevation_value 46 118 178 255 -elevation_value 45 116 177 255 -elevation_value 43 114 176 255 -elevation_value 42 113 174 255 -elevation_value 41 111 173 255 -elevation_value 39 109 171 255 -elevation_value 38 107 169 255 -elevation_value 37 105 168 255 -elevation_value 36 103 166 255 -elevation_value 34 101 164 255 -elevation_value 33 99 162 255 -elevation_value 32 97 160 255 -elevation_value 31 95 158 255 -elevation_value 30 93 156 255 -elevation_value 29 91 154 255 -elevation_value 27 89 151 255 -elevation_value 26 87 149 255 -elevation_value 25 85 147 255 -elevation_value 24 83 144 255 -elevation_value 23 81 142 255 -elevation_value 22 79 139 255 -elevation_value 20 77 137 255 -elevation_value 19 75 134 255 -elevation_value 18 73 131 255 -elevation_value 17 71 129 255 -elevation_value 16 69 126 255 -elevation_value 15 67 123 255 -elevation_value 14 65 120 255 -elevation_value 13 63 117 255 -elevation_value 12 61 114 255 -elevation_value 10 59 112 255 -elevation_value 9 56 109 255 -elevation_value 8 54 106 255 -elevation_value 7 52 103 255 -elevation_value 6 50 100 255 \ No newline at end of file diff --git a/utils/crop_ifgs.py b/utils/crop_ifgs.py index 4cb91d07f..58c782432 100644 --- a/utils/crop_ifgs.py +++ b/utils/crop_ifgs.py @@ -17,7 +17,7 @@ python utility to crop an interferogram example usage: -python pyrate/utils/crop_ifgs.py -i tests/test_data/small_test/tif/geo_060619-061002.tif +python pyrate/utils/crop_ifgs.py -i tests/test_data/small_test/tif/geo_060619-061002_unw.tif -o out.tif -e '150.91 -34.229999976 150.949166651 -34.17' """ diff --git a/utils/developer_hints.txt b/utils/developer_hints.txt index e03085fbd..86284ed83 100644 --- a/utils/developer_hints.txt +++ b/utils/developer_hints.txt @@ -7,9 +7,11 @@ Summary: 1) create new 'release' in Github. this will tag a commit hash. 2) copy link to tar.gz Source Code from github release, add to download_url in setup.py 3) update version number in setup.py (to match the Github release) -4) `python setup.py sdist` - prepare package +4) `python setup.py sdist bdist_wheel` - prepare package 5) `twine upload dist/*` - upload package to PyPI +check out here also: https://packaging.python.org/tutorials/packaging-projects/ + --------------------------------------------- # Clone Github repo cd ~