From 2358841a5c574b08a9ddade8e531c8f4714950af Mon Sep 17 00:00:00 2001 From: Mitchell Wheeler Date: Tue, 31 May 2022 16:56:20 +1000 Subject: [PATCH 1/5] A whole bunch of pylint warning/error/refactor fixes, these bring pylint messages from hundreds down to just 50 - the ones that remain are larger issues mostly, or things I'm not too sure on a nice way to fix just yet --- pyrate/configuration.py | 215 ++++++++++++++------- pyrate/constants.py | 51 +++-- pyrate/conv2tif.py | 37 ++-- pyrate/core/algorithm.py | 33 ++-- pyrate/core/aps.py | 6 +- pyrate/core/covariance.py | 20 +- pyrate/core/dem_error.py | 79 +++++--- pyrate/core/gamma.py | 74 ++++--- pyrate/core/geometry.py | 102 +++++++--- pyrate/core/logger.py | 26 +-- pyrate/core/mpiops.py | 23 ++- pyrate/core/mst.py | 26 +-- pyrate/core/orbital.py | 55 +++--- pyrate/core/phase_closure/closure_check.py | 107 ++++++---- pyrate/core/phase_closure/collect_loops.py | 21 +- pyrate/core/phase_closure/mst_closure.py | 47 +++-- pyrate/core/phase_closure/plot_closure.py | 35 ++-- pyrate/core/phase_closure/sum_closure.py | 83 +++++--- pyrate/core/prepifg_helper.py | 53 +++-- pyrate/core/ref_phs_est.py | 39 ++-- pyrate/core/refpixel.py | 125 +++++++----- pyrate/core/roipac.py | 35 ++-- pyrate/core/shared.py | 189 ++++++++++-------- pyrate/core/stack.py | 67 ++++--- pyrate/core/timeseries.py | 33 ++-- pyrate/correct.py | 45 +++-- pyrate/default_parameters.py | 4 +- pyrate/main.py | 83 +++++--- pyrate/merge.py | 113 ++++++----- pyrate/prepifg.py | 172 ++++++++++------- tests/test_geometry.py | 4 +- tests/test_refpixel.py | 8 +- tests/test_system.py | 4 +- 33 files changed, 1220 insertions(+), 794 deletions(-) diff --git a/pyrate/configuration.py b/pyrate/configuration.py index 47cde352d..3db878210 100644 --- a/pyrate/configuration.py +++ b/pyrate/configuration.py @@ -23,8 +23,8 @@ from typing import Union import pyrate.constants as C -from pyrate.constants import NO_OF_PARALLEL_PROCESSES, sixteen_digits_pattern, twelve_digits_pattern, ORB_ERROR_DIR, \ - DEM_ERROR_DIR, TEMP_MLOOKED_DIR +from pyrate.constants import NO_OF_PARALLEL_PROCESSES, sixteen_digits_pattern, \ + twelve_digits_pattern, ORB_ERROR_DIR, DEM_ERROR_DIR, TEMP_MLOOKED_DIR from pyrate.core import ifgconstants as ifg from pyrate.default_parameters import PYRATE_DEFAULT_CONFIGURATION from pyrate.core.algorithm import factorise_integer @@ -32,61 +32,81 @@ def set_parameter_value(data_type, input_value, default_value, required, input_name): - if len(input_value) < 1: + """ + Convers a user-provided value into a final value, by applying data types and + default values on the input. + + This function will raise an error if an input value is not given, but required. + """ + if input_value is not None and len(input_value) < 1: input_value = None if required: # pragma: no cover - raise ValueError("A required parameter is missing value in input configuration file: " + str(input_name)) + msg = f"A required parameter is missing from the configuration file: {input_name}" + raise ValueError(msg) if input_value is not None: if str(data_type) in "path": return Path(input_value) return data_type(input_value) + return default_value -def validate_parameter_value(input_name, input_value, min_value=None, max_value=None, possible_values=None): +def validate_parameter_value( + input_name: str, + input_value, + min_value=None, + max_value=None, + possible_values=None +): + """Validates that a supplied value sits within a range or set of allowed values.""" + if isinstance(input_value, PurePath): if not Path.exists(input_value): # pragma: no cover - raise ValueError("Given path: " + str(input_value) + " does not exist.") + raise ValueError(f"Given path: {input_value} does not exist.") + if input_value is not None: 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) + ". Provide a value greater than or equal to " + str(min_value) + ".") - if input_value is not None: + msg = f"Invalid value for {input_name} supplied: {input_value}. " \ + f"Provide a value greater than or equal to {min_value}." + raise ValueError(msg) + 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) + ". Provide a value less than or equal to " + str(max_value) + ".") + msg = f"Invalid value for {input_name} supplied: {input_value}. " \ + f"Provide a value less than or equal to {max_value}." + raise ValueError(msg) 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) + ". Provide one of these values: " + str(possible_values) + ".") + msg = f"Invalid value for {input_name} supplied: {input_value}. " \ + f"Provide one of these values: {possible_values}." + raise ValueError(msg) return True def validate_file_list_values(file_list, no_of_epochs): + """Validates that the provided IFG files have enough epochs""" if file_list is None: # pragma: no cover raise ValueError("No value supplied for input file list: " + str(file_list)) files = parse_namelist(file_list) - for f in files: - if not Path(f).exists(): # pragma: no cover - raise ConfigException(f"{f} does not exist") - else: - matches = extract_epochs_from_filename(filename_with_epochs=f) - if len(matches) < no_of_epochs: # pragma: no cover - raise ConfigException(f"the number of epochs in {f} names are less the required number: {no_of_epochs}") + for file in files: + if not Path(file).exists(): # pragma: no cover + raise ConfigException(f"{file} does not exist") + + matches = extract_epochs_from_filename(filename_with_epochs=file) + if len(matches) < no_of_epochs: # pragma: no cover + msg = f"the number of epochs in {file} names are less " \ + f"the required number: {no_of_epochs}" + raise ConfigException(msg) class MultiplePaths: + """A utility class which provides paths related to an interferogram file.""" def __init__(self, file_name: str, params: dict, input_type: InputTypes = InputTypes.IFG): - self.input_type = input_type out_dir = params[C.OUT_DIR] tempdir = params[C.TEMP_MLOOKED_DIR] @@ -98,13 +118,14 @@ def __init__(self, file_name: str, params: dict, input_type: InputTypes = InputT if d is None: # could be 6 digit epoch dates d = re.search(twelve_digits_pattern, b.stem) if d is None: - raise ValueError(f"{input_type.value} filename does not contain two 8- or 6-digit date strings") + msg = f"{input_type.value} filename does not contain two 8 or 6 digit date strings" + raise ValueError(msg) filestr = d.group() + '_' else: filestr = '' - dir_exists = input_type.value in InputTypes.dir_map.value.keys() - anchoring_dir = Path(out_dir).joinpath(InputTypes.dir_map.value[input_type.value]) \ + dir_exists = input_type.value in InputTypes.DIR_MAP.value.keys() + anchoring_dir = Path(out_dir).joinpath(InputTypes.DIR_MAP.value[input_type.value]) \ if dir_exists else Path(out_dir) if b.suffix == ".tif": @@ -113,7 +134,8 @@ def __init__(self, file_name: str, params: dict, input_type: InputTypes = InputT self.sampled_path = anchoring_dir.joinpath(filestr + input_type.value + '.tif') else: self.unwrapped_path = b.as_posix() - converted_path = anchoring_dir.joinpath(b.stem.split('.')[0] + '_' + b.suffix[1:]).with_suffix('.tif') + converted_path = b.stem.split('.')[0] + '_' + b.suffix[1:] + converted_path = (anchoring_dir / converted_path).with_suffix('.tif') self.sampled_path = converted_path.with_name(filestr + input_type.value + '.tif') # tmp_sampled_paths are used after prepifg, during correct steps @@ -123,6 +145,7 @@ def __init__(self, file_name: str, params: dict, input_type: InputTypes = InputT @staticmethod def orb_error_path(ifg_path: Union[str, Path], params) -> Path: + """Returns the orbit error path for a specific interferogram file.""" if isinstance(ifg_path, str): ifg_path = Path(ifg_path) return Path(params[C.OUT_DIR], C.ORB_ERROR_DIR, @@ -135,6 +158,7 @@ def orb_error_path(ifg_path: Union[str, Path], params) -> Path: @staticmethod def dem_error_path(ifg_path: Union[str, Path], params) -> Path: + """Returns the DEM error path for a specific interferogram file.""" if isinstance(ifg_path, str): ifg_path = Path(ifg_path) return Path(params[C.OUT_DIR], C.DEM_ERROR_DIR, @@ -142,6 +166,7 @@ def dem_error_path(ifg_path: Union[str, Path], params) -> Path: @staticmethod def aps_error_path(ifg_path: Union[str, Path], params) -> Path: + """Returns the APS error path for a specific interferogram file.""" if isinstance(ifg_path, str): ifg_path = Path(ifg_path) return Path(params[C.OUT_DIR], C.APS_ERROR_DIR, @@ -156,28 +181,62 @@ def aps_error_path(ifg_path: Union[str, Path], params) -> Path: ]) + '_aps_error.npy') def __str__(self): # pragma: no cover - st = "" + value = "" if self.unwrapped_path is not None: - st += """\nunwrapped_path = """ + self.unwrapped_path + value += """\nunwrapped_path = """ + self.unwrapped_path else: - st += """\nunwrapped_path = None""" - st += """ - converted_path = """ + self.converted_path + """ - sampled_path = """ + self.sampled_path + """ - tmp_sampled_path = """ + self.tmp_sampled_path + """ + value += """\nunwrapped_path = None""" + value += f""" + converted_path = {self.converted_path} + sampled_path = {self.sampled_path} + tmp_sampled_path = {self.tmp_sampled_path} """ - return st + return value class Configuration: + """ + The main configuration class for PyRate, which allows access to the values of + configurable properties. - def __init__(self, config_file_path): + :param config_file_path: The path to the configuration text file to load + :ivar outdir: The PyRate output directory + """ + outdir: str + refchipsize: int + + # pylint: disable=invalid-name + # JUSTIFICATION: these are long-established in the code already, out of scope to fix just yet. + + # TODO: Write documentation for each of the configuration options, at the very least refering + # to some other place it's all documented (currently I'm not sure they are?) + # + # pylint: disable=missing-function-docstring + # JUSTIFICATION: This is a whole job in and of itself, the configuration options need docs... + # and in writing those docs, the typing info about them can be established / QA can be improved + + # pylint: disable=too-many-locals,too-many-branches,too-many-statements + # JUSTIFICATION: cleaning up configuration class is out of scope currently, will probably be + # tackled when things are restructured (shared.py and configuration.py are cluttered). + + def __init__(self, config_file_path: Union[str, Path]): + # Promote config to path object + if not isinstance(config_file_path, Path): + config_file_path = Path(config_file_path) + + # Setup default values (in case they're not in the config file) + self.parallel = False + self.processes = 1 + self.cohfilelist = None + self.basefilelist = None + self.demfile = None + + # Load config file parser = ConfigParser() parser.optionxform = str # mimic header to fulfil the requirement for configparser - with open(config_file_path) as stream: - parser.read_string("[root]\n" + stream.read()) + parser.read_string("[root]\n" + config_file_path.read_text("utf-8")) for key, value in parser["root"].items(): self.__dict__[key] = value @@ -187,7 +246,8 @@ def __init__(self, config_file_path): # custom correct sequence if 'correct' section is provided in config if 'correct' in parser and 'steps' in parser['correct']: - self.__dict__['correct'] = list(filter(None, parser['correct'].get('steps').splitlines())) + filtered_correct = filter(None, parser['correct'].get('steps').splitlines()) + self.__dict__['correct'] = list(filtered_correct) else: self.__dict__['correct'] = [ 'orbfit', @@ -226,8 +286,8 @@ def __init__(self, config_file_path): # bespoke parameter validation if self.refchipsize % 2 != 1: # pragma: no cover if self.refchipsize - 1 > 1: - # Configuration parameters refchipsize must be odd - # values too large (>101) will slow down the process without significant gains in results. + # Configuration parameters refchipsize must be odd, values too large (>101) + # will slow down the process without significant gains in results. self.refchipsize = self.refchipsize - 1 # calculate rows and cols if not supplied @@ -235,10 +295,10 @@ def __init__(self, config_file_path): self.rows, self.cols = int(self.rows), int(self.cols) else: if NO_OF_PARALLEL_PROCESSES > 1: # i.e. mpirun - self.rows, self.cols = [int(num) for num in factorise_integer(NO_OF_PARALLEL_PROCESSES)] + self.rows, self.cols = factorise_integer(NO_OF_PARALLEL_PROCESSES) else: if self.parallel: # i.e. joblib parallelism - self.rows, self.cols = [int(num) for num in factorise_integer(self.processes)] + self.rows, self.cols = factorise_integer(self.processes) else: # i.e. serial self.rows, self.cols = 1, 1 @@ -305,20 +365,25 @@ 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', input_type=InputTypes.COH) + paths = self.__get_files_from_attr('cohfilelist', input_type=InputTypes.COH) + self.coherence_file_paths = paths if self.basefilelist is not None: # if self.processor != 0: # not roipac validate_file_list_values(self.basefilelist, 1) - self.baseline_file_paths = self.__get_files_from_attr('basefilelist', input_type=InputTypes.BASE) + paths = self.__get_files_from_attr('basefilelist', input_type=InputTypes.BASE) + self.baseline_file_paths = paths - self.header_file_paths = self.__get_files_from_attr('hdrfilelist', input_type=InputTypes.HEADER) + paths = self.__get_files_from_attr('hdrfilelist', input_type=InputTypes.HEADER) + self.header_file_paths = paths - self.interferogram_files = self.__get_files_from_attr('ifgfilelist') + paths = self.__get_files_from_attr('ifgfilelist') + self.interferogram_files = paths self.dem_file = MultiplePaths(self.demfile, self.__dict__, input_type=InputTypes.DEM) # backward compatibility for string paths + # pylint: disable=consider-using-dict-items for key in self.__dict__: if isinstance(self.__dict__[key], PurePath): self.__dict__[key] = str(self.__dict__[key]) @@ -355,7 +420,9 @@ def phase_closure_filtered_ifgs_list(params): def refresh_ifg_list(self, params): # update params dict filtered_ifgs_list = self.phase_closure_filtered_ifgs_list(params) files = parse_namelist(filtered_ifgs_list.as_posix()) - params[C.INTERFEROGRAM_FILES] = [MultiplePaths(p, self.__dict__, input_type=InputTypes.IFG) for p in files] + params[C.INTERFEROGRAM_FILES] = [ + MultiplePaths(p, self.__dict__, input_type=InputTypes.IFG) for p in files + ] return params @staticmethod @@ -382,6 +449,7 @@ def closure(self): closure_d = Path(self.phase_closure_dir) class Closure: + """Just a simple data class""" def __init__(self): self.closure = closure_d.joinpath('closure.npy') self.ifgs_breach_count = closure_d.joinpath('ifgs_breach_count.npy') @@ -393,17 +461,25 @@ def __init__(self): @staticmethod def coherence_stats(params): coh_d = Path(params[C.COHERENCE_DIR]) - return {k: coh_d.joinpath(k.lower() + '.tif').as_posix() for k in [ifg.COH_MEDIAN, ifg.COH_MEAN, ifg.COH_STD]} + keys = [ifg.COH_MEDIAN, ifg.COH_MEAN, ifg.COH_STD] + return { k: coh_d.joinpath(k.lower() + '.tif').as_posix() for k in keys } @staticmethod def geometry_files(params): geom_dir = Path(params[C.GEOMETRY_DIR]) - return {k: geom_dir.joinpath(k.lower() + '.tif').as_posix() for k in C.GEOMETRY_OUTPUT_TYPES} + return { + k: geom_dir.joinpath(k.lower() + '.tif').as_posix() for k in C.GEOMETRY_OUTPUT_TYPES + } def write_config_parser_file(conf: ConfigParser, output_conf_file: Union[str, Path]): - """replacement function for write_config_file which uses dict instead of a ConfigParser instance""" - with open(output_conf_file, 'w') as configfile: + """ + Writes a config to a config file. Similar to write_config_file but for ConfigParser. + + :param ConfigParser conf: The config parser instance to write to file. + :param str output_conf_file: output file name + """ + with open(output_conf_file, 'w', encoding="utf-8") as configfile: conf.write(configfile) @@ -417,27 +493,27 @@ def write_config_file(params, output_conf_file): :return: config file :rtype: list """ - with open(output_conf_file, 'w') as f: - for k, v in params.items(): - if v is not None: - if k == 'correct': - f.write(''.join(['[', k, ']' ':\t', '', '\n'])) + with open(output_conf_file, 'w', encoding="utf-8") as f: + for key, val in params.items(): + if val is not None: + if key == 'correct': + f.write(''.join(['[', key, ']', ':\t', '', '\n'])) f.write(''.join(['steps = ', '\n'])) - for vv in v: - f.write(''.join(['\t' + str(vv), '\n'])) - elif isinstance(v, list): + for i in val: + f.write(''.join(['\t' + str(i), '\n'])) + elif isinstance(val, list): continue else: - if isinstance(v, MultiplePaths): - if v.unwrapped_path is None: - vv = v.converted_path + if isinstance(val, MultiplePaths): + if val.unwrapped_path is None: + path = val.converted_path else: - vv = v.unwrapped_path + path = val.unwrapped_path else: - vv = v - f.write(''.join([k, ':\t', str(vv), '\n'])) + path = val + f.write(''.join([key, ':\t', str(path), '\n'])) else: - f.write(''.join([k, ':\t', '', '\n'])) + f.write(''.join([key, ':\t', '', '\n'])) def parse_namelist(nml): @@ -449,7 +525,7 @@ def parse_namelist(nml): :return: list of interferogram file names :rtype: list """ - with open(nml) as f_in: + with open(nml, encoding="utf-8") as f_in: lines = [line.rstrip() for line in f_in] return filter(None, lines) @@ -458,4 +534,3 @@ class ConfigException(Exception): """ Default exception class for configuration errors. """ - diff --git a/pyrate/constants.py b/pyrate/constants.py index 8898550e4..46f703b8d 100644 --- a/pyrate/constants.py +++ b/pyrate/constants.py @@ -1,14 +1,14 @@ -import os import re from pathlib import Path import numpy as np +from pyrate.core.mpiops import comm PYRATEPATH = Path(__file__).parent.parent __version__ = "0.6.0" CLI_DESCRIPTION = """ -PyRate workflow: +PyRate workflow: Step 1: conv2tif Step 2: prepifg @@ -21,8 +21,6 @@ more details. """ -from pyrate.core.mpiops import comm - NO_OF_PARALLEL_PROCESSES = comm.Get_size() CONV2TIF = 'conv2tif' @@ -107,13 +105,20 @@ NO_DATA_VALUE = 'noDataValue' #: FLOAT; No data averaging threshold for prepifg NO_DATA_AVERAGING_THRESHOLD = 'noDataAveragingThreshold' -# BOOL (1/2/3); Re-project data from Line of sight, 1 = vertical, 2 = horizontal, 3 = no conversion +# INT (1/2/3); Re-project data from Line of sight +# 1 = vertical +# 2 = horizontal +# 3 = no conversion # REPROJECTION = 'prjflag' # NOT CURRENTLY USED #: BOOL (0/1): Convert no data values to Nan NAN_CONVERSION = 'nan_conversion' # Prepifg parameters -#: BOOL (1/2/3/4); Method for cropping interferograms, 1 = minimum overlapping area (intersection), 2 = maximum area (union), 3 = customised area, 4 = all ifgs already same size +#: INT (1/2/3/4); Method for cropping interferograms +# 1 = minimum overlapping area (intersection) +# 2 = maximum area (union) +# 3 = customised area +# 4 = all ifgs already same size IFG_CROP_OPT = 'ifgcropopt' #: INT; Multi look factor for interferogram preparation in x dimension IFG_LKSX = 'ifglksx' @@ -141,9 +146,12 @@ REFNY = "refny" #: INT; Dimension of reference pixel search window (in number of pixels) REF_CHIP_SIZE = 'refchipsize' -#: FLOAT; Minimum fraction of observations required in search window for pixel to be a viable reference pixel +#: FLOAT; Minimum fraction of observations required in search window for pixel to be +# a viable reference pixel REF_MIN_FRAC = 'refminfrac' -#: BOOL (1/2); Reference phase estimation method (1: median of the whole interferogram, 2: median within the window surrounding the reference pixel) +#: INT (1/2); Reference phase estimation method +# 1: median of the whole interferogram +# 2: median within the window surrounding the reference pixel) REF_EST_METHOD = 'refest' MAXVAR = 'maxvar' @@ -166,7 +174,8 @@ #: STR; Name of the file list containing the pool of available baseline files BASE_FILE_LIST = 'basefilelist' -#: STR; Name of the file containing the GAMMA lookup table between lat/lon and radar coordinates (row/col) +#: STR; Name of the file containing the GAMMA lookup table between lat/lon and +# radar coordinates (row/col) LT_FILE = 'ltfile' # atmospheric error correction parameters NOT CURRENTLY USED @@ -190,9 +199,12 @@ # orbital error correction/parameters #: BOOL (1/0); Perform orbital error correction (1: yes, 0: no) ORBITAL_FIT = 'orbfit' -#: BOOL (1/2); Method for orbital error correction (1: independent, 2: network) +#: INT (1/2); Method for orbital error correction (1: independent, 2: network) ORBITAL_FIT_METHOD = 'orbfitmethod' -#: BOOL (1/2/3) Polynomial order of orbital error model (1: planar in x and y - 2 parameter model, 2: quadratic in x and y - 5 parameter model, 3: quadratic in x and cubic in y - part-cubic 6 parameter model) +#: INT (1/2/3) Polynomial order of orbital error model +# 1: planar in x and y - 2 parameter model +# 2: quadratic in x and y - 5 parameter model +# 3: quadratic in x and cubic in y - part-cubic 6 parameter model ORBITAL_FIT_DEGREE = 'orbfitdegrees' #: INT; Multi look factor for orbital error calculation in x dimension ORBITAL_FIT_LOOKS_X = 'orbfitlksx' @@ -205,7 +217,9 @@ ORBFIT_INTERCEPT = 'orbfitintercept' # Stacking parameters -#: FLOAT; Threshold ratio between 'model minus observation' residuals and a-priori observation standard deviations for stacking estimate acceptance (otherwise remove furthest outlier and re-iterate) +#: FLOAT; Threshold ratio between 'model minus observation' residuals and a-priori observation +# standard deviations for stacking estimate acceptance +# (otherwise remove furthest outlier and re-iterate) LR_NSIG = 'nsig' #: INT; Number of required observations per pixel for stacking to occur LR_PTHRESH = 'pthr' @@ -230,7 +244,8 @@ SLPF_CUTOFF = 'slpfcutoff' #: INT (1/0); Do spatial interpolation at NaN locations (1 for interpolation, 0 for zero fill) SLPF_NANFILL = 'slpnanfill' -#: #: STR; Method for spatial interpolation (one of: linear, nearest, cubic), only used when slpnanfill=1 +#: #: STR; Method for spatial interpolation (one of: linear, nearest, cubic), +# only used when slpnanfill=1 SLPF_NANFILL_METHOD = 'slpnanfill_method' # DEM error correction parameters @@ -251,7 +266,10 @@ # tsinterp is automatically assigned in the code; not needed in conf file # TIME_SERIES_INTERP = 'tsinterp' -#: BOOL (0/1/2); Use parallelisation/Multi-threading (0: in serial, 1: in parallel by rows, 2: in parallel by pixel) +#: INT (0/1/2); Use parallelisation/Multi-threading +# 0: in serial +# 1: in parallel by rows +# 2: in parallel by pixel PARALLEL = 'parallel' #: INT; Number of processes for multi-threading PROCESSES = 'processes' @@ -271,7 +289,9 @@ PART_CUBIC: 'PART CUBIC'} # geometry outputs -GEOMETRY_OUTPUT_TYPES = ['rdc_azimuth', 'rdc_range', 'look_angle', 'incidence_angle', 'azimuth_angle', 'range_dist'] +GEOMETRY_OUTPUT_TYPES = [ + 'rdc_azimuth', 'rdc_range', 'look_angle', 'incidence_angle', 'azimuth_angle', 'range_dist' +] # sign convention for phase data SIGNAL_POLARITY = 'signal_polarity' @@ -302,4 +322,3 @@ GEOMETRY_DIR = 'geometry_dir' TIMESERIES_DIR = 'timeseries_dir' VELOCITY_DIR = 'velocity_dir' - diff --git a/pyrate/conv2tif.py b/pyrate/conv2tif.py index 1cf72ec2f..082367a05 100644 --- a/pyrate/conv2tif.py +++ b/pyrate/conv2tif.py @@ -14,15 +14,16 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -This Python script converts ROI_PAC or GAMMA format input interferograms +This Python script converts ROI_PAC or GAMMA format input interferograms into geotiff format files """ # -*- coding: utf-8 -*- import os from typing import Tuple, List +from pathlib import Path + from joblib import Parallel, delayed import numpy as np -from pathlib import Path import pyrate.constants as C from pyrate.core.prepifg_helper import PreprocessError @@ -52,7 +53,7 @@ def main(params): if params[C.PROCESSOR] == 2: # if geotif log.warning("'conv2tif' step not required for geotiff!") - return + return None mpi_vs_multiprocess_logging("conv2tif", params) @@ -80,7 +81,7 @@ def do_geotiff(unw_paths: List[MultiplePaths], params: dict) -> List[str]: parallel = params[C.PARALLEL] if parallel: - log.info("Running geotiff conversion in parallel with {} processes".format(params[C.PROCESSES])) + log.info("Running geotiff conversion in parallel with {params[C.PROCESSES]} processes") dest_base_ifgs = Parallel(n_jobs=params[C.PROCESSES], verbose=shared.joblib_log_level( C.LOG_LEVEL))( delayed(_geotiff_multiprocessing)(p, params) for p in unw_paths) @@ -99,19 +100,19 @@ def _geotiff_multiprocessing(unw_path: MultiplePaths, params: dict) -> Tuple[str processor = params[C.PROCESSOR] # roipac or gamma # Create full-res geotiff if not already on disk - if not os.path.exists(dest): - if processor == GAMMA: - header = gamma.gamma_header(unw_path.unwrapped_path, params) - elif processor == ROIPAC: - log.info("Warning: ROI_PAC support will be deprecated in a future PyRate release") - 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[ - C.NO_DATA_VALUE]) - Path(dest).chmod(0o444) # readonly output - return dest, True - else: + if os.path.exists(dest): log.warning(f"Full-res geotiff already exists in {dest}! Returning existing geotiff!") return dest, False + + if processor == GAMMA: + header = gamma.gamma_header(unw_path.unwrapped_path, params) + elif processor == ROIPAC: + log.info("Warning: ROI_PAC support will be deprecated in a future PyRate release") + 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[ + C.NO_DATA_VALUE]) + Path(dest).chmod(0o444) # readonly output + return dest, True diff --git a/pyrate/core/algorithm.py b/pyrate/core/algorithm.py index 94276a71a..2110fedb4 100644 --- a/pyrate/core/algorithm.py +++ b/pyrate/core/algorithm.py @@ -78,6 +78,8 @@ def least_squares_covariance(A, b, v): V = diag(1.0 / v.squeeze()) q, r = qr(A) # Orthogonal-triangular Decomposition efg = dot(q.T, dot(V, q)) # TODO: round it?? + # pylint: disable=invalid-sequence-index + # JUSTIFICATION: Pylint is just wrong here, even if we cast to int(n) it complains... g = efg[n:, n:] # modified to 0 indexing cd = dot(q.T, b) # q.T * b f = efg[:n, n:] # TODO: check +1/indexing @@ -147,15 +149,14 @@ def ifg_date_lookup(ifgs, date_pair): # and not in order if date_pair[0] > date_pair[1]: date_pair = date_pair[1], date_pair[0] - except: - raise ValueError("Bad date_pair arg to ifg_date_lookup()") + except Exception as error: + raise ValueError("Bad date_pair arg to ifg_date_lookup()") from error for i in ifgs: if date_pair == (i.first, i.second): return i - raise ValueError("Cannot find Ifg with first/second" - "image dates of %s" % str(date_pair)) + raise ValueError(f"Cannot find Ifg with first/second image dates of {date_pair}") def ifg_date_index_lookup(ifgs, date_pair): @@ -178,14 +179,14 @@ def ifg_date_index_lookup(ifgs, date_pair): try: if date_pair[0] > date_pair[1]: date_pair = date_pair[1], date_pair[0] - except: - raise ValueError("Bad date_pair arg to ifg_date_lookup()") + except Exception as error: + raise ValueError("Bad date_pair arg to ifg_date_lookup()") from error for i, _ in enumerate(ifgs): if date_pair == (ifgs[i].first, ifgs[i].second): return i - raise ValueError("Cannot find Ifg with first/second image dates of %s" % str(date_pair)) + raise ValueError(f"Cannot find Ifg with first/second image dates of {date_pair}") def get_epochs(ifgs: Union[Iterable, Dict]) -> Tuple[EpochList, int]: @@ -234,10 +235,10 @@ def first_second_ids(dates): """ dset = sorted(set(dates)) - return dict([(date_, i) for i, date_ in enumerate(dset)]) + return { date_:i for i, date_ in enumerate(dset) } -def factorise_integer(n, memo={}, left=2): +def factorise_integer(n, memo=None, left=2): """ Returns two factors a and b of a supplied number n such that a * b = n. The two factors are evaluated to be as close to each other in size as possible @@ -252,26 +253,26 @@ def factorise_integer(n, memo={}, left=2): :rtype: int """ n = int(n) - if (n, left) in memo: + if memo is not None and (n, left) in memo: return memo[(n, left)] if left == 1: return n, [n] i = 2 best = n - bestTuple = [n] + best_tuple = [n] while i * i <= n: if n % i == 0: rem = factorise_integer(n / i, memo, left - 1) if rem[0] + i < best: best = rem[0] + i - bestTuple = [i] + rem[1] + best_tuple = [i] + rem[1] i += 1 # handle edge case when only one processor is available - if bestTuple == [4]: + if best_tuple == [4]: return 2, 2 - if len(bestTuple) == 1: - bestTuple.append(1) + if len(best_tuple) == 1: + best_tuple.append(1) - return int(bestTuple[0]), int(bestTuple[1]) + return int(best_tuple[0]), int(best_tuple[1]) diff --git a/pyrate/core/aps.py b/pyrate/core/aps.py index 9308fab03..050f190a7 100644 --- a/pyrate/core/aps.py +++ b/pyrate/core/aps.py @@ -158,7 +158,7 @@ def _make_aps_corrections(ts_aps: np.ndarray, ifgs: List[Ifg], params: dict) -> :param params: Dictionary of PyRate configuration parameters. """ log.debug('Reconstructing interferometric observations from time series') - # get first and second image indices + # get first and second image indices _ , n = mpiops.run_once(get_epochs, ifgs) index_first, index_second = n[:len(ifgs)], n[len(ifgs):] @@ -305,7 +305,7 @@ def gaussian_spatial_filter(image: np.ndarray, cutoff: float, x_size: float, # Estimate sigma value for Gaussian kernel function in spectral domain # by converting cutoff distance to wavenumber and applying a scaling - # factor based on fixed kernel window size. + # factor based on fixed kernel window size. sigma = np.std(dist) * (1 / cutoff) # Calculate kernel weights wgt = _kernel(dist, sigma) @@ -333,7 +333,7 @@ def temporal_high_pass_filter(tsincr: np.ndarray, epochlist: EpochList, log.info('Applying temporal high-pass filter') threshold = params[C.TLPF_PTHR] cutoff_day = params[C.TLPF_CUTOFF] - if cutoff_day < 1 or type(cutoff_day) != int: + if cutoff_day < 1 or not isinstance(cutoff_day, int): raise ValueError(f'tlpf_cutoff must be an integer greater than or ' f'equal to 1 day. Value provided = {cutoff_day}') diff --git a/pyrate/core/covariance.py b/pyrate/core/covariance.py index 95ca52a2a..d3f519962 100644 --- a/pyrate/core/covariance.py +++ b/pyrate/core/covariance.py @@ -109,7 +109,8 @@ def _save_cvd_data(acg, r_dist, ifg_path, outdir): Function to save numpy array of autocorrelation data to disk """ data = np.column_stack((acg, r_dist)) - data_file = join(outdir, 'cvd_data_{b}.npy'.format(b=basename(ifg_path).split('.')[0])) + suffix = basename(ifg_path).split('.')[0] + data_file = join(outdir, f'cvd_data_{suffix}.npy') np.save(file=data_file, arr=data) @@ -188,12 +189,11 @@ def cvd_from_phase(phase, ifg, r_dist, calc_alpha, save_acg=False, params=None): alphaguess = 2 / (maxbin * bin_width) alpha = fmin(_pendiffexp, x0=alphaguess, args=(cvdav,), disp=False, xtol=1e-6, ftol=1e-6) - log.debug("1st guess alpha {}, converged " - "alpha: {}".format(alphaguess, alpha)) + log.debug(f"1st guess alpha {alphaguess}, converged alpha: {alpha}") # maximum variance usually at the zero lag: max(acg[:len(r_dist)]) return np.max(acg), alpha[0] # alpha unit 1/km - else: - return np.max(acg), None + + return np.max(acg), None class RDist(): @@ -280,7 +280,6 @@ def get_vcmt(ifgs, maxvar): if isinstance(ifgs, dict): ifgs = {k: v for k, v in ifgs.items() if isinstance(v, PrereadIfg)} ifgs = OrderedDict(sorted(ifgs.items())) - # pylint: disable=redefined-variable-type ifgs = ifgs.values() nifgs = len(ifgs) @@ -330,9 +329,14 @@ def _get_r_dist(ifg_path): r_dist = mpiops.run_once(_get_r_dist, ifg_paths[0]) prcs_ifgs = mpiops.array_split(list(enumerate(ifg_paths))) process_maxvar = {} + n_prcs = len(prcs_ifgs) + n_ifgs = len(ifg_paths) + for n, i in prcs_ifgs: - log.debug(f'Calculating maxvar for {n} of process ifgs {len(prcs_ifgs)} of total {len(ifg_paths)}') - process_maxvar[int(n)] = cvd(i, params, r_dist, calc_alpha=True, write_vals=True, save_acg=True)[0] + log.debug(f'Calculating maxvar for {n} of process ifgs {n_prcs} of total {n_ifgs}') + val = cvd(i, params, r_dist, calc_alpha=True, write_vals=True, save_acg=True)[0] + process_maxvar[int(n)] = val + maxvar_d = shared.join_dicts(mpiops.comm.allgather(process_maxvar)) maxvar = [v[1] for v in sorted(maxvar_d.items(), key=lambda s: s[0])] diff --git a/pyrate/core/dem_error.py b/pyrate/core/dem_error.py index 660f33f62..9ce88e575 100644 --- a/pyrate/core/dem_error.py +++ b/pyrate/core/dem_error.py @@ -14,14 +14,15 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -This Python module implements the calculation of correction for residual topographic effects (a.k.a. as DEM errors). +This Python module implements the calculation of correction for residual +topographic effects (a.k.a. as DEM errors). """ # pylint: disable=invalid-name, too-many-locals, too-many-arguments import os -import numpy as np from typing import Tuple, Optional from os.path import join from pathlib import Path +import numpy as np import pyrate.constants as C from pyrate.core import geometry, shared, mpiops, ifgconstants as ifc @@ -64,10 +65,11 @@ def dem_error_calc_wrapper(params: dict) -> None: log.warning("No baseline files supplied: DEM error has not been computed") return - # todo: subtract other corrections (e.g. orbital) from displacement phase before estimating the DEM error + # TODO: subtract other corrections (e.g. orbital) from displacement + # phase before estimating the DEM error - # the following code is a quick way to do the bperp calculation, but is not identical to the GAMMA output - # where the near range of the first SLC is used for each pair. + # the following code is a quick way to do the bperp calculation, but is not identical to + # the GAMMA output where the near range of the first SLC is used for each pair. # calculate look angle for interferograms (using the Near Range of the first SLC) # look_angle = geometry.calc_local_geometry(ifg0, None, rg, lon, lat, params) # bperp = geometry.calc_local_baseline(ifg0, az, look_angle) @@ -116,26 +118,34 @@ def _process_dem_error_per_tile(tile: Tile, params: dict) -> None: lat_parts, lon_parts, dem_parts) log.debug(f"Calculating DEM error for tile {tile.index} during DEM error correction") + tmp_dir = Path(params[C.TMPDIR]) + # mst_tile = np.load(Configuration.mst_path(params, tile.index)) # calculate the DEM error estimate and the correction values for each IFG - # current implementation uses the look angle and range distance matrix of the primary SLC in the last IFG - # todo: check the impact of using the same information from another SLC - dem_error, dem_error_correction, _ = calc_dem_errors(ifg_parts, bperp, look_angle, range_dist, threshold) - # dem_error contains the estimated DEM error for each pixel (i.e. the topographic change relative to the DEM) + # current implementation uses the look angle and range distance matrix of + # the primary SLC in the last IFG + + # TODO: check the impact of using the same information from another SLC + dem_error, dem_error_correction, _ = calc_dem_errors( + ifg_parts, bperp, look_angle, range_dist, threshold + ) + # dem_error contains the estimated DEM error for each pixel + # (i.e. the topographic change relative to the DEM) # size [row, col] + # dem_error_correction contains the correction value for each interferogram # size [num_ifg, row, col] # save tiled data in tmpdir - np.save(file=os.path.join(params[C.TMPDIR], 'dem_error_{}.npy'.format(tile.index)), arr=dem_error) + np.save(file=tmp_dir / f'dem_error_{tile.index}.npy', arr=dem_error) # swap the axes of 3D array to fit the style used in function assemble_tiles tmp_array = np.moveaxis(dem_error_correction, 0, -1) # new dimension is [row, col, num_ifg] # save tiled data into tmpdir - np.save(file=os.path.join(params[C.TMPDIR], 'dem_error_correction_{}.npy'.format(tile.index)), arr=tmp_array) + np.save(file=tmp_dir / f'dem_error_correction_{tile.index}.npy', arr=tmp_array) # Calculate and save the average perpendicular baseline for the tile bperp_avg = np.nanmean(bperp, axis=(1, 2), dtype=np.float64) - np.save(file=os.path.join(params[C.TMPDIR], 'bperp_avg_{}.npy'.format(tile.index)), arr=bperp_avg) + np.save(file=tmp_dir / f'bperp_avg_{tile.index}.npy', arr=bperp_avg) def _calculate_bperp_wrapper(ifg_paths: list, az_parts: np.ndarray, rg_parts: np.ndarray, @@ -161,7 +171,9 @@ def _calculate_bperp_wrapper(ifg_paths: list, az_parts: np.ndarray, rg_parts: np ifg = Ifg(ifg_path) ifg.open(readonly=True) # calculate look angle for interferograms (using the Near Range of the primary SLC) - look_angle, _, _, range_dist = geometry.calc_pixel_geometry(ifg, rg_parts, lon_parts, lat_parts, dem_parts) + look_angle, _, _, range_dist = geometry.calc_pixel_geometry( + ifg, rg_parts, lon_parts, lat_parts, dem_parts + ) bperp[ifg_num, :, :] = geometry.calc_local_baseline(ifg, az_parts, look_angle) return bperp, look_angle, range_dist @@ -169,8 +181,8 @@ def _calculate_bperp_wrapper(ifg_paths: list, az_parts: np.ndarray, rg_parts: np def calc_dem_errors(ifgs: list, bperp: np.ndarray, look_angle: np.ndarray, range_dist: np.ndarray, threshold: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ - Function to calculate the DEM error for each pixel using least-squares adjustment of phase data and - perpendicular baseline. The least-squares adjustment co-estimates the velocities. + Function to calculate the DEM error for each pixel using least-squares adjustment of phase + data and perpendicular baseline. The least-squares adjustment co-estimates the velocities. - *nrows* is the number of rows in the ifgs, and - *ncols* is the number of columns in the ifgs. @@ -180,12 +192,13 @@ def calc_dem_errors(ifgs: list, bperp: np.ndarray, look_angle: np.ndarray, range :param range_dist: Per-pixel range distance measurement :param threshold: minimum number of redundant phase values at pixel (config parameter de_pthr) :return: dem_error: estimated per-pixel dem error (nrows x ncols) - :return: dem_error_correction: DEM error correction for each pixel and interferogram (nifgs x nrows x ncols) + :return: dem_error_correction: DEM error correction for each pixel and interferogram + (nifgs x nrows x ncols) :return: vel: velocity estimate for each pixel (nrows x ncols) """ ifg_data, mst, ncols, nrows, ifg_time_span = _per_tile_setup(ifgs) if threshold < 4: - msg = f"pixel threshold too low (i.e. <4) resulting in non-redundant DEM error estimation" + msg = "pixel threshold too low (i.e. <4) resulting in non-redundant DEM error estimation" raise DEMError(msg) # pixel-by-pixel calculation # preallocate empty arrays for results @@ -204,25 +217,30 @@ def calc_dem_errors(ifgs: list, bperp: np.ndarray, look_angle: np.ndarray, range # If NaN: skip pixel if np.isnan(y).any() or np.isnan(bperp_pix).any(): continue - # using the actual geometry of a particular IFG would be possible but is likely not signif. different + # using the actual geometry of a particular IFG would be possible but is likely not + # significantly different... # geom = bperp_pix / (range_dist[row, col] * np.sin(look_angle[row, col])) time_span = ifg_time_span[sel] # new covariance matrix using actual number of observations m = len(sel) - # Design matrix of least-squares system, velocities are co-estimated as done in StaMPS or MintPy + # Design matrix of least-squares system, velocities are co-estimated as + # done in StaMPS or MintPy A = np.column_stack((np.ones(m), bperp_pix, time_span)) # solve ordinary least-squares system using numpy.linalg.lstsq function - xhat, res, rnk, s = np.linalg.lstsq(A, y, rcond=None) + xhat, _, _, _ = np.linalg.lstsq(A, y, rcond=None) # dem error estimate for the pixel is the second parameter (cf. design matrix) dem_error[row][col] = xhat[1] # velocity estimate for the pixel is the third parameter (cf. design matrix) vel[row][col] = xhat[2] - # calculate correction value for each IFG by multiplying the least-squares estimate with the Bperp value + # calculate correction value for each IFG by multiplying the + # least-squares estimate with the Bperp value dem_error_correction = np.multiply(dem_error, bperp) # calculate metric difference to DEM by multiplying the estimate with the per-pixel geometry # (i.e. range distance and look angle, see Eq. (2.4.12) in Hanssen (2001)) - # also scale by -0.001 since the phase observations are in mm with positive values away from the sensor + # + # also scale by -0.001 since the phase observations are in mm with + # positive values away from the sensor dem_error = np.multiply(dem_error, np.multiply(range_dist, np.sin(look_angle))) * (-0.001) return dem_error, dem_error_correction, vel @@ -233,7 +251,8 @@ def _per_tile_setup(ifgs: list) -> Tuple[np.ndarray, np.ndarray, int, int, np.nd Convenience function for setting up DEM error computation parameters :param ifgs: list of interferogram class objects. :return: ifg_data: phase observations for each pixel and interferogram - :return: mst: an array of booleans representing valid ifg connections (i.e. the minimum spanning tree) + :return: mst: an array of booleans representing valid ifg connections + (i.e. the minimum spanning tree) :return: ncols: number of columns :return: nrows: number of rows :return: ifg_time_span: date difference for each interferogram @@ -258,7 +277,9 @@ def _per_tile_setup(ifgs: list) -> Tuple[np.ndarray, np.ndarray, int, int, np.nd def _write_dem_errors(ifg_paths: list, params: dict, preread_ifgs: dict) -> None: """ - Convenience function for writing DEM error (one file) and DEM error correction for each IFG to disc + Convenience function for writing DEM error (one file) and DEM error correction + for each IFG to disc. + :param ifg_paths: List of interferogram class objects. :param params: Dictionary of PyRate configuration parameters. :param preread_ifgs: Dictionary of interferogram metadata. @@ -290,8 +311,10 @@ def _write_dem_errors(ifg_paths: list, params: dict, preread_ifgs: dict) -> None ifg.open() shared.nan_and_mm_convert(ifg, params) # ensure we have phase data in mm # read dem error correction file from tmpdir - dem_error_correction_ifg = assemble_tiles(shape, params[C.TMPDIR], tiles, out_type='dem_error_correction', - index=idx) + dem_error_correction_ifg = assemble_tiles( + shape, params[C.TMPDIR], tiles, + out_type='dem_error_correction', index=idx + ) # calculate average bperp value across all tiles for the ifg bperp_val = np.nanmean(bperp[:, idx]) dem_error_correction_on_disc = MultiplePaths.dem_error_path(ifg.data_path, params) @@ -305,7 +328,9 @@ def _write_dem_errors(ifg_paths: list, params: dict, preread_ifgs: dict) -> None def __check_and_apply_demerrors_found_on_disc(ifg_paths: list, params: dict) -> bool: """ - Convenience function to check if DEM error correction files have already been produced in a previous run + Convenience function to check if DEM error correction files have already been + produced in a previous run + :param ifg_paths: List of interferogram class objects. :param params: Dictionary of PyRate configuration parameters. :return: bool value: True if dem error files found on disc, otherwise False. diff --git a/pyrate/core/gamma.py b/pyrate/core/gamma.py index cc0c92360..38b30bdf1 100644 --- a/pyrate/core/gamma.py +++ b/pyrate/core/gamma.py @@ -22,6 +22,7 @@ from os.path import split from pathlib import Path from datetime import date, time, timedelta +import struct import numpy as np import pyrate.constants as C @@ -30,7 +31,6 @@ from pyrate.constants import sixteen_digits_pattern, BASELINE_FILE_PATHS, BASE_FILE_DIR from pyrate.core.shared import extract_epochs_from_filename, data_format from pyrate.core.logger import pyratelogger as log -import struct # constants @@ -68,7 +68,7 @@ def _parse_header(path): """Parses all GAMMA header file fields into a dictionary""" - with open(path) as f: + with open(path, "r", encoding="utf-8") as f: text = f.read().splitlines() raw_segs = [line.split() for line in text if ':' in line] @@ -80,11 +80,17 @@ def parse_epoch_header(path): """ Returns dictionary of epoch metadata required for PyRate - :param str path: Full path to GAMMA mli.par file. Note that the mli.par is required as input since the baseline calculations require the input values valid for the GAMMA multi-looked products and also the GAMMA lookup table gives radar coordinates for the multi-looked geometry. + :param str path: Full path to GAMMA mli.par file. + Note that the mli.par is required as input since the baseline calculations require the + input values valid for the GAMMA multi-looked products and also the GAMMA lookup table + gives radar coordinates for the multi-looked geometry. :return: subset: subset of full metadata :rtype: dict """ + # pylint: disable=too-many-statements + # JUSTIFICATION: keeping the validation code together makes it the easiest to understand, + # splitting it apart is premature at this stage, would eventually make sense if it grew too big lookup = _parse_header(path) subset = _parse_date_time(lookup) @@ -177,21 +183,21 @@ def _parse_date_time(lookup): year, month, day, = [int(float(i)) for i in lookup[GAMMA_DATE][:3]] if lookup.get(GAMMA_TIME) is not None: t = lookup[GAMMA_TIME][0] - h, m, s = str(timedelta(seconds=float(t))).split(":") - hour = int(h) - min = int(m) - sec = int(s.split(".")[0]) + hour, mins, sec = str(timedelta(seconds=float(t))).split(":") + hour = int(hour) + mins = int(mins) + sec = int(sec.split(".")[0]) else: # Occasionally GAMMA header has no time information - default to midnight - hour, min, sec = 0, 0, 0 + hour, mins, sec = 0, 0, 0 elif len(lookup[GAMMA_DATE]) == 6: - year, month, day, hour, min, sec = [int(float(i)) for i in lookup[GAMMA_DATE][:6]] + year, month, day, hour, mins, sec = [int(float(i)) for i in lookup[GAMMA_DATE][:6]] else: # pragma: no cover msg = "Date and time information not complete in GAMMA headers" raise GammaException(msg) subset[ifc.FIRST_DATE] = date(year, month, day) - subset[ifc.FIRST_TIME] = time(hour, min, sec) + subset[ifc.FIRST_TIME] = time(hour, mins, sec) return subset @@ -208,7 +214,10 @@ def parse_dem_header(path): lookup = _parse_header(path) # NB: many lookup fields have multiple elements, eg ['1000', 'Hz'] - subset = {ifc.PYRATE_NCOLS: int(lookup[GAMMA_WIDTH][0]), ifc.PYRATE_NROWS: int(lookup[GAMMA_NROWS][0])} + subset = { + ifc.PYRATE_NCOLS: int(lookup[GAMMA_WIDTH][0]), + ifc.PYRATE_NROWS: int(lookup[GAMMA_NROWS][0]) + } expected = ['decimal', 'degrees'] for k in [GAMMA_CORNER_LAT, GAMMA_CORNER_LONG, GAMMA_X_STEP, GAMMA_Y_STEP]: @@ -287,7 +296,11 @@ def combine_headers(hdr0, hdr1, dem_hdr, base_hdr=None): :return: chdr: combined metadata :rtype: dict """ - if not all([isinstance(a, dict) for a in [hdr0, hdr1, dem_hdr]]): + # pylint: disable=too-many-locals,too-many-branches,too-many-statements + # JUSTIFICATION: keeping the validation code together makes it the easiest to understand, + # splitting it apart is premature at this stage, would eventually make sense if it grew too big + + if not all(isinstance(a, dict) for a in [hdr0, hdr1, dem_hdr]): raise GammaException('Header args need to be dicts') if base_hdr and not isinstance(base_hdr, dict): @@ -296,7 +309,7 @@ def combine_headers(hdr0, hdr1, dem_hdr, base_hdr=None): date0, date1 = hdr0[ifc.FIRST_DATE], hdr1[ifc.FIRST_DATE] if date0 == date1: raise GammaException("Can't combine headers for the same day") - elif date1 < date0: + if date1 < date0: raise GammaException("Wrong date order") chdr = {ifc.PYRATE_TIME_SPAN: (date1 - date0).days / ifc.DAYS_PER_YEAR, @@ -461,7 +474,9 @@ def manage_headers(dem_header_file, header_paths, baseline_paths=None): else: combined_header = combine_headers(hdrs[0], hdrs[1], dem_header) elif len(header_paths) > 2: - msg = f'There are too many parameter files for one interferogram; there should only be two. {len(header_paths)} parameter files have been given: {header_paths}.' + msg = 'There are too many parameter files for one interferogram;' \ + ' there should only be two. {len(header_paths)} parameter files' \ + f' have been given: {header_paths}.' raise GammaException(msg) else: # probably have DEM or incidence file @@ -476,7 +491,8 @@ def get_header_paths(input_file, slc_file_list): Function that matches input GAMMA file names with GAMMA header file names :param str input_file: input GAMMA image file. - :param slc_file_list: file listing the pool of available header files (GAMMA: slc.par, ROI_PAC: .rsc) + :param slc_file_list: file listing the pool of available header files + Example: slc.par (for GAMMA) or .rsc (for ROI_PAC) :return: list of matching header files :rtype: list """ @@ -490,14 +506,14 @@ def get_header_paths(input_file, slc_file_list): def gamma_header(ifg_file_path, params): """ Function to obtain combined Gamma headers for image file - + Args: ifg_file_path: Path to interferogram file to find headers for. params: PyRate parameters dictionary. Returns: A combined header dictionary containing metadata from matching - gamma headers and DEM header. + gamma headers and DEM header. """ dem_hdr_path = params[C.DEM_HEADER_FILE] header_paths = get_header_paths(ifg_file_path, params[C.HDR_FILE_LIST]) @@ -517,7 +533,6 @@ def gamma_header(ifg_file_path, params): def read_lookup_table(head, data_path, xlooks, ylooks, xmin, xmax, ymin, ymax): - # pylint: disable = too - many - statements """ Creates a copy of input lookup table file in a numpy array and applies the ifg ML factors @@ -544,7 +559,8 @@ def read_lookup_table(head, data_path, xlooks, ylooks, xmin, xmax, ymin, ymax): ifg_proc = head.meta_data[ifc.PYRATE_INSAR_PROCESSOR] # get dimensions of lookup table file - bytes_per_col, fmtstr = data_format(ifg_proc, True, ncols_lt*2) # float complex data set containing value tupels + # float complex data set containing value tuples + bytes_per_col, fmtstr = data_format(ifg_proc, True, ncols_lt*2) # check if lookup table has the correct size small_size = _check_raw_data(bytes_per_col * 2, data_path, ncols_lt, nrows_lt) @@ -576,12 +592,13 @@ def read_lookup_table(head, data_path, xlooks, ylooks, xmin, xmax, ymin, ymax): idx_start = ymin + int((ylooks - 1) / 2) row_idx = np.arange(idx_start, ymax, ylooks) - # read the binary lookup table file and save the range/azimuth value pair for each position in the cropped and - # multi-looked data set + # read the binary lookup table file and save the range/azimuth value pair for each + # position in the cropped and multi-looked data set log.debug(f"Reading lookup table file {data_path}") with open(data_path, 'rb') as f: for y in range(nrows_lt): # loop through all lines in file - # this could potentially be made quicker by skipping unwanted bytes in the f.read command? + # this could potentially be made quicker by skipping unwanted bytes + # in the f.read command? data = struct.unpack(fmtstr, f.read(row_bytes)) # but only read data from lines in row index: if y in row_idx: @@ -606,8 +623,10 @@ def _check_raw_data(bytes_per_col, data_path, ncols, nrows): # test data set doesn't currently fit the lookup table size, stop further calculation # todo: delete this if statement once a new test data set has been introduced return True - else: - raise GammaException(msg % (data_path, size, act_size)) + + raise GammaException(msg % (data_path, size, act_size)) + + return False class GammaException(Exception): @@ -635,10 +654,13 @@ def baseline_paths_for(path: str, params: dict) -> str: _, filename = split(path) try: epoch = re.search(sixteen_digits_pattern, filename).group(0) - except: # catch cases where filename does not have two epochs, e.g. DEM file + # pylint: disable=bare-except + # JUSTIFICATION: intentionally catch when filename does not have two epochs, e.g. DEM file + except: return None - base_file_paths = [f.unwrapped_path for f in params[BASELINE_FILE_PATHS] if epoch in f.unwrapped_path] + base_file_paths = params[BASELINE_FILE_PATHS] + base_file_paths = [f.unwrapped_path for f in base_file_paths if epoch in f.unwrapped_path] if len(base_file_paths) > 1: raise ConfigException(f"'{BASE_FILE_DIR}': found more than one baseline " diff --git a/pyrate/core/geometry.py b/pyrate/core/geometry.py index c0d11b877..b360d8317 100644 --- a/pyrate/core/geometry.py +++ b/pyrate/core/geometry.py @@ -14,14 +14,16 @@ # See the License for the specific language governing permissions and # limitations under the License. """ -This Python module implements the calculation and output of the per-pixel vector of the radar viewing geometry -(i.e. local angles, incidence angles and azimuth angles) as well as the calculation of per-pixel baseline values -used for correcting interferograms for residual topographic effects (a.k.a. DEM errors). +This Python module implements the calculation and output of the per-pixel vector of the radar +viewing geometry (i.e. local angles, incidence angles and azimuth angles) as well as the +calculation of per-pixel baseline values used for correcting interferograms for residual +topographic effects (a.k.a. DEM errors). """ # pylint: disable=invalid-name, too-many-locals, too-many-arguments -import numpy as np from typing import Tuple, Union +import numpy as np + import pyrate.constants as C from pyrate.core import ifgconstants as ifc from pyrate.core.gamma import read_lookup_table @@ -71,13 +73,14 @@ def calc_radar_coords(ifg: Ifg, params: dict, xmin: int, xmax: int, lookup_table = params[C.LT_FILE] if lookup_table is None: - msg = f"No lookup table file supplied: Geometry cannot be computed" + msg = "No lookup table file supplied: Geometry cannot be computed" raise FileNotFoundError(msg) # PyRate IFG multi-looking factors ifglksx = params[C.IFG_LKSX] ifglksy = params[C.IFG_LKSY] - # transform float lookup table file to np array, min/max pixel coordinates are required for cropping + # transform float lookup table file to np array, + # min/max pixel coordinates are required for cropping lt_az, lt_rg = read_lookup_table(ifg, lookup_table, ifglksx, ifglksy, xmin, xmax, ymin, ymax) # replace 0.0 with NaN lt_az[lt_az == 0.0] = np.nan @@ -86,42 +89,63 @@ def calc_radar_coords(ifg: Ifg, params: dict, xmin: int, xmax: int, return lt_az, lt_rg -def get_sat_positions(lat: np.ndarray, lon: np.ndarray, look_angle: np.ndarray, inc_angle: np.ndarray, - heading: np.float64, look_dir: np.float64) -> Tuple[np.ndarray, np.ndarray]: +def get_sat_positions( + lat: np.ndarray, + lon: np.ndarray, + look_angle: np.ndarray, + inc_angle: np.ndarray, + heading: np.float64, + look_dir: np.float64 +) -> Tuple[np.ndarray, np.ndarray]: """ Function to calculate the lon/lat position of the satellite for each pixel. :param lat: Ground latitude for each pixel (decimal degrees). :param lon: Ground point longitude for each pixel (decimal degrees). :param look_angle: Look angle (between nadir and look vector) for each pixel (radians). - :param inc_angle: Local incidence angle (between vertical and look vector) for each pixel (radians). + :param inc_angle: + Local incidence angle (between vertical and look vector) for each pixel (radians). :param heading: Satellite flight heading (radians). :param look_dir: Look direction w.r.t. satellite heading; +ve = right looking (radians). :return: sat_lat: Satellite position latitude for each pixel (decimal degrees). :return: sat_lon: Satellite position longitude for each pixel (decimal degrees). """ - # note that the accuracy of satellite lat/lon positions code could be improved by calculating satellite positions - # for each azimuth row from orbital state vectors given in .par file, using the following workflow: + # note that the accuracy of satellite lat/lon positions code could be improved by + # calculating satellite positions for each azimuth row from orbital state vectors given + # in .par file, using the following workflow: # 1. read orbital state vectors and start/stop times from mli.par # 2. for each pixel get the corresponding radar row (from matrix az) # 3. get the corresponding radar time for that row (using linear interpolation) - # 4. calculate the satellite XYZ position for that time by interpolating the time and velocity state vectors + # 4. calculate the satellite XYZ position for that time by interpolating the + # time and velocity state vectors # angle at the Earth's center between se and re epsilon = np.pi - look_angle - (np.pi - inc_angle) - # azimuth of satellite look vector (satellite heading + look direction (+90 deg for right-looking SAR) + # azimuth of satellite look vector (satellite heading + look direction) + # (note: +90 deg for right-looking SAR) sat_azi = heading + look_dir - # the following equations are adapted from Section 4.4 (page 4-16) in EARTH-REFERENCED AIRCRAFT NAVIGATION AND - # SURVEILLANCE ANALYSIS (https://ntlrepository.blob.core.windows.net/lib/59000/59300/59358/DOT-VNTSC-FAA-16-12.pdf) - sat_lon = np.divide(np.arcsin(-(np.multiply(np.sin(epsilon), np.sin(sat_azi)))), np.cos(lat)) + lon # Eq. 103 - temp = np.multiply(np.divide(np.cos(0.5 * (sat_azi + sat_lon - lon)), np.cos(0.5 * (sat_azi - sat_lon + lon))), \ - np.tan(0.5 * (np.pi / 2 + lat - epsilon))) # Eq. 104 + # the following equations are adapted from Section 4.4 (page 4-16) + # in EARTH-REFERENCED AIRCRAFT NAVIGATION AND SURVEILLANCE ANALYSIS: + # https://ntlrepository.blob.core.windows.net/lib/59000/59300/59358/DOT-VNTSC-FAA-16-12.pdf + + # Eq. 103 + sat_lon = (np.arcsin(-(np.multiply(np.sin(epsilon), np.sin(sat_azi)))) / np.cos(lat)) + lon + # Eq. 104 + temp = np.multiply( + np.divide(np.cos(0.5 * (sat_azi + sat_lon - lon)), np.cos(0.5 * (sat_azi - sat_lon + lon))), + np.tan(0.5 * (np.pi / 2 + lat - epsilon)) + ) sat_lat = -np.pi / 2 + 2 * np.arctan(temp) return np.real(sat_lat), np.real(sat_lon) -def calc_pixel_geometry(ifg: Union[Ifg, IfgPart], rg: np.ndarray, lon: np.ndarray, lat: np.ndarray, - dem_height: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: +def calc_pixel_geometry( + ifg: Union[Ifg, IfgPart], + rg: np.ndarray, + lon: np.ndarray, + lat: np.ndarray, + dem_height: np.ndarray +) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Function to calculate angular satellite to ground geometries and distance for each pixel. :param ifg: pyrate.core.shared.Ifg Class object. @@ -129,8 +153,10 @@ def calc_pixel_geometry(ifg: Union[Ifg, IfgPart], rg: np.ndarray, lon: np.ndarra :param lon: Longitude for each pixel (decimal degrees). :param lat: Latitude for each pixel (decimal degrees). :param dem_height: Height from DEM for each pixel (metres). - :return: look_angle: Look angle (between nadir and look vector) for each pixel (radians). - :return: incidence_angle: Local incidence angle (between vertical and look vector) for each pixel (radians). + :return: look_angle: + Look angle (between nadir and look vector) for each pixel (radians). + :return: incidence_angle: + Local incidence angle (between vertical and look vector) for each pixel (radians). :return: azimuth_angle: Geodetic azimuth for each pixel (radians). :return: range_dist: Distance from satellite to ground for each pixel (metres). """ @@ -142,7 +168,7 @@ def calc_pixel_geometry(ifg: Union[Ifg, IfgPart], rg: np.ndarray, lon: np.ndarra near_range = float(ifg.meta_data[ifc.PYRATE_NEAR_RANGE_METRES]) rps = float(ifg.meta_data[ifc.PYRATE_RANGE_PIX_METRES]) heading = float(ifg.meta_data[ifc.PYRATE_HEADING_DEGREES]) - # direction of look vector w.r.t. satellite heading. + # direction of look vector w.r.t. satellite heading. # Gamma convention: +ve = right; -ve = left. look_dir = float(ifg.meta_data[ifc.PYRATE_AZIMUTH_DEGREES]) @@ -161,16 +187,19 @@ def calc_pixel_geometry(ifg: Union[Ifg, IfgPart], rg: np.ndarray, lon: np.ndarra # look angle at pixel ij -> law of cosines in "satellite - Earth centre - ground pixel" triangle # see e.g. Section 2 in https://www.cs.uaf.edu/~olawlor/ref/asf/sar_equations_2006_08_17.pdf - look_angle = np.arccos(np.divide(se ** 2 + np.square(range_dist) - np.square(re), 2 * se * range_dist)) + look_angle = np.divide(se ** 2 + np.square(range_dist) - np.square(re), 2 * se * range_dist) + look_angle = np.arccos(look_angle) # add per-pixel height to the earth radius to obtain a more accurate ground pixel position for # incidence angle calculation re = re + dem_height - # incidence angle at pixel ij -> law of cosines in "satellite - Earth centre - ground pixel" triangle + # incidence angle at pixel ij: + # law of cosines in "satellite - Earth centre - ground pixel" triangle # see e.g. Section 2 in https://www.cs.uaf.edu/~olawlor/ref/asf/sar_equations_2006_08_17.pdf - incidence_angle = np.pi - np.arccos(np.divide(np.square(range_dist) + np.square(re) - se ** 2, - 2 * np.multiply(range_dist, re))) + incidence_angle = np.pi - np.arccos( + np.divide(np.square(range_dist) + np.square(re) - se ** 2, 2 * np.multiply(range_dist, re)) + ) # calculate satellite positions for each pixel sat_lat, sat_lon = get_sat_positions(lat, lon, look_angle, incidence_angle, heading, look_dir) @@ -203,14 +232,18 @@ def calc_local_baseline(ifg: Ifg, az: np.ndarray, look_angle: np.ndarray) -> np. baserate_N = float(ifg.meta_data[ifc.PYRATE_BASELINE_RATE_N]) # calculate per pixel baseline vectors across track (C) and normal to the track (N) - mean_az = az_n / 2 - 0.5 # mean azimuth line - prf = prf / az_looks # Pulse Repetition Frequency needs to be adjusted according to GAMMA azimuth looks + + # mean azimuth line + mean_az = az_n / 2 - 0.5 + # Pulse Repetition Frequency needs to be adjusted according to GAMMA azimuth looks + prf = prf / az_looks base_C_local = base_C + baserate_C * (az - mean_az) / prf base_N_local = base_N + baserate_N * (az - mean_az) / prf # calculate the per-pixel perpendicular baseline (see Eq. 3.5 in Baehr, 2012 available here: # http://www.dgk.badw.de/fileadmin/user_upload/Files/DGK/docs/c-719.pdf) - bperp = np.multiply(base_C_local, np.cos(look_angle)) - np.multiply(base_N_local, np.sin(look_angle)) + bperp = np.multiply(base_C_local, np.cos(look_angle)) + bperp -= np.multiply(base_N_local, np.sin(look_angle)) return bperp @@ -245,10 +278,13 @@ def vincinv(lat1: np.ndarray, lon1: np.ndarray, lat2: np.ndarray, lon2: np.ndarr omega = lon # Iterate until the change in lambda, lambda_sigma, is insignificant # (< 1e-12) or after 1000 iterations have been completed - for i in range(1000): + for _ in range(1000): # Eq. 74 + # pylint: disable=line-too-long + # JUSTIFICATION: slightly over but a lot easier to read this way. sin_sigma = np.sqrt( - (np.cos(u2) * np.sin(lon)) ** 2 + (np.cos(u1) * np.sin(u2) - np.sin(u1) * np.cos(u2) * np.cos(lon)) ** 2) + (np.cos(u2) * np.sin(lon)) ** 2 + (np.cos(u1) * np.sin(u2) - np.sin(u1) * np.cos(u2) * np.cos(lon)) ** 2 + ) # Eq. 75 cos_sigma = np.sin(u1) * np.sin(u2) + np.cos(u1) * np.cos(u2) * np.cos(lon) # Eq. 76 @@ -260,6 +296,8 @@ def vincinv(lat1: np.ndarray, lon1: np.ndarray, lat2: np.ndarray, lon2: np.ndarr # Eq. 79 c = (f / 16) * np.cos(alpha) ** 2 * (4 + f * (4 - 3 * np.cos(alpha) ** 2)) # Eq. 80 + # pylint: disable=line-too-long + # JUSTIFICATION: slightly over but a lot easier to read this way. new_lon = omega + (1 - c) * f * np.sin(alpha) * ( sigma + c * np.sin(sigma) * (cos_two_sigma_m + c * np.cos(sigma) * (-1 + 2 * (cos_two_sigma_m ** 2))) ) diff --git a/pyrate/core/logger.py b/pyrate/core/logger.py index 7b3752b6a..b749433b1 100644 --- a/pyrate/core/logger.py +++ b/pyrate/core/logger.py @@ -25,24 +25,26 @@ pyratelogger = logging.getLogger(__name__) -formatter = logging.Formatter("%(asctime)s %(module)s:%(lineno)d %(process)d %(levelname)s " + str(rank) + "/" + str( - size-1)+" %(message)s", "%H:%M:%S") - +formatter = logging.Formatter( + f"%(asctime)s %(module)s:%(lineno)d %(process)d %(levelname)s {rank}/{size-1} %(message)s", + "%H:%M:%S" +) def configure_stage_log(verbosity, step_name, log_file_name='pyrate.log.'): - log_file_name = run_once(str.__add__, log_file_name, step_name + '.' + datetime.now().isoformat()) + timestamp = datetime.now().isoformat() + log_file_name = run_once(str.__add__, log_file_name, step_name + '.' + timestamp) - ch = MPIStreamHandler() - ch.setLevel(verbosity) - ch.setFormatter(formatter) + stream_handler = MPIStreamHandler() + stream_handler.setLevel(verbosity) + stream_handler.setFormatter(formatter) - fh = logging.FileHandler(log_file_name) - fh.setLevel(verbosity) - fh.setFormatter(formatter) + file_handler = logging.FileHandler(log_file_name) + file_handler.setLevel(verbosity) + file_handler.setFormatter(formatter) - pyratelogger.addHandler(ch) - pyratelogger.addHandler(fh) + pyratelogger.addHandler(stream_handler) + pyratelogger.addHandler(file_handler) def warn_with_traceback(message, category, filename, lineno, line=None): diff --git a/pyrate/core/mpiops.py b/pyrate/core/mpiops.py index 4584b13be..051042ff6 100644 --- a/pyrate/core/mpiops.py +++ b/pyrate/core/mpiops.py @@ -18,13 +18,11 @@ """ # pylint: disable=no-member # pylint: disable=invalid-name -import logging -import pickle from typing import Callable, Any, Iterable import numpy as np -"""For MPI compatibility""" +# For MPI compatibility try: from mpi4py import MPI MPI_INSTALLED = True @@ -41,6 +39,9 @@ # the rank of the node. rank = comm.Get_rank() except ImportError: + # pylint: disable=missing-function-docstring,missing-class-docstring,unused-argument + # JUSTIFICATION: This is a mock for a thirdparty API which has documentation already + MPI_INSTALLED = False size = 1 rank = 0 @@ -87,10 +88,13 @@ def bcast(arr, root=0): class MPIException(Exception): - pass + """Base class for MPI related exceptions""" def validate_mpi(): + """ + Validates that MPI has been installed in the running environment, raising an exception if not. + """ if not MPI_INSTALLED: raise MPIException("MPI needs to be installed in order to use this module") @@ -113,6 +117,7 @@ def run_once(f: Callable, *args, **kwargs) -> Any: # and these processes never exit, i.e., this MPI call hangs in case process 0 errors. try: f_result = f(*args, **kwargs) + # pylint: disable=broad-except except Exception as e: f_result = e else: @@ -120,10 +125,10 @@ def run_once(f: Callable, *args, **kwargs) -> Any: result = comm.bcast(f_result, root=0) if isinstance(result, Exception): raise result - else: - return result - else: - return f(*args, **kwargs) + + return result + + return f(*args, **kwargs) def array_split(arr: Iterable, process: int = None) -> Iterable: @@ -145,11 +150,13 @@ def array_split(arr: Iterable, process: int = None) -> Iterable: def sum_vars(x, y, dtype): + """Simple wrapper for `numpy.sum()` for usage with MPI""" s = np.sum([x, y], axis=0) return s def sum_axis_0(x, y, dtype): + """Simple wrapper for `numpy.sum(numpy.stack())` for usage with MPI""" s = np.sum(np.stack((x, y)), axis=0) return s diff --git a/pyrate/core/mst.py b/pyrate/core/mst.py index 0e44fef82..75d9c6cfd 100644 --- a/pyrate/core/mst.py +++ b/pyrate/core/mst.py @@ -19,7 +19,6 @@ functionality for selecting interferometric observations. """ # pylint: disable=invalid-name -from pathlib import Path from itertools import product from numpy import array, nan, isnan, float32, empty, sum as nsum import numpy as np @@ -95,18 +94,20 @@ def mst_parallel(ifgs, params): result = empty(shape=(no_ifgs, no_y, no_x), dtype=np.bool) if params[C.PARALLEL]: - log.info('Calculating MST using {} tiles in parallel using {} ' \ - 'processes'.format(no_tiles, ncpus)) + log.info(f'Calculating MST using {no_tiles} tiles in parallel using {ncpus} processes') t_msts = Parallel(n_jobs=params[C.PROCESSES], verbose=joblib_log_level(C.LOG_LEVEL))( delayed(mst_multiprocessing)(t, ifg_paths, params=params) for t in tiles ) for k, tile in enumerate(tiles): - result[:, tile.top_left_y:tile.bottom_right_y, tile.top_left_x: tile.bottom_right_x] = t_msts[k] + y_slice = slice(tile.top_left_y, tile.bottom_right_y) + x_slice = slice(tile.top_left_x, tile.bottom_right_x) + result[:, y_slice, x_slice] = t_msts[k] else: - log.info('Calculating MST using {} tiles in serial'.format(no_tiles)) + log.info(f'Calculating MST using {no_tiles} tiles in serial') for k, tile in enumerate(tiles): - result[:, tile.top_left_y:tile.bottom_right_y, tile.top_left_x: tile.bottom_right_x] = \ - mst_multiprocessing(tile, ifg_paths, params=params) + y_slice = slice(tile.top_left_y, tile.bottom_right_y) + x_slice = slice(tile.top_left_x, tile.bottom_right_x) + result[:, y_slice, x_slice] = mst_multiprocessing(tile, ifg_paths, params=params) return result @@ -155,15 +156,15 @@ def mst_boolean_array(ifgs): """ # The MSTs are stripped of connecting edge info, leaving just the ifgs. nifgs = len(ifgs) - ny, nx = ifgs[0].phase_data.shape - result = empty(shape=(nifgs, ny, nx), dtype=np.bool) + height, width = ifgs[0].phase_data.shape + result = empty(shape=(nifgs, height, width), dtype=np.bool) for y, x, mst in mst_matrix_networkx(ifgs): # mst is a list of datetime.date tuples if isinstance(mst, EdgeView): ifg_sub = [ifg_date_index_lookup(ifgs, d) for d in mst] - ifg_sub_bool = [True if i in ifg_sub else False - for i in range(nifgs)] # boolean conversion + # boolean conversion + ifg_sub_bool = [(i in ifg_sub) for i in range(nifgs)] result[:, y, x] = np.array(ifg_sub_bool) else: result[:, y, x] = np.zeros(nifgs, dtype=bool) @@ -235,7 +236,8 @@ def mst_matrix_networkx(ifgs): if nan_count == 0: yield y, x, edges continue - elif nan_count == nifgs: + + if nan_count == nifgs: yield y, x, nan continue diff --git a/pyrate/core/orbital.py b/pyrate/core/orbital.py index 98d49eddb..050f2cac4 100644 --- a/pyrate/core/orbital.py +++ b/pyrate/core/orbital.py @@ -18,9 +18,8 @@ """ # pylint: disable=invalid-name import tempfile -from typing import Optional, List, Dict, Iterable +from typing import Optional, List from collections import OrderedDict -from pathlib import Path from numpy import empty, isnan, reshape, float32, squeeze from numpy import dot, vstack, zeros, meshgrid import numpy as np @@ -29,7 +28,7 @@ import pyrate.constants as C from pyrate.core.algorithm import first_second_ids, get_all_epochs from pyrate.core import shared, ifgconstants as ifc, prepifg_helper, mst, mpiops -from pyrate.core.shared import nanmedian, Ifg, InputTypes, iterable_split +from pyrate.core.shared import nanmedian, Ifg, iterable_split from pyrate.core.logger import pyratelogger as log from pyrate.prepifg import find_header from pyrate.configuration import MultiplePaths @@ -86,14 +85,14 @@ def remove_orbital_error(ifgs: List, params: dict) -> None: orbfitlksy = params[C.ORBITAL_FIT_LOOKS_Y] # Sanity check of the orbital params - if type(orbfitlksx) != int or type(orbfitlksy) != int: - msg = f"Multi-look factors for orbital correction should be of type: int" + if not isinstance(orbfitlksx, int) or not isinstance(orbfitlksy, int): + msg = "Multi-look factors for orbital correction should be of type: int" raise OrbitalError(msg) if degree not in [PLANAR, QUADRATIC, PART_CUBIC]: - msg = "Invalid degree of %s for orbital correction" % C.ORB_DEGREE_NAMES.get(degree) + msg = f"Invalid degree of {C.ORB_DEGREE_NAMES.get(degree)} for orbital correction" raise OrbitalError(msg) if method not in [NETWORK_METHOD, INDEPENDENT_METHOD]: - msg = "Invalid method of %s for orbital correction" % C.ORB_METHOD_NAMES.get(method) + msg = f"Invalid method of {C.ORB_METHOD_NAMES.get(method)} for orbital correction" raise OrbitalError(msg) # Give informative log messages based on selected options @@ -129,7 +128,9 @@ def remove_orbital_error(ifgs: List, params: dict) -> None: def __create_multilooked_datasets(params): exts, ifg_paths, multi_paths = __extents_from_params(params) - mlooked_datasets = [_create_mlooked_dataset(m, i, exts, params) for m, i in zip(multi_paths, ifg_paths)] + mlooked_datasets = [ + _create_mlooked_dataset(m, i, exts, params) for m, i in zip(multi_paths, ifg_paths) + ] mlooked = [Ifg(m) for m in mlooked_datasets] for m in mlooked: @@ -159,8 +160,9 @@ def _create_mlooked_dataset(multi_path, ifg_path, exts, params): xlooks = params[C.ORBITAL_FIT_LOOKS_X] ylooks = params[C.ORBITAL_FIT_LOOKS_Y] out_path = tempfile.mktemp() - log.debug(f'Multi-looking {ifg_path} with factors X = {xlooks} and Y = {ylooks} for orbital correction') - resampled_data, out_ds = prepifg_helper.prepare_ifg( + log.debug(f'Multi-looking {ifg_path} with factors' \ + 'X = {xlooks} and Y = {ylooks} for orbital correction') + _, out_ds = prepifg_helper.prepare_ifg( ifg_path, xlooks, ylooks, exts, thresh, crop_opt, header, False, out_path ) return out_ds @@ -175,8 +177,8 @@ def _validate_mlooked(mlooked, ifgs): msg = "Mismatching # ifgs and # multilooked ifgs" raise OrbitalError(msg) - if not all([hasattr(i, 'phase_data') for i in mlooked]): - msg = "Mismatching types in multilooked ifgs arg:\n%s" % mlooked + if not all(hasattr(i, 'phase_data') for i in mlooked): + msg = f"Mismatching types in multilooked ifgs arg:\n{mlooked}" raise OrbitalError(msg) @@ -192,8 +194,7 @@ def _get_num_params(degree, intercept: Optional[bool] = False): elif degree == PART_CUBIC: nparams = 6 else: - msg = "Invalid orbital model degree: %s" \ - % C.ORB_DEGREE_NAMES.get(degree) + msg = f"Invalid orbital model degree: {C.ORB_DEGREE_NAMES.get(degree)}" raise OrbitalError(msg) # NB: independent method only, network method handles intercept terms differently @@ -255,7 +256,9 @@ def independent_orbital_correction(ifg_path, params): mlooked_dm = get_design_matrix(ifg, degree, intercept=intercept) # invert to obtain the correction image (forward model) at full-res - orbital_correction = __orb_correction(fullres_dm, mlooked_dm, fullres_phase, vphase, offset=offset) + orbital_correction = __orb_correction( + fullres_dm, mlooked_dm, fullres_phase, vphase, offset=offset + ) # save correction to disc if not orb_on_disc.parent.exists(): @@ -309,7 +312,8 @@ def network_orbital_correction(ifg_paths, params, m_ifgs: Optional[List] = None) :param list ifg_paths: List of Ifg class objects reduced to a minimum spanning tree network :param dict params: dictionary of configuration parameters - :param list m_ifgs: list of multilooked Ifg class objects (sequence must be multilooked versions of 'ifgs' arg) + :param list m_ifgs: The list of multilooked Ifg class objects + (sequence must be multilooked versions of 'ifg_paths' arg) :return: None - interferogram phase data is updated and saved to disk """ @@ -374,11 +378,11 @@ def calc_network_orb_correction(src_ifgs, degree, scale, nepochs, intercept=Fals :param intercept: whether to include a constant offset to fit to each ifg. This intercept is discarded and not returned. - :return coefs: a list of coefficient lists, indexed by epoch. The coefficient lists are in the following order: - - PLANAR - x, y - QUADRATIC - x^2, y^2, x*y, x, y - PART_CUBIC - x*y^2, x^2, y^2, x*y, x, y + :return coefs: a list of coefficient lists, indexed by epoch. + The coefficient lists are in the following order: + PLANAR - x, y + QUADRATIC - x^2, y^2, x*y, x, y + PART_CUBIC - x*y^2, x^2, y^2, x*y, x, y """ vphase = vstack([i.phase_data.reshape((i.num_cells, 1)) for i in src_ifgs]) vphase = squeeze(vphase) @@ -434,8 +438,11 @@ def _save_orbital_error_corrected_phase(ifg, params): orbital fit correction """ # set orbfit tags after orbital error correction - ifg.dataset.SetMetadataItem(ifc.PYRATE_ORB_METHOD, __methods_as_string(params[C.ORBITAL_FIT_METHOD])) - ifg.dataset.SetMetadataItem(ifc.PYRATE_ORB_DEG, __degrees_as_string(params[C.ORBITAL_FIT_DEGREE])) + fit_method = __methods_as_string(params[C.ORBITAL_FIT_METHOD]) + fit_degree = __degrees_as_string(params[C.ORBITAL_FIT_DEGREE]) + + ifg.dataset.SetMetadataItem(ifc.PYRATE_ORB_METHOD, fit_method) + ifg.dataset.SetMetadataItem(ifc.PYRATE_ORB_DEG, fit_degree) ifg.dataset.SetMetadataItem(ifc.PYRATE_ORB_XLOOKS, str(params[C.ORBITAL_FIT_LOOKS_X])) ifg.dataset.SetMetadataItem(ifc.PYRATE_ORB_YLOOKS, str(params[C.ORBITAL_FIT_LOOKS_Y])) ifg.dataset.SetMetadataItem(ifc.PYRATE_ORBITAL_ERROR, ifc.ORB_REMOVED) @@ -540,7 +547,7 @@ def get_network_design_matrix(ifgs, degree, scale, intercept=True): nifgs = len(ifgs) if nifgs < 1: # can feasibly do correction on a single Ifg/2 epochs - raise OrbitalError("Invalid number of Ifgs: %s" % nifgs) + raise OrbitalError(f"Invalid number of Ifgs: {nifgs}") # init sparse network design matrix nepochs = len(set(get_all_epochs(ifgs))) diff --git a/pyrate/core/phase_closure/closure_check.py b/pyrate/core/phase_closure/closure_check.py index 857571b44..496b39c55 100644 --- a/pyrate/core/phase_closure/closure_check.py +++ b/pyrate/core/phase_closure/closure_check.py @@ -15,13 +15,14 @@ # limitations under the License. from collections import defaultdict -from typing import List, Dict, Tuple, Any +from typing import List, Tuple, Any from nptyping import NDArray, UInt16, Float32 import numpy as np import pyrate.constants as C from pyrate.core import mpiops -from pyrate.core.phase_closure.mst_closure import sort_loops_based_on_weights_and_date, WeightedLoop, Edge +from pyrate.core.phase_closure.mst_closure import sort_loops_based_on_weights_and_date, \ + WeightedLoop, Edge from pyrate.configuration import Configuration, MultiplePaths from pyrate.core.phase_closure.sum_closure import sum_phase_closures from pyrate.core.phase_closure.plot_closure import plot_closure @@ -29,9 +30,11 @@ from pyrate.core.logger import pyratelogger as log -def mask_pixels_with_unwrapping_errors(ifgs_breach_count: NDArray[(Any, Any, Any), UInt16], - num_occurrences_each_ifg: NDArray[(Any,), UInt16], - params: dict) -> None: +def mask_pixels_with_unwrapping_errors( + ifgs_breach_count: NDArray[(Any, Any, Any), UInt16], + num_occurrences_each_ifg: NDArray[(Any,), UInt16], + params: dict +): """ Find pixels in the phase data that breach closure_thr, and mask (assign NaNs) to those pixels in those ifgs. @@ -41,7 +44,9 @@ def mask_pixels_with_unwrapping_errors(ifgs_breach_count: NDArray[(Any, Any, Any """ log.debug("Masking phase data of retained ifgs") - for i, m_p in enumerate(params[C.INTERFEROGRAM_FILES]): + ifg_files = params[C.INTERFEROGRAM_FILES] + + for i, m_p in enumerate(ifg_files): pix_index = ifgs_breach_count[:, :, i] == num_occurrences_each_ifg[i] ifg = Ifg(m_p.tmp_sampled_path) ifg.open() @@ -49,14 +54,15 @@ def mask_pixels_with_unwrapping_errors(ifgs_breach_count: NDArray[(Any, Any, Any ifg.phase_data[pix_index] = np.nan ifg.write_modified_phase() - log.info(f"Masked phase data of {i + 1} retained ifgs after phase closure") - return None + log.info(f"Masked phase data of {len(ifg_files)} retained ifgs after phase closure") -def __drop_ifgs_if_not_part_of_any_loop(ifg_files: List[str], loops: List[WeightedLoop], params: dict) -> List[str]: +def __drop_ifgs_if_not_part_of_any_loop( + ifg_files: List[str], loops: List[WeightedLoop], params: dict +) -> List[str]: """ - Check if an ifg is part of any of the loops, otherwise drop it from the list of interferograms for further PyRate - processing. + Check if an ifg is part of any of the loops, otherwise drop it from the list of interferograms + for further PyRate processing. """ loop_ifgs = set() for weighted_loop in loops: @@ -77,7 +83,9 @@ def __drop_ifgs_if_not_part_of_any_loop(ifg_files: List[str], loops: List[Weight return selected_ifg_files -def __drop_ifgs_exceeding_threshold(orig_ifg_files: List[str], ifgs_breach_count, num_occurences_each_ifg, params): +def __drop_ifgs_exceeding_threshold( + orig_ifg_files: List[str], ifgs_breach_count, num_occurences_each_ifg, params +): """ Function to identify and drop ifgs, based on two thresholds. We demand two thresholds to be breached before an ifg is dropped: @@ -85,19 +93,21 @@ def __drop_ifgs_exceeding_threshold(orig_ifg_files: List[str], ifgs_breach_count enough loops to accurately check for unwrapping errors? 2. The second threshold is an average check of pixels breached taking all loops into account. It is evaluated as follows: - (i) ifgs_breach_count contains the number of loops where this pixel in this ifg had a closure exceeding closure_thr. - (b) sum(ifgs_breach_count[:, :, i]) is the number of pixels in ifg exceeding closure_thr over all loops - (c) divide by loop_count_of_this_ifg and num of cells (nrows x ncols) for a normalised measure of threshold. + (i) ifgs_breach_count contains the number of loops where this pixel in this ifg had + a closure exceeding closure_thr. + (b) sum(ifgs_breach_count[:, :, i]) is the number of pixels in ifg exceeding closure_thr + over all loops + (c) divide by loop_count_of_this_ifg and num of cells (nrows x ncols) for a normalised + measure of threshold. """ orig_ifg_files.sort() - nrows, ncols, n_ifgs = ifgs_breach_count.shape + nrows, ncols, _ = ifgs_breach_count.shape selected_ifg_files = [] for i, ifg_file in enumerate(orig_ifg_files): loop_count_of_this_ifg = num_occurences_each_ifg[i] if loop_count_of_this_ifg: # if the ifg participated in at least one loop - ifg_remove_threshold_breached = \ - np.sum(ifgs_breach_count[:, :, i] == loop_count_of_this_ifg) / (nrows * ncols) > params[ - C.IFG_DROP_THR] + ifg_breached = np.sum(ifgs_breach_count[:, :, i] == loop_count_of_this_ifg) + ifg_remove_threshold_breached = ifg_breached / (nrows * ncols) > params[C.IFG_DROP_THR] if not ( # min loops count # check 1 @@ -110,13 +120,15 @@ def __drop_ifgs_exceeding_threshold(orig_ifg_files: List[str], ifgs_breach_count return selected_ifg_files -def iterative_closure_check(config, interactive_plot=True) -> \ - Tuple[List[str], NDArray[(Any, Any, Any), UInt16], NDArray[(Any,), UInt16]]: +def iterative_closure_check( + config, interactive_plot=True +) -> Tuple[List[str], NDArray[(Any, Any, Any), UInt16], NDArray[(Any,), UInt16]]: """ This function iterates the closure check until a stable list of interferogram files is returned. :param config: Configuration class instance :param interactive_plot: bool, whether to plot sum closures of loops - :return: stable list of ifg files, their ifgs_breach_count, and number of occurrences of ifgs in loops + :return: stable list of ifg files, their ifgs_breach_count, and + number of occurrences of ifgs in loops """ params = config.__dict__ ifg_files = [ifg_path.tmp_sampled_path for ifg_path in params[C.INTERFEROGRAM_FILES]] @@ -126,7 +138,8 @@ def iterative_closure_check(config, interactive_plot=True) -> \ log.info(f"Closure check iteration #{i}: working on {len(ifg_files)} ifgs") rets = __wrap_closure_check(config) if rets is None: - return + return None + new_ifg_files, closure, ifgs_breach_count, num_occurences_each_ifg, loops = rets if interactive_plot: if mpiops.rank == 0: @@ -134,19 +147,25 @@ def iterative_closure_check(config, interactive_plot=True) -> \ thr=params[C.CLOSURE_THR], iteration=i) if len(ifg_files) == len(new_ifg_files): break - else: - i += 1 - ifg_files = new_ifg_files # exit condition could be some other check like number_of_loops + + i += 1 + ifg_files = new_ifg_files + # exit condition could be some other check like number_of_loops mpiops.comm.barrier() - log.info(f"Stable list of ifgs achieved after iteration #{i}. {len(ifg_files)} ifgs are retained") + log.info( + f"Stable list of ifgs achieved after iteration #{i}. {len(ifg_files)} ifgs are retained" + ) return ifg_files, ifgs_breach_count, num_occurences_each_ifg -def discard_loops_containing_max_ifg_count(loops: List[WeightedLoop], params) -> List[WeightedLoop]: +def discard_loops_containing_max_ifg_count( + loops: List[WeightedLoop], params +) -> List[WeightedLoop]: """ - This function will discard loops when each ifg participating in a loop has met the max loop count criteria. + This function will discard loops when each ifg participating in a loop has met the + max loop count criteria. :param loops: list of loops :param params: params dict @@ -158,8 +177,8 @@ def discard_loops_containing_max_ifg_count(loops: List[WeightedLoop], params) -> edge_appearances = np.array([ifg_counter[e] for e in loop.edges]) if not np.all(edge_appearances > params[C.MAX_LOOP_REDUNDANCY]): selected_loops.append(loop) - for e in loop.edges: - ifg_counter[e] += 1 + for edge in loop.edges: + ifg_counter[edge] += 1 else: log.debug(f"Loop {loop.loop} ignored: all constituent ifgs have been in a loop " f"{params[C.MAX_LOOP_REDUNDANCY]} times or more") @@ -190,18 +209,25 @@ def __wrap_closure_check(config: Configuration) -> \ if len(sorted_signed_loops) < 1: return None - retained_loops = mpiops.run_once(discard_loops_containing_max_ifg_count, - sorted_signed_loops, params) - ifgs_with_loops = mpiops.run_once(__drop_ifgs_if_not_part_of_any_loop, ifg_files, retained_loops, params) + retained_loops = mpiops.run_once( + discard_loops_containing_max_ifg_count, + sorted_signed_loops, params + ) + + ifgs_with_loops = mpiops.run_once( + __drop_ifgs_if_not_part_of_any_loop, + ifg_files, retained_loops, params + ) msg = f"After applying MAX_LOOP_REDUNDANCY = {params[C.MAX_LOOP_REDUNDANCY]} criteria, " \ f"{len(retained_loops)} loops are retained" if len(retained_loops) < 1: return None - else: - log.info(msg) - closure, ifgs_breach_count, num_occurences_each_ifg = sum_phase_closures(ifgs_with_loops, retained_loops, params) + log.info(msg) + + closure_sum = sum_phase_closures(ifgs_with_loops, retained_loops, params) + closure, ifgs_breach_count, num_occurences_each_ifg = closure_sum if mpiops.rank == 0: closure_ins = config.closure() @@ -210,8 +236,10 @@ def __wrap_closure_check(config: Configuration) -> \ np.save(closure_ins.num_occurences_each_ifg, num_occurences_each_ifg) np.save(closure_ins.loops, retained_loops, allow_pickle=True) - selected_ifg_files = mpiops.run_once(__drop_ifgs_exceeding_threshold, - ifgs_with_loops, ifgs_breach_count, num_occurences_each_ifg, params) + selected_ifg_files = mpiops.run_once( + __drop_ifgs_exceeding_threshold, + ifgs_with_loops, ifgs_breach_count, num_occurences_each_ifg, params + ) # update the ifg list in the parameters dictionary params[C.INTERFEROGRAM_FILES] = \ @@ -219,7 +247,6 @@ def __wrap_closure_check(config: Configuration) -> \ return selected_ifg_files, closure, ifgs_breach_count, num_occurences_each_ifg, retained_loops -# TODO: Consider whether this helper function is better homed in a generic module and used more widely def update_ifg_list(ifg_files: List[str], multi_paths: List[MultiplePaths]) -> List[MultiplePaths]: """ Function to extract full paths for a subsetted list of interferograms diff --git a/pyrate/core/phase_closure/collect_loops.py b/pyrate/core/phase_closure/collect_loops.py index 0ab015383..a709bfb5f 100644 --- a/pyrate/core/phase_closure/collect_loops.py +++ b/pyrate/core/phase_closure/collect_loops.py @@ -24,7 +24,7 @@ def dfs(graph, marked, n, vert, start, count, loop, all_loops): Adapted from https://www.geeksforgeeks.org/print-all-the-cycles-in-an-undirected-graph/ """ # Number of vertices - V = graph.shape[0] + V = graph.shape[0] # mark the vertex vert as visited marked[vert] = True @@ -39,9 +39,8 @@ def dfs(graph, marked, n, vert, start, count, loop, all_loops): if graph[vert][start] == 1: count += 1 all_loops.append(loop) - return count - else: - return count + + return count # For searching every possible path of length (n-1) for i in range(V): @@ -111,11 +110,11 @@ def dedupe_loops(loops: List[List]) -> List: """ seen_sets = set() filtered = [] - for l in loops: - loop = l[:] - l.sort() - l = tuple(l) - if l not in seen_sets: - seen_sets.add(l) - filtered.append(loop) + for loop in loops: + orig_loop = loop[:] + loop.sort() + loop = tuple(loop) + if loop not in seen_sets: + seen_sets.add(loop) + filtered.append(orig_loop) return filtered diff --git a/pyrate/core/phase_closure/mst_closure.py b/pyrate/core/phase_closure/mst_closure.py index cc59a8210..dbbb34d8a 100644 --- a/pyrate/core/phase_closure/mst_closure.py +++ b/pyrate/core/phase_closure/mst_closure.py @@ -16,7 +16,7 @@ from collections import namedtuple -from typing import List, Union +from typing import List from datetime import date import numpy as np import networkx as nx @@ -29,6 +29,7 @@ class SignedEdge: + """Represents an edge (line between two points) with an ssociated sign value""" def __init__(self, edge: Edge, sign: int): self.edge = edge @@ -41,6 +42,7 @@ def __repr__(self): class SignedWeightedEdge(SignedEdge): + """Represents an signed edge value (`SignedEdge`) with an associated weight value""" def __init__(self, signed_edge: SignedEdge, weight: int): super().__init__(signed_edge.edge, sign=signed_edge.sign) @@ -50,7 +52,8 @@ def __init__(self, signed_edge: SignedEdge, weight: int): class WeightedLoop: """ - Loop with weight equal to the sum of the edge weights, where edge weights are the ifg duration in days + Loop with weight equal to the sum of the edge weights, + where edge weights are the ifg duration in days """ def __init__(self, loop: List[SignedWeightedEdge]): @@ -58,31 +61,34 @@ def __init__(self, loop: List[SignedWeightedEdge]): @property def weight(self): + """Returns the weight value of the whole loop (sum of loop edge weights)""" return sum([swe.weight for swe in self.loop]) @property def earliest_date(self): + """Returns the earliest date in the loop""" return min({swe.first for swe in self.loop}) @property def primary_dates(self): + """Returns the primary dates for each edge in the loop as a string""" first_dates = [swe.first for swe in self.loop] first_dates.sort() - st = ''.join([str(d) for d in first_dates]) - return st + return ''.join([str(d) for d in first_dates]) @property def secondary_dates(self): + """Returns the secondary dates for each edge in the loop as a string""" first_dates = [swe.second for swe in self.loop] first_dates.sort() - st = ''.join([str(d) for d in first_dates]) - return st + return ''.join([str(d) for d in first_dates]) def __len__(self): return len(self.loop) @property def edges(self): + """Returns this loops edge values""" return [Edge(swe.first, swe.edge.second) for swe in self.loop] @@ -103,36 +109,36 @@ def __find_closed_loops(edges: List[Edge], max_loop_length: int) -> List[List[da loops.extend(loops_) node_list = g.nodes() - node_list_dict = {i: n for i, n in enumerate(node_list)} + node_list_dict = dict(enumerate(node_list)) loop_subset = [] - for l in loops: - loop = [] - for ll in l: - loop.append(node_list_dict[ll]) - loop_subset.append(loop) + for loop in loops: + new_loop = [node_list_dict[i] for i in loop] + loop_subset.append(new_loop) log.debug(f"Total number of loops is {len(loop_subset)}") return loop_subset -def __add_signs_and_weights_to_loops(loops: List[List[date]], available_edges: List[Edge]) -> List[WeightedLoop]: +def __add_signs_and_weights_to_loops( + loops: List[List[date]], available_edges: List[Edge] +) -> List[WeightedLoop]: """ add signs and weights to loops. Additionally, sort the loops (change order of ifgs appearing in loop) by weight and date """ weighted_signed_loops = [] available_edges = set(available_edges) # hash it once for O(1) lookup - for i, l in enumerate(loops): + for loop in loops: weighted_signed_loop = [] - l.append(l[0]) # add the closure loop - for ii, ll in enumerate(l[:-1]): - if l[ii+1] > ll: - edge = Edge(ll, l[ii+1]) + loop.append(loop[0]) # add the closure loop + for i, loop_idx in enumerate(loop[:-1]): + if loop[i+1] > loop_idx: + edge = Edge(loop_idx, loop[i+1]) assert edge in available_edges signed_edge = SignedEdge(edge, 1) # opposite direction of ifg else: - edge = Edge(l[ii+1], ll) + edge = Edge(loop[i+1], loop_idx) assert edge in available_edges signed_edge = SignedEdge(edge, -1) # in direction of ifg weighted_signed_edge = SignedWeightedEdge( @@ -163,7 +169,8 @@ def __find_signed_closed_loops(params: dict) -> List[WeightedLoop]: ifg_files.sort() log.debug(f"The number of ifgs in the list is {len(ifg_files)}") available_edges = __setup_edges(ifg_files) - all_loops = __find_closed_loops(available_edges, max_loop_length=params[C.MAX_LOOP_LENGTH]) # find loops with weights + # find loops with weights + all_loops = __find_closed_loops(available_edges, max_loop_length=params[C.MAX_LOOP_LENGTH]) signed_weighted_loops = __add_signs_and_weights_to_loops(all_loops, available_edges) return signed_weighted_loops diff --git a/pyrate/core/phase_closure/plot_closure.py b/pyrate/core/phase_closure/plot_closure.py index 8815a5be6..39a611985 100644 --- a/pyrate/core/phase_closure/plot_closure.py +++ b/pyrate/core/phase_closure/plot_closure.py @@ -16,32 +16,42 @@ from typing import List -import numpy as np from pathlib import Path +import numpy as np from pyrate.core.phase_closure.mst_closure import WeightedLoop from pyrate.core.logger import pyratelogger as log from pyrate.configuration import Configuration -def plot_closure(closure: np.ndarray, loops: List[WeightedLoop], - config: Configuration, thr: float, iteration: int): +def plot_closure( + closure: np.ndarray, + loops: List[WeightedLoop], + config: Configuration, + thr: float, + iteration: int +): + """ + Produces a set of graphial plots for the closure loops. + + The plots will be saved into the `config.phase_closure_dir` directory. + """ thr = thr * np.pi try: import matplotlib.pyplot as plt import matplotlib as mpl from mpl_toolkits.axes_grid1 import make_axes_locatable cmap = mpl.cm.Spectral - except ImportError as e: - log.warn(ImportError(e)) + except ImportError as error: + log.warn(ImportError(error)) log.warn("Required plotting packages are not found in environment. " "Closure loop plot will not be generated!!!") return - nrows, ncols, n_loops = closure.shape + _, _, n_loops = closure.shape # 49 ifgs per fig - plt_rows = 7 + plt_rows = 7 plt_cols = 7 plots_per_fig = plt_rows * plt_cols n_figs = n_loops // plots_per_fig + (n_loops % plots_per_fig > 0) @@ -52,14 +62,14 @@ def plot_closure(closure: np.ndarray, loops: List[WeightedLoop], this_fig_plots = 0 for p_r in range(plt_rows): for p_c in range(plt_cols): - if all_fig_plots == n_loops + 1: + if all_fig_plots == n_loops + 1: break - + ax = fig.add_subplot(plt_rows, plt_cols, plt_cols * p_r + p_c + 1) data = closure[:, :, plt_cols * p_r + p_c + fig_i * plots_per_fig] loop = loops[plt_cols * p_r + p_c + fig_i * plots_per_fig] title = ',\n'.join([repr(l) for l in loop.loop]) - im = ax.imshow(data, vmin=-thr, vmax=thr, cmap=cmap) + image = ax.imshow(data, vmin=-thr, vmax=thr, cmap=cmap) plt.tick_params(axis='both', which='major', labelsize=12) text = ax.set_title(title) text.set_fontsize(16) @@ -67,13 +77,14 @@ def plot_closure(closure: np.ndarray, loops: List[WeightedLoop], divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(im, cax=cax) + plt.colorbar(image, cax=cax) this_fig_plots += 1 all_fig_plots += 1 fig.tight_layout() - closure_plot_file = Path(config.phase_closure_dir).joinpath(f'closure_loops_iteration_{iteration}_fig_{fig_i}.png') + closure_plot_file = f'closure_loops_iteration_{iteration}_fig_{fig_i}.png' + closure_plot_file = Path(config.phase_closure_dir) / closure_plot_file plt.savefig(closure_plot_file, dpi=100) plt.close(fig) log.info(f'{this_fig_plots} closure loops plotted in {closure_plot_file}') diff --git a/pyrate/core/phase_closure/sum_closure.py b/pyrate/core/phase_closure/sum_closure.py index 819adbd39..d7ff12c15 100644 --- a/pyrate/core/phase_closure/sum_closure.py +++ b/pyrate/core/phase_closure/sum_closure.py @@ -14,7 +14,6 @@ # See the License for the specific language governing permissions and # limitations under the License. -import resource from collections import namedtuple from typing import List, Dict, Tuple, Any from nptyping import NDArray, Float32, UInt16 @@ -61,11 +60,12 @@ def _func(ifg, index): return ret_combined -def sum_phase_closures(ifg_files: List[str], loops: List[WeightedLoop], params: dict) -> \ - Tuple[NDArray[(Any, Any, Any), Float32], NDArray[(Any, Any, Any), UInt16], NDArray[(Any,), UInt16]]: +def sum_phase_closures( + ifg_files: List[str], loops: List[WeightedLoop], params: dict +) -> Tuple[NDArray[(Any, Any, Any), Float32], NDArray[(Any, Any, Any), UInt16], NDArray[(Any,), UInt16]]: """ Compute the closure sum for each pixel in each loop, and count the number of times a pixel - contributes to a failed closure loop (where the summed closure is above/below the + contributes to a failed closure loop (where the summed closure is above/below the CLOSURE_THR threshold). :param ifg_files: list of ifg files :param loops: list of loops @@ -80,24 +80,29 @@ def sum_phase_closures(ifg_files: List[str], loops: List[WeightedLoop], params: ifgs = [v.IfgPhase for v in edge_to_indexed_ifgs.values()] n_ifgs = len(ifgs) + ifg_shape = ifgs[0].phase_data.shape + if params[C.PARALLEL]: # rets = Parallel(n_jobs=params[cf.PROCESSES], verbose=joblib_log_level(cf.LOG_LEVEL))( - # delayed(__compute_ifgs_breach_count)(ifg0, n_ifgs, weighted_loop, edge_to_indexed_ifgs, params) + # delayed(__compute_ifgs_breach_count)( + # ifg0, n_ifgs, weighted_loop, edge_to_indexed_ifgs, params + # ) # for weighted_loop in loops # ) # for k, r in enumerate(rets): # closure_dict[k], ifgs_breach_count_dict[k] = r # TODO: enable multiprocessing - needs pickle error workaround - closure = np.zeros(shape=(* ifgs[0].phase_data.shape, len(loops)), dtype=np.float32) - ifgs_breach_count = np.zeros(shape=(ifgs[0].phase_data.shape + (n_ifgs,)), dtype=np.uint16) + closure = np.zeros(shape=(*ifg_shape, len(loops)), dtype=np.float32) + ifgs_breach_count = np.zeros(shape=(ifg_shape + (n_ifgs,)), dtype=np.uint16) for k, weighted_loop in enumerate(loops): - closure[:, :, k], ifgs_breach_count_l = __compute_ifgs_breach_count(weighted_loop, edge_to_indexed_ifgs, - params) + closure[:, :, k], ifgs_breach_count_l = __compute_ifgs_breach_count( + weighted_loop, edge_to_indexed_ifgs, params + ) ifgs_breach_count += ifgs_breach_count_l else: process_loops = mpiops.array_split(loops) - closure_process = np.zeros(shape=(* ifgs[0].phase_data.shape, len(process_loops)), dtype=np.float32) - ifgs_breach_count_process = np.zeros(shape=(ifgs[0].phase_data.shape + (n_ifgs,)), dtype=np.uint16) + closure_process = np.zeros(shape=(* ifg_shape, len(process_loops)), dtype=np.float32) + ifgs_breach_count_process = np.zeros(shape=(ifg_shape + (n_ifgs,)), dtype=np.uint16) for k, weighted_loop in enumerate(process_loops): closure_process[:, :, k], ifgs_breach_count_l = \ __compute_ifgs_breach_count(weighted_loop, edge_to_indexed_ifgs, params) @@ -111,40 +116,50 @@ def sum_phase_closures(ifg_files: List[str], loops: List[WeightedLoop], params: total_gb = mpiops.comm.allreduce(closure_process.nbytes / 1e9, op=mpiops.MPI.SUM) log.debug(f"Memory usage to compute closure_process was {total_gb} GB") if mpiops.rank == 0: - ifgs_breach_count = np.zeros(shape=(ifgs[0].phase_data.shape + (n_ifgs,)), dtype=np.uint16) + ifgs_breach_count = np.zeros(shape=(ifg_shape + (n_ifgs,)), dtype=np.uint16) # closure - closure = np.zeros(shape=(* ifgs[0].phase_data.shape, len(loops)), dtype=np.float32) + closure = np.zeros(shape=(* ifg_shape, len(loops)), dtype=np.float32) main_process_indices = mpiops.array_split(range(len(loops))).astype(np.uint16) closure[:, :, main_process_indices] = closure_process for rank in range(1, mpiops.size): rank_indices = mpiops.array_split(range(len(loops)), rank).astype(np.uint16) - this_rank_closure = np.zeros(shape=(* ifgs[0].phase_data.shape, len(rank_indices)), dtype=np.float32) - mpiops.comm.Recv(this_rank_closure, source=rank, tag=rank) - closure[:, :, rank_indices] = this_rank_closure + rank_closure = np.zeros(shape=(*ifg_shape, len(rank_indices)), dtype=np.float32) + mpiops.comm.Recv(rank_closure, source=rank, tag=rank) + closure[:, :, rank_indices] = rank_closure else: closure = None ifgs_breach_count = None mpiops.comm.Send(closure_process, dest=0, tag=mpiops.rank) if mpiops.MPI_INSTALLED: - mpiops.comm.Reduce([ifgs_breach_count_process, mpiops.MPI.UINT16_T], - [ifgs_breach_count, mpiops.MPI.UINT16_T], op=mpiops.MPI.SUM, root=0) # global + mpiops.comm.Reduce( + [ifgs_breach_count_process, mpiops.MPI.UINT16_T], + [ifgs_breach_count, mpiops.MPI.UINT16_T], + op=mpiops.MPI.SUM, + root=0 + ) # global else: - ifgs_breach_count = mpiops.comm.reduce(ifgs_breach_count_process, op=mpiops.sum0_op, root=0) + ifgs_breach_count = mpiops.comm.reduce( + ifgs_breach_count_process, + op=mpiops.sum0_op, + root=0 + ) - log.debug(f"successfully summed phase closure breach array") + log.debug("successfully summed phase closure breach array") - num_occurrences_each_ifg = None + num_occurrences = None if mpiops.rank == 0: - num_occurrences_each_ifg = _find_num_occurrences_each_ifg(loops, edge_to_indexed_ifgs, n_ifgs) + num_occurrences = _find_num_occurrences_each_ifg(loops, edge_to_indexed_ifgs, n_ifgs) - return closure, ifgs_breach_count, num_occurrences_each_ifg + return closure, ifgs_breach_count, num_occurrences -def _find_num_occurrences_each_ifg(loops: List[WeightedLoop], - edge_to_indexed_ifgs: Dict[Edge, IndexedIfg], - n_ifgs: int) -> NDArray[(Any,), UInt16]: +def _find_num_occurrences_each_ifg( + loops: List[WeightedLoop], + edge_to_indexed_ifgs: Dict[Edge, IndexedIfg], + n_ifgs: int +) -> NDArray[(Any,), UInt16]: """find how many times each ifg appears in total in all loops""" num_occurrences_each_ifg = np.zeros(shape=n_ifgs, dtype=np.uint16) for weighted_loop in loops: @@ -155,9 +170,11 @@ def _find_num_occurrences_each_ifg(loops: List[WeightedLoop], return num_occurrences_each_ifg -def __compute_ifgs_breach_count(weighted_loop: WeightedLoop, - edge_to_indexed_ifgs: Dict[Edge, IndexedIfg], params: dict) \ - -> Tuple[NDArray[(Any, Any), Float32], NDArray[(Any, Any, Any), UInt16]]: +def __compute_ifgs_breach_count( + weighted_loop: WeightedLoop, + edge_to_indexed_ifgs: Dict[Edge, IndexedIfg], + params: dict +) -> Tuple[NDArray[(Any, Any), Float32], NDArray[(Any, Any, Any), UInt16]]: """Compute summed `closure` of each loop, and compute `ifgs_breach_count` for each pixel.""" n_ifgs = len(edge_to_indexed_ifgs) indexed_ifg = list(edge_to_indexed_ifgs.values())[0] @@ -176,14 +193,16 @@ def __compute_ifgs_breach_count(weighted_loop: WeightedLoop, if use_median: closure -= np.nanmedian(closure) # optionally subtract the median closure phase - # this will deal with nans in `closure`, i.e., nans are not selected in indices_breaching_threshold + # this will deal with nans in `closure`, + # i.e., nans are not selected in indices_breaching_threshold indices_breaching_threshold = np.absolute(closure) > closure_thr for signed_edge in weighted_loop.loop: ifg_index = edge_to_indexed_ifgs[signed_edge.edge].index # the variable 'ifgs_breach_count' is increased by 1 for that pixel # make sure we are not incrementing the nan positions in the closure - # as we don't know the phase of these pixels and also they were converted to zero before closure check - # Therefore, we leave them out of ifgs_breach_count, i.e., we don't increment their ifgs_breach_count values + # as we don't know the phase of these pixels and also they were converted to + # zero before closure check. Therefore, we leave them out of ifgs_breach_count, + # i.e., we don't increment their ifgs_breach_count values ifgs_breach_count[indices_breaching_threshold, ifg_index] += 1 return closure, ifgs_breach_count diff --git a/pyrate/core/prepifg_helper.py b/pyrate/core/prepifg_helper.py index 34cbe9a37..a780711e7 100644 --- a/pyrate/core/prepifg_helper.py +++ b/pyrate/core/prepifg_helper.py @@ -30,7 +30,8 @@ from typing import List, Tuple, Union, Callable from numpy import array, nan, isnan, nanmean, float32, zeros, sum as nsum -from pyrate.constants import sixteen_digits_pattern, COHERENCE_FILE_PATHS, IFG_CROP_OPT, IFG_LKSX, IFG_LKSY +from pyrate.constants import sixteen_digits_pattern, COHERENCE_FILE_PATHS, IFG_CROP_OPT, \ + IFG_LKSX, IFG_LKSY from pyrate.configuration import ConfigException from pyrate.core.gdal_python import crop_resample_average from pyrate.core.shared import dem_or_ifg, Ifg, DEM @@ -48,8 +49,13 @@ GRID_TOL = 1e-6 -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]: +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. @@ -64,15 +70,15 @@ def get_analysis_extent(crop_opt: int, rasters: List[Union[Ifg, DEM]], xlooks: i """ if crop_opt not in CROP_OPTIONS: - raise PreprocessError("Unrecognised crop option: %s" % crop_opt) + raise PreprocessError(f"Unrecognised crop option: {crop_opt}") if crop_opt == CUSTOM_CROP: if not user_exts: raise PreprocessError('No custom cropping extents specified') - elif len(user_exts) != 4: # check for required numbers + if len(user_exts) != 4: # check for required numbers raise PreprocessError('Custom extents must have all 4 values') - elif len(user_exts) == 4: # check for non floats - if not all([_is_number(z) for z in user_exts]): + if len(user_exts) == 4: # check for non floats + if not all(_is_number(z) for z in user_exts): raise PreprocessError('Custom extents must be 4 numbers') _check_looks(xlooks, ylooks) @@ -101,12 +107,12 @@ def _check_looks(xlooks, ylooks): if not (isinstance(xlooks, Number) and isinstance(ylooks, Number)): # pragma: no cover - msg = "Non-numeric looks parameter(s), x: %s, y: %s" % (xlooks, ylooks) + msg = f"Non-numeric looks parameter(s), x: {xlooks}, y: {ylooks}" raise PreprocessError(msg) if not (xlooks > 0 and ylooks > 0): # pragma: no cover - msg = "Invalid looks parameter(s), x: %s, y: %s. " \ - "Looks must be an integer greater than zero" % (xlooks, ylooks) + msg = f"Invalid looks parameter(s), x: {xlooks}, y: {ylooks}. " \ + "Looks must be an integer greater than zero" raise PreprocessError(msg) if xlooks != ylooks: @@ -126,7 +132,7 @@ def _check_resolution(ifgs: List[Union[Ifg, DEM]]): i.close() values = array(values) if not (values == values[0]).all(): # pragma: no cover - msg = "Grid resolution does not match for %s" % var + msg = f"Grid resolution does not match for {var}" raise PreprocessError(msg) @@ -152,8 +158,10 @@ def _get_extents(ifgs, crop_opt, user_exts=None): return extents -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): +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. Returns are only given if write_to_disk=False @@ -195,8 +203,15 @@ def prepare_ifg(raster_path, xlooks, ylooks, exts, thresh, crop_opt, header, wri # # reproject() driver_type = 'GTiff' if write_to_disk else 'MEM' resampled_data, out_ds = crop_resample_average( - input_tif=raster.data_path, extents=exts, new_res=resolution, output_file=out_path, thresh=thresh, - out_driver_type=driver_type, hdr=header, coherence_path=coherence_path, coherence_thresh=coherence_thresh + input_tif=raster.data_path, + extents=exts, + new_res=resolution, + output_file=out_path, + thresh=thresh, + out_driver_type=driver_type, + hdr=header, + coherence_path=coherence_path, + coherence_thresh=coherence_thresh ) return resampled_data, out_ds @@ -422,9 +437,13 @@ def coherence_paths_for(path: str, params: dict, tif=False) -> str: _, filename = split(path) epoch = re.search(sixteen_digits_pattern, filename).group(0) if tif: - coh_file_paths = [f.converted_path for f in params[COHERENCE_FILE_PATHS] if epoch in f.converted_path] + coh_file_paths = [ + f.converted_path for f in params[COHERENCE_FILE_PATHS] if epoch in f.converted_path + ] else: - coh_file_paths = [f.unwrapped_path for f in params[COHERENCE_FILE_PATHS] if epoch in f.unwrapped_path] + coh_file_paths = [ + f.unwrapped_path for f in params[COHERENCE_FILE_PATHS] if epoch in f.unwrapped_path + ] if len(coh_file_paths) > 1: raise ConfigException(f"found more than one coherence " diff --git a/pyrate/core/ref_phs_est.py b/pyrate/core/ref_phs_est.py index 50ddcf134..73980fe25 100644 --- a/pyrate/core/ref_phs_est.py +++ b/pyrate/core/ref_phs_est.py @@ -17,8 +17,6 @@ """ This Python module implements a reference phase estimation algorithm. """ -from pathlib import Path -from typing import List from joblib import Parallel, delayed import numpy as np @@ -69,13 +67,15 @@ def _inner(ifg_paths): else: ref_phs = np.zeros(len(ifgs)) for n, ifg in enumerate(ifgs): - ref_phs[n] = _est_ref_phs_patch_median(phase_data[n], half_chip_size, refpx, refpy, thresh) + ref_phs[n] = _est_ref_phs_patch_median( + phase_data[n], half_chip_size, refpx, refpy, thresh + ) return ref_phs - + process_ifgs_paths = mpiops.array_split(ifg_paths) ref_phs = _inner(process_ifgs_paths) - return ref_phs + return ref_phs def _est_ref_phs_patch_median(phase_data, half_chip_size, refpx, refpy, thresh): @@ -87,10 +87,12 @@ def _est_ref_phs_patch_median(phase_data, half_chip_size, refpx, refpy, thresh): patch = np.reshape(patch, newshape=(-1, 1), order='F') nanfrac = np.sum(~np.isnan(patch)) if nanfrac < thresh: - raise ReferencePhaseError('The data window at the reference pixel ' - 'does not have enough valid observations. ' - 'Actual = {}, Threshold = {}.'.format( - nanfrac, thresh)) + msg = 'The data window at the reference pixel ' \ + 'does not have enough valid observations. ' \ + f'Actual = {nanfrac}, Threshold = {thresh}.' + + raise ReferencePhaseError(msg) + ref_ph = nanmedian(patch) return ref_ph @@ -125,9 +127,7 @@ def _process_phase_sum(ifg_paths): return ifg_phase_data_sum def _inner(proc_ifgs, phase_data_sum): - if isinstance(proc_ifgs[0], Ifg): - proc_ifgs = proc_ifgs - else: + if not isinstance(proc_ifgs[0], Ifg): proc_ifgs = [Ifg(ifg_path) for ifg_path in proc_ifgs] for ifg in proc_ifgs: @@ -184,15 +184,14 @@ def __inner(ifg, ref_ph): ifg.close() log.info("Correcting ifgs by subtracting reference phase") - for i, rp in zip(mpiops.array_split(ifgs), mpiops.array_split(ref_phs)): - __inner(i, rp) + for ifg, ref_phase in zip(mpiops.array_split(ifgs), mpiops.array_split(ref_phs)): + __inner(ifg, ref_phase) class ReferencePhaseError(Exception): """ Generic class for errors in reference phase estimation. """ - pass def ref_phase_est_wrapper(params): @@ -203,14 +202,14 @@ def ref_phase_est_wrapper(params): refpx, refpy = params[C.REFX_FOUND], params[C.REFY_FOUND] if len(ifg_paths) < 2: raise ReferencePhaseError( - "At least two interferograms required for reference phase correction ({len_ifg_paths} " - "provided).".format(len_ifg_paths=len(ifg_paths)) + f"At least two interferograms required for reference" + f" phase correction ({len(ifg_paths)} provided)." ) # this is not going to be true as we now start with fresh multilooked ifg copies - remove? if mpiops.run_once(shared.check_correction_status, ifg_paths, ifc.PYRATE_REF_PHASE): log.warning('Reference phase correction already applied to ifgs; returning') - return + return None ifgs = [Ifg(ifg_path) for ifg_path in ifg_paths] # Save reference phase numpy arrays to disk. @@ -228,7 +227,9 @@ def ref_phase_est_wrapper(params): log.info("Calculating reference phase as median of interferogram") ref_phs = est_ref_phase_ifg_median(ifg_paths, params) elif params[C.REF_EST_METHOD] == 2: - log.info('Calculating reference phase in a patch surrounding pixel (x, y): ({}, {})'.format(refpx, refpy)) + log.info( + f'Calculating reference phase in a patch surrounding pixel (x, y): ({refpx}, {refpy})' + ) ref_phs = est_ref_phase_patch_median(ifg_paths, params, refpx, refpy) else: raise ReferencePhaseError("No such option, set parameter 'refest' to '1' or '2'.") diff --git a/pyrate/core/refpixel.py b/pyrate/core/refpixel.py index 3ad6fe23b..db55ded1e 100644 --- a/pyrate/core/refpixel.py +++ b/pyrate/core/refpixel.py @@ -43,13 +43,16 @@ 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) + pyrate_refpix_lon, pyrate_refpix_lat = mpiops.run_once( + convert_pixel_coord_to_lat_lon, + 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) + log.debug(f"Updating metadata for: {ifg_file}") ifg = Ifg(ifg_file) log.debug("Open dataset") ifg.open(readonly=True) @@ -75,40 +78,58 @@ def update_refpix_metadata(ifg_paths, refx, refy, transform, params): ifg.close() -def convert_pixel_value_to_geographic_coordinate(refx, refy, transform): +def convert_pixel_coord_to_lat_lon(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. + refy: The pixel y coordinate. transform: The geotransform array of the image. Returns: Tuple of lon, lat geographic coordinate. """ - lon = lon_from_pixel_coordinate(refx, transform) + lon = lon_from_pixel_coord(refx, transform) - lat = lat_from_pixel_coordinate(refy, transform) + lat = lat_from_pixel_coord(refy, transform) return lon, lat -def lat_from_pixel_coordinate(refy, transform): - yOrigin = transform[3] - pixelHeight = -transform[5] - lat = yOrigin - refy * pixelHeight +def lat_from_pixel_coord(refy, transform): + """ + Gets the pixel Y coordinate for a latitude coordinate, given a geotransform for an image. + + Args: + refy: The pixel y coordinate. + transform: The geotransform array of the image. + Returns: + The resulting latitude coordinate. + """ + y_origin = transform[3] + pixel_height = -transform[5] + lat = y_origin - refy * pixel_height return lat -def lon_from_pixel_coordinate(refx, transform): - xOrigin = transform[0] - pixelWidth = transform[1] - lon = refx * pixelWidth + xOrigin +def lon_from_pixel_coord(refx, transform): + """ + Gets the pixel X coordinate for a longitude coordinate, given a geotransform for an image. + + Args: + refy: The pixel x coordinate. + transform: The geotransform array of the image. + Returns: + The resulting longitude coordinate. + """ + x_origin = transform[0] + pixel_width = transform[1] + lon = refx * pixel_width + x_origin return lon -def convert_geographic_coordinate_to_pixel_value(lon, lat, transform): +def convert_lat_lon_to_pixel_coord(lon, lat, transform): """ Converts a latitude/longitude coordinate to a pixel coordinate given the geotransform of the image. @@ -120,13 +141,13 @@ def convert_geographic_coordinate_to_pixel_value(lon, lat, transform): Tuple of refx, refy pixel coordinates. """ - xOrigin = transform[0] - yOrigin = transform[3] - pixelWidth = transform[1] - pixelHeight = -transform[5] + x_origin = transform[0] + y_origin = transform[3] + pixel_width = transform[1] + pixel_height = -transform[5] - refx = round((lon - xOrigin) / pixelWidth) - refy = round((yOrigin - lat) / pixelHeight) + refx = round((lon - x_origin) / pixel_width) + refy = round((y_origin - lat) / pixel_height) return int(refx), int(refy) @@ -186,7 +207,7 @@ def find_min_mean(mean_sds, grid): :param list grid: List of ref pixel coordinates tuples :return: Tuple of (refy, refx) with minimum mean - :rtype: tuple + :rtype: tuple """ log.debug('Ranking ref pixel candidates based on mean values') try: @@ -201,10 +222,10 @@ def ref_pixel_setup(ifgs_or_paths, params): """ Sets up the grid for reference pixel computation and saves numpy files to disk for later use during ref pixel computation. - + :param list ifgs_or_paths: List of interferogram filenames or Ifg objects :param dict params: Dictionary of configuration parameters - + :return: half_patch_size: size of patch :rtype: float :return: thresh @@ -263,12 +284,14 @@ def save_ref_pixel_blocks(grid, half_patch_size, ifg_paths, params): ifg.nodata_value = params[C.NO_DATA_VALUE] ifg.convert_to_nans() ifg.convert_to_mm() + + b = os.path.basename(pth).split('.')[0] + for y, x in grid: data = ifg.phase_data[y - half_patch_size:y + half_patch_size + 1, x - half_patch_size:x + half_patch_size + 1] - data_file = join(outdir, 'ref_phase_data_{b}_{y}_{x}.npy'.format( - b=os.path.basename(pth).split('.')[0], y=y, x=x)) + data_file = join(outdir, f'ref_phase_data_{b}_{y}_{x}.npy') np.save(file=data_file, arr=data) ifg.close() log.debug('Saved ref pixel blocks') @@ -299,10 +322,8 @@ def _ref_pixel_multi(g, half_patch_size, phase_data_or_ifg_paths, data = [] output_dir = params[C.TMPDIR] for p in phase_data_or_ifg_paths: - data_file = os.path.join(output_dir, - 'ref_phase_data_{b}_{y}_{x}.npy'.format( - b=os.path.basename(p).split('.')[0], - y=y, x=x)) + b = os.path.basename(p).split('.')[0] + data_file = os.path.join(output_dir, f'ref_phase_data_{b}_{y}_{x}.npy') data.append(np.load(file=data_file)) else: # phase_data_or_ifg is phase_data list data = [p[y - half_patch_size:y + half_patch_size + 1, @@ -312,8 +333,8 @@ def _ref_pixel_multi(g, half_patch_size, phase_data_or_ifg_paths, if all(valid): # ignore if 1+ ifgs have too many incoherent cells sd = [std(i[~isnan(i)]) for i in data] return mean(sd) - else: - return np.nan + + return np.nan def _step(dim, ref, radius): @@ -432,19 +453,21 @@ def ref_pixel_calc_wrapper(params: dict) -> Tuple[int, int]: ref_pixel_file = Configuration.ref_pixel_path(params) - def __reuse_ref_pixel_file_if_exists(): - if ref_pixel_file.exists(): - refx, refy = np.load(ref_pixel_file) - log.info('Reusing pre-calculated ref-pixel values: ({}, {}) from file {}'.format( - refx, refy, ref_pixel_file.as_posix())) - log.warning("Reusing ref-pixel values from previous run!!!") - params[C.REFX_FOUND], params[C.REFY_FOUND] = int(refx), int(refy) - return int(refx), int(refy) - else: + def _reuse_ref_pixel_file_if_exists(): + if not ref_pixel_file.exists(): return None, None + refx, refy = np.load(ref_pixel_file) + ref_path = ref_pixel_file.as_posix() + log.info( + f'Reusing pre-calculated ref-pixel values: ({refx}, {refy}) from file {ref_path}' + ) + log.warning("Reusing ref-pixel values from previous run!!!") + params[C.REFX_FOUND], params[C.REFY_FOUND] = int(refx), int(refy) + return int(refx), int(refy) + # read and return - refx, refy = mpiops.run_once(__reuse_ref_pixel_file_if_exists) + refx, refy = mpiops.run_once(_reuse_ref_pixel_file_if_exists) if (refx is not None) and (refy is not None): update_refpix_metadata(ifg_paths, int(refx), int(refy), transform, params) return refx, refy @@ -465,22 +488,22 @@ def __reuse_ref_pixel_file_if_exists(): if isinstance(refpixel_returned, ValueError): raise RefPixelError( "Reference pixel calculation returned an all nan slice!\n" - "Cannot continue downstream computation. Please change reference pixel algorithm used before " - "continuing.") + "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 = convert_pixel_value_to_geographic_coordinate(refx, refy, transform) - log.info('Selected reference pixel coordinate (lon, lat): ({}, {})'.format(lon, lat)) + log.info(f'Selected reference pixel coordinate (x, y): ({refx}, {refy})') + lon, lat = convert_pixel_coord_to_lat_lon(refx, refy, transform) + log.info(f'Selected reference pixel coordinate (lon, lat): ({lon}, {lat})') else: - log.info('Using reference pixel from config file (lon, lat): ({}, {})'.format(lon, lat)) + log.info(f'Using reference pixel from config file (lon, lat): ({lon}, {lat})') log.warning("Ensure user supplied reference pixel values are in lon/lat") - refx, refy = convert_geographic_coordinate_to_pixel_value(lon, lat, transform) - log.info('Converted reference pixel coordinate (x, y): ({}, {})'.format(refx, refy)) + refx, refy = convert_lat_lon_to_pixel_coord(lon, lat, transform) + log.info(f'Converted reference pixel coordinate (x, y): ({refx}, {refy})') np.save(file=ref_pixel_file, arr=[int(refx), int(refy)]) update_refpix_metadata(ifg_paths, refx, refy, transform, params) - log.debug("refpx, refpy: "+str(refx) + " " + str(refy)) + log.debug(f"refpx, refpy: {refx} {refy}") ifg.close() params[C.REFX_FOUND], params[C.REFY_FOUND] = int(refx), int(refy) return int(refx), int(refy) diff --git a/pyrate/core/roipac.py b/pyrate/core/roipac.py index d669656e8..98cba8c1d 100644 --- a/pyrate/core/roipac.py +++ b/pyrate/core/roipac.py @@ -115,9 +115,9 @@ def parse_header(hdr_file): if is_dem and DATUM not in headers: msg = 'No "DATUM" parameter in DEM header/resource file' raise RoipacException(msg) - except ValueError: + except ValueError as error: msg = "Unable to parse content of %s. Is it a ROIPAC header file?" - raise RoipacException(msg % hdr_file) + raise RoipacException(msg % hdr_file) from error for k in headers.keys(): if k in INT_HEADERS: @@ -173,14 +173,17 @@ def _parse_dates_from(filename): p = re.compile(r'\d{6}-\d{6}') # match 2 sets of 6 digits separated by '-' m = p.search(filename) - if m: - s = m.group() - min_date_len = 13 # assumes "nnnnnn-nnnnnn" format - if len(s) == min_date_len: - return parse_date(s) - else: # pragma: no cover - msg = "Filename does not include first/second image dates: %s" - raise RoipacException(msg % filename) + if not m: # pragma: no cover + msg = f"Filename does not include first/second image dates: {filename}" + raise RoipacException(msg) + + s = m.group() + min_date_len = 13 # assumes "nnnnnn-nnnnnn" format + if len(s) != min_date_len: # pragma: no cover + msg = f"Filename does not match expected yymmdd-yymmdd format: {filename}" + raise RoipacException(msg) + + return parse_date(s) def manage_header(header_file, projection): @@ -201,13 +204,13 @@ def manage_header(header_file, projection): return header -def roipac_header(file_path, params): +def roipac_header(file_path: Path, params): """ Function to obtain a header for roipac interferogram file or converted geotiff. """ rsc_file = params[C.DEM_HEADER_FILE] - p = Path(file_path) + if rsc_file is not None: projection = parse_header(rsc_file)[ifc.PYRATE_DATUM] else: @@ -216,15 +219,15 @@ def roipac_header(file_path, params): header_file = os.path.join(params[C.DEM_HEADER_FILE]) elif file_path.endswith('unw_ifg.tif') or file_path.endswith('unw.tif'): # TODO: improve this - interferogram_epoches = extract_epochs_from_filename(p.name) + interferogram_epoches = extract_epochs_from_filename(Path(file_path).name) for header_path in params[C.HEADER_FILE_PATHS]: - h = Path(header_path.unwrapped_path) - header_epochs = extract_epochs_from_filename(h.name) + unwrapped_path = Path(header_path.unwrapped_path) + header_epochs = extract_epochs_from_filename(unwrapped_path.name) if set(header_epochs).__eq__(set(interferogram_epoches)): header_file = header_path.unwrapped_path break else: - header_file = "%s%s" % (file_path, ROI_PAC_HEADER_FILE_EXT) + header_file = f"{file_path}{ROI_PAC_HEADER_FILE_EXT}" header = manage_header(header_file, projection) diff --git a/pyrate/core/shared.py b/pyrate/core/shared.py index 7dd48c59f..f36ff9cef 100644 --- a/pyrate/core/shared.py +++ b/pyrate/core/shared.py @@ -20,11 +20,10 @@ """ # pylint: disable=too-many-lines import re -from typing import List, Union, Optional, Iterable, Callable, Dict +from typing import List, Union, Iterable, Callable import errno import math -from joblib import Parallel, delayed from math import floor import os from os.path import basename, join @@ -33,16 +32,17 @@ from datetime import date from itertools import product from enum import Enum +from joblib import Parallel, delayed import numpy as np from numpy import where, nan, isnan, sum as nsum, isclose import pyproj import pkg_resources -import pyrate.constants as C - from osgeo import osr, gdal from osgeo.gdalconst import GA_Update, GA_ReadOnly +import pyrate.constants as C + from pyrate.core import ifgconstants as ifc, mpiops from pyrate.core.logger import pyratelogger as log @@ -66,13 +66,14 @@ class InputTypes(Enum): + """An enum of the types of input files & directories""" IFG = 'ifg' COH = 'coh' BASE = 'base' LT = 'lt' DEM = 'dem' HEADER = 'header' - dir_map = { + DIR_MAP = { IFG: C.INTERFEROGRAM_DIR, COH: C.COHERENCE_DIR, DEM: C.GEOMETRY_DIR, @@ -82,11 +83,11 @@ class InputTypes(Enum): def joblib_log_level(level: str) -> int: """ Convert python log level to joblib int verbosity. - """ + """ if level == 'INFO': return 0 - else: - return 60 + + return 60 def mkdir_p(path): """ @@ -105,7 +106,7 @@ def mkdir_p(path): raise -class RasterBase(object): +class RasterBase: """ Base class for PyRate GeoTIFF based raster datasets. """ @@ -127,22 +128,24 @@ def __init__(self, path: Union[gdal.Dataset, str, Path]): def __str__(self): name = self.__class__.__name__ - return "%s('%s')" % (name, self.data_path) + return f"{name}('{self.data_path}')" def __repr__(self): name = self.__class__.__name__ - return "%s('%s')" % (name, self.data_path) + return f"{name}('{self.data_path}')" def open(self, readonly=None): """ Opens generic raster dataset. """ if self.dataset is not None: - msg = "open() already called for %s" % self + msg = f"open() already called for {self}" raise RasterException(msg) if not os.path.exists(self.data_path): - raise IOError('The file {path} does not exist. Consider first running prepifg'.format(path=self.data_path)) + raise IOError( + f'The file {self.data_path} does not exist. Consider first running prepifg' + ) # unless read only, by default open files as writeable if readonly not in [True, False, None]: @@ -154,7 +157,7 @@ def open(self, readonly=None): flag = GA_ReadOnly if self._readonly else GA_Update self.dataset = gdal.Open(self.data_path, flag) if self.dataset is None: - raise RasterException("Error opening %s" % self.data_path) + raise RasterException(f"Error opening {self.data_path}") self.add_geographic_data() @@ -168,7 +171,9 @@ def add_geographic_data(self): self.lat_centre = self.y_first + (self.y_step * self.y_centre) self.long_centre = self.x_first + (self.x_step * self.x_centre) # use cell size from centre of scene - self.x_size, self.y_size = cell_size(self.lat_centre, self.long_centre, self.x_step, self.y_step) + self.x_size, self.y_size = cell_size( + self.lat_centre, self.long_centre, self.x_step, self.y_step + ) @property def ncols(self): @@ -238,11 +243,11 @@ def num_cells(self): """ Total number of pixels in raster dataset """ - if self.is_open: - return self.dataset.RasterXSize * self.dataset.RasterYSize - else: + if not self.is_open: raise RasterException('Dataset not open') + return self.dataset.RasterXSize * self.dataset.RasterYSize + @property def is_open(self): """ @@ -272,10 +277,10 @@ def _get_band(self, band): :param int band: number of band, starting at 1 """ - if self.dataset is not None: - return self.dataset.GetRasterBand(band) - else: - raise RasterException("Raster %s has not been opened" % self.data_path) + if self.dataset is None: + raise RasterException(f"Raster {self.data_path} has not been opened") + + return self.dataset.GetRasterBand(band) class Ifg(RasterBase): @@ -336,7 +341,7 @@ def _to_date(datestr): self.first, self.second = [_to_date(s) for s in datestrs] self.time_span = (self.second - self.first).days/ifc.DAYS_PER_YEAR else: - msg = 'Missing first and/or second date in %s' % self.data_path + msg = f'Missing first and/or second date in {self.data_path}' raise IfgException(msg) def convert_to_nans(self): @@ -349,21 +354,21 @@ def convert_to_nans(self): 'Use ifg.nodata_value = NoDataValue to set nodata_value' log.warning(msg) raise RasterException(msg) + if ((self.dataset.GetMetadataItem(ifc.NAN_STATUS) == ifc.NAN_CONVERTED) or self.nan_converted): self.phase_data = self.phase_data self.nan_converted = True - msg = '{}: ignored as previous nan ' \ - 'conversion detected'.format(self.data_path) + msg = f'{self.data_path}: ignored as previous nan conversion detected' log.debug(msg) return - else: - self.phase_data = where( - isclose(self.phase_data, self._nodata_value, atol=1e-6), - nan, - self.phase_data) - self.meta_data[ifc.NAN_STATUS] = ifc.NAN_CONVERTED - self.nan_converted = True + + self.phase_data = where( + isclose(self.phase_data, self._nodata_value, atol=1e-6), + nan, + self.phase_data) + self.meta_data[ifc.NAN_STATUS] = ifc.NAN_CONVERTED + self.nan_converted = True # self.write_modified_phase(self.phase_data) @@ -408,10 +413,8 @@ def convert_to_mm(self): """ if self.dataset.GetMetadataItem(ifc.DATA_UNITS) == MILLIMETRES: self.mm_converted = True - msg = '{}: ignored as phase units are already ' \ - 'millimetres'.format(self.data_path) + msg = f'{self.data_path}: ignored as phase units are already millimetres' log.debug(msg) - return elif self.dataset.GetMetadataItem(ifc.DATA_UNITS) == RADIANS: self.phase_data = convert_radians_to_mm(self.phase_data, self.wavelength) self.meta_data[ifc.DATA_UNITS] = MILLIMETRES @@ -420,9 +423,8 @@ def convert_to_mm(self): # otherwise NaN's don't write to bytecode properly # and numpy complains # self.dataset.FlushCache() - msg = '{}: converted phase units to millimetres'.format(self.data_path) + msg = f'{self.data_path}: converted phase units to millimetres' log.debug(msg) - return else: # pragma: no cover msg = 'Phase units are not millimetres or radians' raise IfgException(msg) @@ -437,15 +439,12 @@ def convert_to_radians(self): self.phase_data = convert_mm_to_radians(self.phase_data, wavelength=self.wavelength) self.meta_data[ifc.DATA_UNITS] = RADIANS self.mm_converted = False - msg = '{}: converted phase units to radians'.format(self.data_path) + msg = f'{self.data_path}: converted phase units to radians' log.debug(msg) - return elif self.meta_data[ifc.DATA_UNITS] == RADIANS: self.mm_converted = False - msg = '{}: ignored as phase units are already ' \ - 'radians'.format(self.data_path) + msg = f'{self.data_path}: ignored as phase units are already radians' log.debug(msg) - return else: # pragma: no cover msg = 'Phase units are not millimetres or radians' raise IfgException(msg) @@ -516,6 +515,7 @@ def write_modified_phase(self, data=None): self.dataset.FlushCache() def add_metadata(self, **kwargs): + """Adds metadata to the interferogram's geotiff file""" if (not self.is_open) or self.is_read_only: raise IOError("Ifg not open or readonly. Cannot write!") @@ -550,7 +550,7 @@ def __str__(self): return "Convenience Tile class containing tile co-ordinates" -class IfgPart(object): +class IfgPart: """ Create a tile (subset) of an Ifg data object """ @@ -571,7 +571,7 @@ def __init__(self, ifg_or_path, tile: Tile, ifg_dict=None, params=None): self.first = ifg.first self.second = ifg.second self.time_span = ifg.time_span - phase_file = 'phase_data_{}_{}.npy'.format(basename(ifg_or_path).split('.')[0], tile.index) + phase_file = f'phase_data_{basename(ifg_or_path).split(".")[0]}_{tile.index}.npy' self.phase_data = np.load(join(params[C.TMPDIR], phase_file)) else: # check if Ifg was sent. @@ -674,6 +674,7 @@ def azimuth_data(self): class TileMixin: + """A mixin class that makes the inheriting class callable to subscript by a tile region""" def __call__(self, tile: Tile): t = tile return self.data[t.top_left_y:t.bottom_right_y, t.top_left_x:t.bottom_right_x] @@ -714,7 +715,7 @@ def data(self): class MemGeometry(TileMixin): - + """A simple class which holds geometry data in memory""" def __init__(self, data): """ Set phase data value @@ -723,6 +724,9 @@ def __init__(self, data): @property def data(self): + """ + Returns the geometry band as an array. + """ return self._data @@ -741,7 +745,7 @@ class RasterException(Exception): """ -class EpochList(object): +class EpochList: """ Metadata container for epoch related information. """ @@ -755,10 +759,10 @@ def __init__(self, dates=None, repeat=None, spans=None): self.spans = spans # time span from earliest ifg def __str__(self): - return "EpochList: %s" % str(self.dates) + return f"EpochList: {str(self.dates)}" def __repr__(self): - return "EpochList: %s" % repr(self.dates) + return f"EpochList: {repr(self.dates)}" def convert_radians_to_mm(data, wavelength): @@ -801,8 +805,9 @@ def nanmedian(x): version = [int(i) for i in pkg_resources.get_distribution("numpy").version.split('.')[:2]] if version[0] == 1 and version[1] > 9: return np.nanmedian(x) - else: # pragma: no cover - return np.median(x[~np.isnan(x)]) + + # pragma: no cover + return np.median(x[~np.isnan(x)]) def _is_interferogram(hdr): @@ -870,23 +875,32 @@ def write_fullres_geotiff(header, data_path, dest, nodata): _check_raw_data(bytes_per_col, data_path, ncols, nrows) # position and projection data - gt = [header[ifc.PYRATE_LONG], header[ifc.PYRATE_X_STEP], 0, header[ifc.PYRATE_LAT], 0, header[ifc.PYRATE_Y_STEP]] + gt = [ + header[ifc.PYRATE_LONG], header[ifc.PYRATE_X_STEP], 0, + header[ifc.PYRATE_LAT], 0, header[ifc.PYRATE_Y_STEP] + ] srs = osr.SpatialReference() res = srs.SetWellKnownGeogCS(header[ifc.PYRATE_DATUM]) if res: - msg = 'Unrecognised projection: %s' % header[ifc.PYRATE_DATUM] + msg = f'Unrecognised projection: {header[ifc.PYRATE_DATUM]}' raise GeotiffException(msg) wkt = srs.ExportToWkt() - dtype = 'float32' if (_is_interferogram(header) or _is_incidence(header) or _is_coherence(header) or \ - _is_baseline(header) or _is_lookuptable(header)) else 'int16' + is_float = _is_interferogram(header) + is_float |= _is_incidence(header) + is_float |= _is_coherence(header) + is_float |= _is_baseline(header) + is_float |= _is_lookuptable(header) + dtype = 'float32' if is_float else 'int16' # get subset of metadata relevant to PyRate md = collate_metadata(header) # create GDAL object ds = gdal_dataset( - dest, ncols, nrows, driver="GTiff", bands=1, dtype=dtype, metadata=md, crs=wkt, geotransform=gt, + dest, ncols, nrows, + driver="GTiff", bands=1, dtype=dtype, + metadata=md, crs=wkt, geotransform=gt, creation_opts=["compress=packbits"] ) @@ -947,7 +961,7 @@ def collate_metadata(header): :return: dict of relevant metadata for PyRate """ - md = dict() + md = {} def __common_ifg_coh_update(header, md): for k in [ifc.PYRATE_WAVELENGTH_METRES, ifc.PYRATE_TIME_SPAN, @@ -1006,7 +1020,7 @@ def data_format(ifg_proc, is_ifg, ncols): fmtstr = '<' + ('h' * ncols) # roipac DEM is little endian signed int16 bytes_per_col = 2 else: # pragma: no cover - msg = 'Unrecognised InSAR Processor: %s' % ifg_proc + msg = f'Unrecognised InSAR Processor: {ifg_proc}' raise GeotiffException(msg) return bytes_per_col, fmtstr @@ -1080,10 +1094,14 @@ def write_output_geotiff(md, gt, wkt, data, dest, nodata): ds.SetMetadataItem(ifc.DATA_TYPE, str(md[ifc.DATA_TYPE])) # set other metadata - 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, - ifc.EPOCH_DATE, C.LOS_PROJECTION.upper(), C.SIGNAL_POLARITY.upper(), - C.VELERROR_NSIG.upper()]: + metadata = [ + 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, + ifc.EPOCH_DATE, C.LOS_PROJECTION.upper(), C.SIGNAL_POLARITY.upper(), + C.VELERROR_NSIG.upper() + ] + + for k in metadata: if k in md: ds.SetMetadataItem(k, str(md[k])) @@ -1108,9 +1126,9 @@ def write_geotiff(data, outds, nodata): """ # only support "2 <= dims <= 3" if data.ndim == 3: - count, height, width = data.shape + _, _, _ = data.shape elif data.ndim == 2: - height, width = data.shape + _, _ = data.shape else: msg = "Only support dimensions of '2 <= dims <= 3'." raise GeotiffException(msg) @@ -1156,9 +1174,11 @@ def create_tiles(shape, nrows=2, ncols=2): if ncols > no_x or nrows > no_y: raise ValueError('nrows/cols must be greater than ifg dimensions') - col_arr = np.array_split(range(no_x), ncols) - row_arr = np.array_split(range(no_y), nrows) - return [Tile(i, (r[0], c[0]), (r[-1]+1, c[-1]+1)) for i, (r, c) in enumerate(product(row_arr, col_arr))] + cols = np.array_split(range(no_x), ncols) + rows = np.array_split(range(no_y), nrows) + return [ + Tile(i, (r[0], c[0]), (r[-1]+1, c[-1]+1)) for i, (r, c) in enumerate(product(rows, cols)) + ] def get_tiles(ifg_path, rows, cols) -> List[Tile]: @@ -1280,7 +1300,7 @@ def save_numpy_phase(ifg_paths, params): for t in tiles: p_data = phase_data[t.top_left_y:t.bottom_right_y, t.top_left_x:t.bottom_right_x] - phase_file = 'phase_data_{}_{}.npy'.format(bname, t.index) + phase_file = f'phase_data_{bname}_{t.index}.npy' np.save(file=join(outdir, phase_file), arr=p_data) ifg.close() @@ -1341,7 +1361,7 @@ def check_correction_status(ifgs, meta): # pragma: no cover def close_all(ifgs): for ifg in ifgs: ifg.close() - + if not isinstance(ifgs[0], Ifg): ifgs = [Ifg(ifg_path) for ifg_path in ifgs] @@ -1353,20 +1373,21 @@ def close_all(ifgs): if all(flags): log.info('Skipped: interferograms already corrected') return True - elif not all(flags) and any(flags): + + if not all(flags) and any(flags): log.debug('Detected mix of corrected and uncorrected interferograms') - for i, flag in zip(ifgs, flags): + for flag in flags: if flag: - msg = '{}: correction detected'.format(i.data_path) + msg = '{i.data_path}: correction detected' else: - msg = '{}: correction NOT detected'.format(i.data_path) + msg = '{i.data_path}: correction NOT detected' log.debug(msg) close_all(ifgs) raise CorrectionStatusError(msg) - else: - log.debug('Calculating corrections') - close_all(ifgs) - return False + + log.debug('Calculating corrections') + close_all(ifgs) + return False class CorrectionStatusError(Exception): @@ -1376,6 +1397,7 @@ class CorrectionStatusError(Exception): def extract_epochs_from_filename(filename_with_epochs: str) -> List[str]: + """Extract the 8 or 6 digit epochs from IFG filenames""" src_epochs = re.findall(r"(\d{8})", str(filename_with_epochs)) if not len(src_epochs) > 0: src_epochs = re.findall(r"(\d{6})", str(filename_with_epochs)) @@ -1383,6 +1405,7 @@ def extract_epochs_from_filename(filename_with_epochs: str) -> List[str]: def mpi_vs_multiprocess_logging(step, params): + """Logs the state of job processing (MPI vs. parallel)""" if mpiops.size > 1: # Over-ride input options if this is an MPI job log.info(f"Running '{step}' step with MPI using {mpiops.size} processes") log.warning("Disabling joblib parallel processing (setting parallel = 0)") @@ -1405,12 +1428,11 @@ def dem_or_ifg(data_path: str) -> Union[Ifg, DEM]: """ ds = gdal.Open(data_path) md = ds.GetMetadata() - if ifc.FIRST_DATE in md: # ifg - return Ifg(data_path) - else: - return DEM(data_path) + is_ifg = ifc.FIRST_DATE in md ds = None # close the dataset + return Ifg(data_path) if is_ifg else DEM(data_path) + def join_dicts(dicts: List[dict]) -> dict: """ @@ -1425,9 +1447,10 @@ def join_dicts(dicts: List[dict]) -> dict: def iterable_split(func: Callable, iterable: Iterable, params: dict, *args, **kwargs) -> np.ndarray: """ # TODO: a faster version using buffer-provider objects via the uppercase communication method - A faster version of iterable/tiles_split is possible when the return values from each process is of the same size - and will be addressed in future. In this case a buffer-provider object can be sent between processes using the - uppercase communication (like Gather instead of gather) methods which can be significantly faster. + A faster version of iterable/tiles_split is possible when the return values from each process + is of the same size and will be addressed in future. In this case a buffer-provider object can + be sent between processes using the uppercase communication (like Gather instead of gather) + methods which can be significantly faster. """ if params[C.PARALLEL]: ret_combined = {} diff --git a/pyrate/core/stack.py b/pyrate/core/stack.py index c2443b785..dffff262d 100644 --- a/pyrate/core/stack.py +++ b/pyrate/core/stack.py @@ -21,7 +21,7 @@ import os import pickle as cp from scipy.linalg import solve, cholesky, qr, inv -from numpy import nan, isnan, sqrt, diag, delete, array, float32, size +from numpy import nan, isnan, sqrt, diag, delete, array, float32 import numpy as np import pyrate.constants as C @@ -49,13 +49,34 @@ def stack_rate_array(ifgs, params, vcmt, mst=None): :return: samples: Number of observations used in rate calculation for each pixel :rtype: ndarray """ - nsig, pthresh, cols, error, mst, obs, rate, rows, samples, span = _stack_setup(ifgs, mst, params) + + # stack rate parameters from config file + # n-sigma ratio used to threshold 'model minus observation' residuals + nsig = params[C.LR_NSIG] + # Pixel threshold; minimum number of observations for a pixel + pthresh = params[C.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) # pixel-by-pixel calculation. # nested loops to loop over the 2 image dimensions for i in range(rows): for j in range(cols): - rate[i, j], error[i, j], samples[i, j] = stack_rate_pixel(obs[:, i, j], mst[:, i, j], vcmt, span, nsig, pthresh) + value = stack_rate_pixel(obs[:, i, j], mst[:, i, j], vcmt, span, nsig, pthresh) + + rate[i, j], error[i, j], samples[i, j] = value return rate, params[C.VELERROR_NSIG]*error, samples @@ -84,7 +105,7 @@ def mask_rate(rate, error, maxsig): error[mask] = nan # calculate percentage of masked pixels nummasked = int(np.count_nonzero(mask)/orig*100) - log.info('Percentage of pixels masked = {}%'.format(nummasked)) + log.info(f'Percentage of pixels masked = {nummasked}%') return rate, error @@ -108,13 +129,17 @@ def stack_rate_pixel(obs, mst, vcmt, span, nsig, pthresh): :return: samples: Number of observations used in the rate estimation for the pixel :rtype: int """ - # 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: + # pylint: disable=invalid-name + # JUSTIFICATION: This code is implememnting a mathematical algorithm where these + # variable names likely match up to standard variable names used in mathematics or + # the relevant paper it was originally published in. + # select ifg observations ifgv = obs[ind] @@ -165,32 +190,6 @@ def stack_rate_pixel(obs, mst, vcmt, span, nsig, pthresh): return np.nan, np.nan, default_no_samples -def _stack_setup(ifgs, mst, params): - """ - Convenience function for stack rate setup - """ - # stack rate parameters from config file - # n-sigma ratio used to threshold 'model minus observation' residuals - nsig = params[C.LR_NSIG] - # Pixel threshold; minimum number of observations for a pixel - pthresh = params[C.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, rate, rows, samples, span - - def stack_calc_wrapper(params): """ Wrapper for stacking on a set of tiles. @@ -217,6 +216,6 @@ def _stacking_for_tile(tile, params): ifg_parts = [shared.IfgPart(p, tile, preread_ifgs, params) for p in ifg_paths] mst_tile = np.load(Configuration.mst_path(params, tile.index)) rate, error, samples = stack_rate_array(ifg_parts, params, vcmt, mst_tile) - np.save(file=os.path.join(output_dir, 'stack_rate_{}.npy'.format(tile.index)), arr=rate) - np.save(file=os.path.join(output_dir, 'stack_error_{}.npy'.format(tile.index)), arr=error) - np.save(file=os.path.join(output_dir, 'stack_samples_{}.npy'.format(tile.index)), arr=samples) + np.save(file=os.path.join(output_dir, f'stack_rate_{tile.index}.npy'), arr=rate) + np.save(file=os.path.join(output_dir, f'stack_error_{tile.index}.npy'), arr=error) + np.save(file=os.path.join(output_dir, f'stack_samples_{tile.index}.npy'), arr=samples) diff --git a/pyrate/core/timeseries.py b/pyrate/core/timeseries.py index ad50edd0a..71f21f37d 100644 --- a/pyrate/core/timeseries.py +++ b/pyrate/core/timeseries.py @@ -163,6 +163,8 @@ def _remove_rank_def_rows(b_mat, nvelpar, ifgv, sel): """ _, _, e_var = qr(b_mat, mode='economic', pivoting=True) licols = e_var[matrix_rank(b_mat):nvelpar] + # pylint: disable=unbalanced-tuple-unpacking + # JUSTIFICATION: pylint is wrong here (due to numpy's python stub code being wrong) [rmrow, _] = where(b_mat[:, licols] != 0) b_mat = delete(b_mat, rmrow, axis=0) ifgv = delete(ifgv, rmrow) @@ -312,19 +314,19 @@ def linear_rate_pixel(y, t): # Mask to exclude nan elements mask = ~isnan(y) # remove nan elements from both arrays - y = y[mask] + y = y[mask] try: t = t[mask] - except IndexError: - raise TimeSeriesError("linear_rate_pixel: y and t are not equal length") + except IndexError as error: + raise TimeSeriesError("linear_rate_pixel: y and t are not equal length") from error # break out of func if not enough time series obs for line fitting nsamp = len(y) if nsamp < 2: return nan, nan, nan, nan, nan - # compute linear regression of tscuml - linrate, intercept, r_value, p_value, std_err = linregress(t, y) + # compute linear regression of tscuml + linrate, intercept, r_value, _, std_err = linregress(t, y) return linrate, intercept, r_value**2, std_err, int(nsamp) @@ -350,9 +352,7 @@ def linear_rate_array(tscuml, ifgs, params): :return: samples: Number of observations used in linear regression for each pixel :rtype: ndarray """ - b0_mat, interp, p_thresh, sm_factor, sm_order, ts_method, ifg_data, mst, \ - ncols, nrows, nvelpar, span, tsvel_matrix = \ - _time_series_setup(ifgs, params) + _, _, _, _, _, _, _, _, ncols, nrows, _, _, _ = _time_series_setup(ifgs, params) epochlist = get_epochs(ifgs)[0] # get cumulative time per epoch @@ -383,7 +383,7 @@ def _missing_option_error(option): """ Convenience function for raising similar missing option errors. """ - msg = "Missing '%s' option in config file" % option + msg = f"Missing '{option}' option in config file" raise ConfigException(msg) @@ -422,16 +422,15 @@ def __calc_time_series_for_tile(tile, params): ifg_parts = [shared.IfgPart(p, tile, preread_ifgs, params) for p in ifg_paths] mst_tile = np.load(Configuration.mst_path(params, tile.index)) tsincr, tscuml, _ = time_series(ifg_parts, params, vcmt, mst_tile) - np.save(file=os.path.join(output_dir, 'tscuml_{}.npy'.format(tile.index)), arr=tscuml) + np.save(file=os.path.join(output_dir, f'tscuml_{tile.index}.npy'), arr=tscuml) # optional save of tsincr npy tiles if params["savetsincr"] == 1: - np.save(file=os.path.join(output_dir, 'tsincr_{}.npy'.format(tile.index)), arr=tsincr) + np.save(file=os.path.join(output_dir, f'tsincr_{tile.index}.npy'), arr=tsincr) tscuml = np.insert(tscuml, 0, 0, axis=2) # add zero epoch to tscuml 3D array log.info('Calculating linear regression of cumulative time series') linrate, intercept, r_squared, std_err, samples = linear_rate_array(tscuml, ifg_parts, params) - np.save(file=os.path.join(output_dir, 'linear_rate_{}.npy'.format(tile.index)), arr=linrate) - np.save(file=os.path.join(output_dir, 'linear_intercept_{}.npy'.format(tile.index)), arr=intercept) - np.save(file=os.path.join(output_dir, 'linear_rsquared_{}.npy'.format(tile.index)), arr=r_squared) - np.save(file=os.path.join(output_dir, 'linear_error_{}.npy'.format(tile.index)), arr=std_err) - np.save(file=os.path.join(output_dir, 'linear_samples_{}.npy'.format(tile.index)), arr=samples) - + np.save(file=os.path.join(output_dir, f'linear_rate_{tile.index}.npy'), arr=linrate) + np.save(file=os.path.join(output_dir, f'linear_intercept_{tile.index}.npy'), arr=intercept) + np.save(file=os.path.join(output_dir, f'linear_rsquared_{tile.index}.npy'), arr=r_squared) + np.save(file=os.path.join(output_dir, f'linear_error_{tile.index}.npy'), arr=std_err) + np.save(file=os.path.join(output_dir, f'linear_samples_{tile.index}.npy'), arr=samples) diff --git a/pyrate/correct.py b/pyrate/correct.py index 069c416df..0ba88db16 100644 --- a/pyrate/correct.py +++ b/pyrate/correct.py @@ -21,8 +21,6 @@ import os from pathlib import Path import pickle as cp -from typing import List -import sys import pyrate.constants as C from pyrate.core import (shared, algorithm, mpiops) @@ -31,14 +29,14 @@ from pyrate.core.mst import mst_calc_wrapper from pyrate.core.orbital import orb_fit_calc_wrapper from pyrate.core.dem_error import dem_error_calc_wrapper -from pyrate.core.phase_closure.closure_check import iterative_closure_check, mask_pixels_with_unwrapping_errors, \ - update_ifg_list +from pyrate.core.phase_closure.closure_check import iterative_closure_check, \ + mask_pixels_with_unwrapping_errors, update_ifg_list from pyrate.core.ref_phs_est import ref_phase_est_wrapper from pyrate.core.refpixel import ref_pixel_calc_wrapper -from pyrate.core.shared import PrereadIfg, Ifg, get_tiles, mpi_vs_multiprocess_logging, join_dicts, \ - nan_and_mm_convert, save_numpy_phase +from pyrate.core.shared import PrereadIfg, Ifg, get_tiles, mpi_vs_multiprocess_logging, \ + join_dicts, nan_and_mm_convert, save_numpy_phase from pyrate.core.logger import pyratelogger as log -from pyrate.configuration import Configuration, MultiplePaths, ConfigException +from pyrate.configuration import Configuration, ConfigException MAIN_PROCESS = 0 @@ -56,7 +54,7 @@ def _create_ifg_dict(params): interferograms that are used later in workflow :rtype: dict """ - dest_tifs = [ifg_path for ifg_path in params[C.INTERFEROGRAM_FILES]] + dest_tifs = list(params[C.INTERFEROGRAM_FILES]) ifgs_dict = {} process_tifs = mpiops.array_split(dest_tifs) for d in process_tifs: @@ -78,13 +76,16 @@ def _create_ifg_dict(params): ifg.close() ifgs_dict = join_dicts(mpiops.comm.allgather(ifgs_dict)) - ifgs_dict = mpiops.run_once(__save_ifgs_dict_with_headers_and_epochs, dest_tifs, ifgs_dict, params, process_tifs) + ifgs_dict = mpiops.run_once( + _save_ifgs_dict_with_headers_and_epochs, + dest_tifs, ifgs_dict, params, process_tifs + ) params[C.PREREAD_IFGS] = ifgs_dict return ifgs_dict -def __save_ifgs_dict_with_headers_and_epochs(dest_tifs, ifgs_dict, params, process_tifs): +def _save_ifgs_dict_with_headers_and_epochs(dest_tifs, ifgs_dict, params, process_tifs): tmpdir = params[C.TMPDIR] if not os.path.exists(tmpdir): shared.mkdir_p(tmpdir) @@ -94,13 +95,14 @@ def __save_ifgs_dict_with_headers_and_epochs(dest_tifs, ifgs_dict, params, proce # add some extra information that's also useful later gt, md, wkt = shared.get_geotiff_header_info(process_tifs[0].tmp_sampled_path) epochlist = algorithm.get_epochs(ifgs_dict)[0] - log.info('Found {} unique epochs in the {} interferogram network'.format(len(epochlist.dates), nifgs)) + log.info(f'Found {len(epochlist.dates)} unique epochs in the {nifgs} interferogram network') ifgs_dict['epochlist'] = epochlist ifgs_dict['gt'] = gt ifgs_dict['md'] = md ifgs_dict['wkt'] = wkt # dump ifgs_dict file for later use - cp.dump(ifgs_dict, open(preread_ifgs_file, 'wb')) + with open(preread_ifgs_file, 'wb') as file: + cp.dump(ifgs_dict, file) for k in ['gt', 'epochlist', 'md', 'wkt']: ifgs_dict.pop(k) @@ -116,9 +118,10 @@ def _copy_mlooked(params): log.info("Copying input files into tempdir for manipulation during 'correct' steps") mpaths = params[C.INTERFEROGRAM_FILES] process_mpaths = mpiops.array_split(mpaths) - for p in process_mpaths: - shutil.copy(p.sampled_path, p.tmp_sampled_path) - Path(p.tmp_sampled_path).chmod(0o664) # assign write permission as prepifg output is readonly + for path in process_mpaths: + shutil.copy(path.sampled_path, path.tmp_sampled_path) + # assign write permission as prepifg output is readonly + Path(path.tmp_sampled_path).chmod(0o664) def main(config): @@ -155,21 +158,21 @@ def phase_closure_wrapper(params: dict, config: Configuration) -> dict: This wrapper will run the iterative phase closure check to return a stable list of checked interferograms, and then mask pixels in interferograms that exceed the unwrapping error threshold. - :param params: Dictionary of PyRate configuration parameters. + :param params: Dictionary of PyRate configuration parameters. :param config: Configuration class instance. :return: params: Updated dictionary of PyRate configuration parameters. """ if not params[C.PHASE_CLOSURE]: log.info("Phase closure correction is not required!") - return + return None rets = iterative_closure_check(config) if rets is None: log.info("Zero loops are returned from the iterative closure check.") log.warning("Abandoning phase closure correction without modifying the interferograms.") - return - + return None + ifg_files, ifgs_breach_count, num_occurences_each_ifg = rets # update params with closure checked ifg list @@ -177,7 +180,7 @@ def phase_closure_wrapper(params: dict, config: Configuration) -> dict: mpiops.run_once(update_ifg_list, ifg_files, params[C.INTERFEROGRAM_FILES]) if mpiops.rank == 0: - with open(config.phase_closure_filtered_ifgs_list(params), 'w') as f: + with open(config.phase_closure_filtered_ifgs_list(params), 'w', encoding="utf-8") as f: lines = [p.converted_path + '\n' for p in params[C.INTERFEROGRAM_FILES]] f.writelines(lines) @@ -236,6 +239,6 @@ def correct_ifgs(config: Configuration) -> None: def __validate_correct_steps(params): for step in params['correct']: - if step not in correct_steps.keys(): + if step not in correct_steps: raise ConfigException(f"{step} is not a supported 'correct' step. \n" f"Supported steps are {correct_steps.keys()}") diff --git a/pyrate/default_parameters.py b/pyrate/default_parameters.py index 0a41c3161..34237931f 100644 --- a/pyrate/default_parameters.py +++ b/pyrate/default_parameters.py @@ -551,7 +551,9 @@ }, "correct": { "DataType": list, - "DefaultValue": ['orbfit', 'refphase', 'demerror', 'mst', 'apscorrect', 'maxvar', 'timeseries', 'stack'], + "DefaultValue": [ + 'orbfit', 'refphase', 'demerror', 'mst', 'apscorrect', 'maxvar', 'timeseries', 'stack' + ], "MinValue": None, "MaxValue": None, "PossibleValues": None, diff --git a/pyrate/main.py b/pyrate/main.py index e19bf0970..febe13358 100644 --- a/pyrate/main.py +++ b/pyrate/main.py @@ -50,51 +50,87 @@ def update_params_due_to_ifg_selection(config): return params +# pylint: disable=too-many-statements +# JUSTIFICATION: extracting statements into separate functions would make it harder to understand. def main(): + """The main PyRate runner program, refer to `docs/index.rst` and `--help`""" start_time = time.time() - parser = argparse.ArgumentParser(prog='pyrate', description=CLI_DESCRIPTION, add_help=True, - formatter_class=RawTextHelpFormatter) - parser.add_argument('-v', '--verbosity', type=str, default='INFO', choices=['DEBUG', 'INFO', 'WARNING', 'ERROR'], - help="Increase output verbosity") + parser = argparse.ArgumentParser( + prog='pyrate', description=CLI_DESCRIPTION, add_help=True, + formatter_class=RawTextHelpFormatter + ) + parser.add_argument( + '-v', '--verbosity', type=str, + default='INFO', choices=['DEBUG', 'INFO', 'WARNING', 'ERROR'], + help="Increase output verbosity" + ) subparsers = parser.add_subparsers(dest='command') subparsers.required = True - parser_conv2tif = subparsers.add_parser('conv2tif', help=' Convert interferograms to geotiff.', - add_help=True) + parser_conv2tif = subparsers.add_parser( + 'conv2tif', + help=' Convert interferograms to geotiff.', + add_help=True + ) parser_prepifg = subparsers.add_parser( - 'prepifg', help='Perform multilooking, cropping and coherence masking to interferogram geotiffs.', - add_help=True) + 'prepifg', + help='Perform multilooking, cropping and coherence masking to interferogram geotiffs.', + add_help=True + ) parser_correct = subparsers.add_parser( - 'correct', help='Calculate and apply corrections to interferogram phase data.', - add_help=True) + 'correct', + help='Calculate and apply corrections to interferogram phase data.', + add_help=True + ) parser_ts = subparsers.add_parser( - 'timeseries', help=' Timeseries inversion of interferogram phase data.', add_help=True + 'timeseries', + help=' Timeseries inversion of interferogram phase data.', + add_help=True ) - parser_stack = subparsers.add_parser('stack', help=' Stacking of interferogram phase data.', - add_help=True) + parser_stack = subparsers.add_parser( + 'stack', + help=' Stacking of interferogram phase data.', + add_help=True + ) parser_merge = subparsers.add_parser( - 'merge', help="Reassemble computed tiles and save as geotiffs.", - add_help=True) + 'merge', + help="Reassemble computed tiles and save as geotiffs.", + add_help=True + ) parser_workflow = subparsers.add_parser( - 'workflow', help=" Sequentially run all the PyRate processing steps.", - add_help=True) - for p in [parser_conv2tif, parser_prepifg, parser_correct, parser_merge, parser_ts, parser_stack, parser_workflow]: - p.add_argument('-f', '--config_file', action="store", type=str, default=None, + 'workflow', + help=" Sequentially run all the PyRate processing steps.", + add_help=True + ) + + sub_parsers = [ + parser_conv2tif, + parser_prepifg, + parser_correct, + parser_merge, + parser_ts, + parser_stack, + parser_workflow + ] + + for sub_parser in sub_parsers: + sub_parser.add_argument('-f', '--config_file', action="store", type=str, default=None, help="Pass configuration file", required=False) args = parser.parse_args() params = mpiops.run_once(_params_from_conf, args.config_file) - configure_stage_log(args.verbosity, args.command, Path(params[C.OUT_DIR]).joinpath('pyrate.log.').as_posix()) + log_path = Path(params[C.OUT_DIR]).joinpath('pyrate.log.').as_posix() + configure_stage_log(args.verbosity, args.command, log_path) log.debug("Starting PyRate") log.debug("Arguments supplied at command line: ") @@ -102,7 +138,7 @@ def main(): if args.verbosity: log.setLevel(args.verbosity) - log.info("Verbosity set to " + str(args.verbosity) + ".") + log.info(f"Verbosity set to {args.verbosity}.") if args.command == "conv2tif": conv2tif.main(params) @@ -154,10 +190,12 @@ def main(): params = mpiops.run_once(_params_from_conf, args.config_file) merge.main(params) - log.info("--- Runtime = %s seconds ---" % (time.time() - start_time)) + elapsed_secs = time.time() - start_time + log.info(f"--- Runtime = {elapsed_secs} seconds ---") def timeseries(config: Configuration) -> None: + """The runner command for calculating the timeseries of a file set""" params = config.__dict__ mpi_vs_multiprocess_logging("timeseries", params) params = update_params_due_to_ifg_selection(config=config) @@ -165,6 +203,7 @@ def timeseries(config: Configuration) -> None: def stack(config: Configuration) -> None: + """The runner command for stacking a set of files""" params = config.__dict__ mpi_vs_multiprocess_logging("stack", params) params = update_params_due_to_ifg_selection(config=config) diff --git a/pyrate/merge.py b/pyrate/merge.py index cba6be3da..e91f20415 100644 --- a/pyrate/merge.py +++ b/pyrate/merge.py @@ -18,12 +18,12 @@ stack rate and time series outputs and save as geotiff files """ from os.path import join, isfile, exists -import pickle -import numpy as np -from osgeo import gdal import subprocess from pathlib import Path from typing import Optional, Tuple +import pickle +import numpy as np +from osgeo import gdal import pyrate.constants as C from pyrate.core import shared, stack, ifgconstants as ifc, mpiops @@ -40,11 +40,11 @@ def main(params: dict) -> None: single geotiff files """ if params[C.SIGNAL_POLARITY] == 1: - log.info(f"Saving output products with the same sign convention as input data") + log.info("Saving output products with the same sign convention as input data") elif params[C.SIGNAL_POLARITY] == -1: - log.info(f"Saving output products with reversed sign convention compared to the input data") + log.info("Saving output products with reversed sign convention compared to the input data") else: - log.warning(f"Check the value of signal_polarity parameter") + log.warning("Check the value of signal_polarity parameter") out_types = [] @@ -58,7 +58,7 @@ def main(params: dict) -> None: if params["savetsincr"] == 1: _merge_timeseries(params, 'tsincr') else: - log.warning('Not merging time series products; {} does not exist'.format(tsfile)) + log.warning(f'Not merging time series products; {tsfile} does not exist') stfile = join(params[C.TMPDIR], 'stack_rate_0.npy') if exists(stfile): @@ -66,7 +66,7 @@ def main(params: dict) -> None: mpiops.run_once(_merge_stack, params) out_types += ['stack_rate', 'stack_error'] else: - log.warning('Not merging stack products; {} does not exist'.format(stfile)) + log.warning(f'Not merging stack products; {stfile} does not exist') if len(out_types) > 0: process_out_types = mpiops.array_split(out_types) @@ -91,14 +91,16 @@ def _merge_stack(params: dict) -> None: # mask pixels according to threshold if params[C.LR_MAXSIG] > 0: - log.info(f"Masking stack_rate and stack_error maps where error is greater than {params[C.LR_MAXSIG]} millimetres") + msg = f"Masking stack_rate and stack_error maps where error is greater " \ + f"than {params[C.LR_MAXSIG]} millimetres" + log.info(msg) rate, error = stack.mask_rate(rate, error, params[C.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, out, ot, savenpy=params["savenpy"]) + for out, otype in zip([rate, error, samples], ['stack_rate', 'stack_error', 'stack_samples']): + __save_merged_files(ifgs_dict, params, out, otype, savenpy=params["savenpy"]) def _merge_linrate(params: dict) -> None: @@ -122,7 +124,7 @@ def _merge_timeseries(params: dict, tstype: str) -> None: """ Merge tiled time series outputs """ - log.info('Merging {} time series outputs'.format(tstype)) + log.info(f'Merging {tstype} time series outputs') shape, tiles, ifgs_dict = __merge_setup(params) # load the first time series file to determine the number of time series tifs @@ -135,15 +137,17 @@ def _merge_timeseries(params: dict, tstype: str) -> None: # 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 {} {} time series tifs of ' - 'total {}'.format(mpiops.rank, len(process_tifs), tstype, no_ts_tifs)) + log.info(f'Process {mpiops.rank} writing {len(process_tifs)} {tstype} ' \ + f'time series tifs of total {no_ts_tifs}') for i in process_tifs: ts_arr = assemble_tiles(shape, params[C.TMPDIR], tiles, out_type=tstype, index=i) - __save_merged_files(ifgs_dict, params, ts_arr, out_type=tstype, index=i, savenpy=params["savenpy"]) + __save_merged_files( + ifgs_dict, params, ts_arr, out_type=tstype, index=i, savenpy=params["savenpy"] + ) mpiops.comm.barrier() - log.debug('Process {} finished writing {} {} time series tifs of ' - 'total {}'.format(mpiops.rank, len(process_tifs), tstype, no_ts_tifs)) + log.debug(f'Process {mpiops.rank} finished writing {len(process_tifs)} {tstype} time series ' + f'tifs of total {no_ts_tifs}') def create_png_and_kml_from_tif(output_folder_path: str, output_type: str) -> None: @@ -183,21 +187,22 @@ def create_png_and_kml_from_tif(output_folder_path: str, output_type: str) -> No """ - with open(kml_file_path, "w") as f: + with open(kml_file_path, "w", encoding="utf-8") as f: f.write(kml_file_content) # Get raster statistics 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) + # 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 in ('stack_rate', 'linear_rate'): # minimum value might be negative maximum = max(abs(minimum), abs(maximum)) minimum = -1 * maximum - # colours: blue -> white -> red (white==0) + # 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) @@ -212,56 +217,63 @@ def create_png_and_kml_from_tif(output_folder_path: str, output_type: str) -> No g = np.flipud(g) * 255 b = np.flipud(b) * 255 if output_type in ('stack_error', 'linear_error', 'linear_rsquared'): - # colours: white -> red (minimum error -> maximum 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)) + cmap_path = join(output_folder_path, f"colourmap_{output_type}.txt") + log.info(f'Saving colour map to file {cmap_path}; min/max values: {minimum:.2f}/{maximum:.2f}') # write colourmap as text file - with open(color_map_path, "w") as f: + with open(cmap_path, "w", encoding="utf-8") as f: 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])) + f.write(f"{value} {r[i]} {g[i]} {b[i]} 255\n") # create PNG image file input_tif_path = join(output_folder_path, f"{output_type}.tif") output_png_path = join(output_folder_path, f"{output_type}.png") subprocess.check_call(["gdaldem", "color-relief", "-of", "PNG", input_tif_path, "-alpha", - color_map_path, output_png_path, "-nearest_color_entry"]) + cmap_path, output_png_path, "-nearest_color_entry"]) log.debug(f'Finished creating quicklook image for {output_type}') -def assemble_tiles(s: Tuple, dir: str, tiles: Tile, out_type: str, index: Optional[int] = None) -> np.ndarray: +def assemble_tiles( + shape: Tuple, + tile_dir: str, + tiles: Tile, + out_type: str, + index: Optional[int] = None +) -> np.ndarray: """ Function to reassemble tiles from numpy files in to a merged array :param s: shape for merged array. - :param dir: path to directory containing numpy tile files. + :param tile_dir: path to directory containing numpy tile files. :param tiles: pyrate.core.shared.Tile Class object. :param out_type: product type string, used to construct numpy tile file name. :param index: array third dimension index to extract from 3D time series array tiles. :return: merged_array: array assembled from all tiles. """ - log.debug('Re-assembling tiles for {}'.format(out_type)) + log.debug(f'Re-assembling tiles for {out_type}') # pre-allocate dest array - merged_array = np.empty(shape=s, dtype=np.float32) + merged_array = np.empty(shape=shape, dtype=np.float32) # 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) + for tile in tiles: + tile_file = Path(join(tile_dir, out_type + '_' + str(tile.index) + '.npy')) + tile_data = np.load(file=tile_file) + slice_y = slice(tile.top_left_y, tile.bottom_right_y) + slice_x = slice(tile.top_left_x, tile.bottom_right_x) 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 + merged_array[slice_y, slice_x] = tile_data else: # 3D array - merged_array[t.top_left_y:t.bottom_right_y, t.top_left_x:t.bottom_right_x] = tile[:, :, index] + merged_array[slice_y, slice_x] = tile_data[:, :, index] - log.debug('Finished assembling tiles for {}'.format(out_type)) + log.debug(f'Finished assembling tiles for {out_type}') return merged_array @@ -282,8 +294,8 @@ def assemble_tiles(s: Tuple, dir: str, tiles: Tile, out_type: str, index: Option los_projection_out_types = {'tsincr', 'tscuml', 'linear_rate', 'stack_rate'} los_projection_divisors = { ifc.LINE_OF_SIGHT: lambda data: 1, - ifc.PSEUDO_VERTICAL: lambda data: np.cos(data), - ifc.PSEUDO_HORIZONTAL: lambda data: np.sin(data) + ifc.PSEUDO_VERTICAL: np.cos, + ifc.PSEUDO_HORIZONTAL: np.sin } @@ -291,10 +303,9 @@ def __save_merged_files(ifgs_dict, params, array, out_type, index=None, savenpy= """ Convenience function to save PyRate geotiff and numpy array files """ - outdir = params[C.OUT_DIR] ts_dir = params[C.TIMESERIES_DIR] vel_dir = params[C.VELOCITY_DIR] - log.debug('Saving PyRate outputs {}'.format(out_type)) + log.debug(f'Saving PyRate outputs {out_type}') gt, md, wkt = ifgs_dict['gt'], ifgs_dict['md'], ifgs_dict['wkt'] epochlist = ifgs_dict['epochlist'] @@ -312,18 +323,21 @@ def __save_merged_files(ifgs_dict, params, array, out_type, index=None, savenpy= md[ifc.DATA_TYPE] = out_type_md_dict[out_type] - if out_type in los_projection_out_types: # apply LOS projection and sign conversion for these outputs + # apply LOS projection and sign conversion for these outputs + los_proj = params[C.LOS_PROJECTION] + + if out_type in los_projection_out_types: incidence_path = Path(Configuration.geometry_files(params)['incidence_angle']) if incidence_path.exists(): # We can do LOS projection - if params[C.LOS_PROJECTION] == ifc.LINE_OF_SIGHT: + if los_proj == ifc.LINE_OF_SIGHT: log.info(f"Retaining Line-of-sight signal projection for file {dest}") - elif params[C.LOS_PROJECTION] != ifc.LINE_OF_SIGHT: + elif los_proj != ifc.LINE_OF_SIGHT: log.info(f"Projecting Line-of-sight signal into " - f"{ifc.LOS_PROJECTION_OPTION[params[C.LOS_PROJECTION]]} for file {dest}") + f"{ifc.LOS_PROJECTION_OPTION[los_proj]} for file {dest}") incidence = shared.Geometry(incidence_path) incidence.open() - array /= los_projection_divisors[params[C.LOS_PROJECTION]](incidence.data) * params[C.SIGNAL_POLARITY] - md[C.LOS_PROJECTION.upper()] = ifc.LOS_PROJECTION_OPTION[params[C.LOS_PROJECTION]] + array /= los_projection_divisors[los_proj](incidence.data) * params[C.SIGNAL_POLARITY] + md[C.LOS_PROJECTION.upper()] = ifc.LOS_PROJECTION_OPTION[los_proj] md[C.SIGNAL_POLARITY.upper()] = params[C.SIGNAL_POLARITY] if out_type in error_out_types: # add n-sigma value to error file metadata @@ -344,7 +358,7 @@ def __save_merged_files(ifgs_dict, params, array, out_type, index=None, savenpy= if savenpy: np.save(file=npy_file, arr=array) - log.debug('Finished saving {}'.format(out_type)) + log.debug(f'Finished saving {out_type}') def __merge_setup(params): @@ -353,7 +367,8 @@ def __merge_setup(params): """ # load previously saved preread_ifgs dict preread_ifgs_file = Configuration.preread_ifgs(params) - ifgs_dict = pickle.load(open(preread_ifgs_file, 'rb')) + with open(preread_ifgs_file, 'rb') as file: + ifgs_dict = pickle.load(file) ifgs = [v for v in ifgs_dict.values() if isinstance(v, shared.PrereadIfg)] shape = ifgs[0].shape tiles = Configuration.get_tiles(params) diff --git a/pyrate/prepifg.py b/pyrate/prepifg.py index 71dc1c5a9..603f0e86f 100644 --- a/pyrate/prepifg.py +++ b/pyrate/prepifg.py @@ -21,17 +21,18 @@ 1. A python/numpy version, which is also the default version, and 2. a `largetifs` option which can be activated using the config option largetifs: 1 -The python/numpy version is recommended when the both the input interferogram and multilooked interferogram can fit -into the memory allocated to the process. - -When dealing with relatively large (compared to memory available to the process) interferogram, the largetif option -can be used. This option uses a slightly modified version of the `gdal_calc.py` included with pyrate. Arbitrarily -large interferograms can be multilooked using largetifs. - -The `largetifs` option uses system calls to gdal utilities and avoid loading large chunks of memory. Our tests -(tests/test_prepifg_largetifs_vs_python.py::test_prepifg_largetifs_vs_python) indicate that for two different small -datasets included in our test suite, the `largetifs` option exactly match the output of the numpy version. However, -on large practical datasets we have observed numerical differences in the multilooked output in the 3rd and 4th decimal +The python/numpy version is recommended when the both the input interferogram and multilooked +interferogram can fit into the memory allocated to the process. + +When dealing with relatively large (compared to memory available to the process) interferogram, +the largetif option can be used. This option uses a slightly modified version of the `gdal_calc.py` +included with pyrate. Arbitrarily large interferograms can be multilooked using largetifs. + +The `largetifs` option uses system calls to gdal utilities and avoid loading large chunks of +memory. Our tests (tests/test_prepifg_largetifs_vs_python.py::test_prepifg_largetifs_vs_python) +indicate that for two different small datasets included in our test suite, the `largetifs` option +exactly match the output of the numpy version. However, on large practical datasets we have +observed numerical differences in the multilooked output in the 3rd and 4th decimal places (in sub mm range). """ @@ -51,8 +52,7 @@ from pyrate.core.logger import pyratelogger as log from pyrate.configuration import MultiplePaths, Configuration from pyrate.core.shared import Ifg, DEM -from pyrate.core.refpixel import convert_geographic_coordinate_to_pixel_value -from pyrate.core import ifgconstants as ifg +from pyrate.core.refpixel import convert_lat_lon_to_pixel_coord GAMMA = 1 @@ -125,7 +125,8 @@ def __calc_coherence_stats(params, ifg_path): phase_data = np.stack([i.phase_data for i in ifgs]) coh_stats = Configuration.coherence_stats(params) - for stat_func, out_type in zip([np.nanmedian, np.nanmean, np.nanstd], [ifg.COH_MEDIAN, ifg.COH_MEAN, ifg.COH_STD]): + stats = zip([np.nanmedian, np.nanmean, np.nanstd], [ifc.COH_MEDIAN, ifc.COH_MEAN, ifc.COH_STD]) + for stat_func, out_type in stats: with warnings.catch_warnings(): warnings.simplefilter("ignore", category=RuntimeWarning) arr = stat_func(phase_data, axis=0) @@ -134,7 +135,9 @@ def __calc_coherence_stats(params, ifg_path): __save_geom_files(ifg_path, dest, arr, out_type) -def do_prepifg(multi_paths: List[MultiplePaths], exts: Tuple[float, float, float, float], 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. @@ -146,8 +149,10 @@ def do_prepifg(multi_paths: List[MultiplePaths], exts: Tuple[float, float, float 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.") + raise FileNotFoundError( + f"Can not find geotiff: {f}. " + "Ensure you have converted your interferograms to geotiffs." + ) if params[C.LARGE_TIFS]: log.info("Using gdal system calls to execute 'prepifg' step") @@ -158,7 +163,9 @@ def do_prepifg(multi_paths: List[MultiplePaths], exts: Tuple[float, float, float res_str = ' '.join([str(e) for e in res_str]) if parallel: Parallel(n_jobs=params[C.PROCESSES], verbose=50)( - delayed(__prepifg_system)(exts, gtiff_path, params, res_str) for gtiff_path in multi_paths) + delayed(__prepifg_system)( + exts, gtiff_path, params, res_str) for gtiff_path in multi_paths + ) else: for m_path in multi_paths: __prepifg_system(exts, m_path, params, res_str) @@ -180,42 +187,48 @@ def do_prepifg(multi_paths: List[MultiplePaths], exts: Tuple[float, float, float def __prepifg_system(exts, gtiff, params, res): thresh = params[C.NO_DATA_AVERAGING_THRESHOLD] - p, c, l = _prepifg_multiprocessing(gtiff, exts, params) - log.info("Multilooking {p} into {l}".format(p=p, l=l)) + ifg_path, coh, sampled = _prepifg_multiprocessing(gtiff, exts, params) + log.info(f"Multilooking {ifg_path} into {sampled}") extents = ' '.join([str(e) for e in exts]) - - if isinstance(prepifg_helper.dem_or_ifg(p), shared.DEM): - check_call('gdalwarp {co} -te\t{extents}\t-tr\t{res}\t-r\taverage \t{p}\t{l}\n'.format( - co=COMMON_OPTIONS, extents=extents, res=res, p=p, l=l), shell=True) - __update_meta_data(p, c, l, params) + opts=COMMON_OPTIONS + + if isinstance(prepifg_helper.dem_or_ifg(ifg_path), shared.DEM): + check_call( + f'gdalwarp {opts} -te\t{extents}\t-tr\t{res}\t-r\taverage \t{ifg_path}\t{sampled}\n', + shell=True + ) + __update_meta_data(ifg_path, coh, sampled, params) return - p_unset = Path(params[C.OUT_DIR]).joinpath(Path(p).name).with_suffix('.unset.tif') - # change nodataval from zero, also leave input geotifs unchanged if one supplies conv2tif output/geotifs - check_call('gdal_translate {co} -a_nodata nan\t{p}\t{q}'.format(co=COMMON_OPTIONS, p=p, q=p_unset), - shell=True) + p_unset = Path(params[C.OUT_DIR]).joinpath(Path(ifg_path).name).with_suffix('.unset.tif') + # change nodataval from zero, also leave input geotifs unchanged if one supplies + # conv2tif output/geotifs + check_call(f'gdal_translate {opts} -a_nodata nan\t{ifg_path}\t{p_unset}', shell=True) # calculate nan-fraction # TODO: use output options and datatypes to reduce size of the next two tifs - nan_frac = Path(l).with_suffix('.nanfrac.tif') - nan_frac_avg = Path(l).with_suffix('.nanfrac.avg.tif') + nan_frac = Path(sampled).with_suffix('.nanfrac.tif') + nan_frac_avg = Path(sampled).with_suffix('.nanfrac.avg.tif') corrected_p = Path(p_unset).with_suffix('.corrected.tif') - if c is not None: + if coh is not None: # find all the nans - log.info(f"applying coherence + nodata masking on {p}") - check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} -A {p_unset} -B {c} --outfile={nan_frac}\t' - f'--calc=\"logical_or((B<{params[C.COH_THRESH]}), isclose(A,0,atol=0.000001))\"\t' - f'--NoDataValue=nan', shell=True) + log.info(f"applying coherence + nodata masking on {ifg_path}") + check_call( + f'{GDAL_CALC} {COMMON_OPTIONS2} -A {p_unset} -B {coh} --outfile={nan_frac}\t' + f'--calc=\"logical_or((B<{params[C.COH_THRESH]}), isclose(A,0,atol=0.000001))\"\t' + f'--NoDataValue=nan', + shell=True + ) # coh masking - check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} --overwrite -A {p_unset} -B {c}\t' + check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} --overwrite -A {p_unset} -B {coh}\t' f'--calc=\"A*(B>={params[C.COH_THRESH]}) - ' f'99999*logical_or((B<{params[C.COH_THRESH]}), isclose(A,0,atol=0.000001))\"\t' f'--outfile={corrected_p}\t' f'--NoDataValue=nan', shell=True) else: - log.info(f"applying nodata masking on {p}") + log.info(f"applying nodata masking on {ifg_path}") check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} --overwrite -A {p_unset}\t' f'--calc=\"isclose(A,0,atol=0.000001)\"\t' f'--outfile={nan_frac}\t' @@ -226,19 +239,23 @@ def __prepifg_system(exts, gtiff, params, res): f'--NoDataValue=nan', shell=True) # crop resample/average multilooking of nan-fraction - check_call('gdalwarp {co} -te\t{extents}\t-tr\t{res}\t-r\taverage\t{p}\t{out_file}'.format( - co=COMMON_OPTIONS, extents=extents, res=res, p=nan_frac, out_file=nan_frac_avg), shell=True) + check_call( + f'gdalwarp {opts} -te\t{extents}\t-tr\t{res}\t-r\taverage\t{nan_frac}\t{nan_frac_avg}', + shell=True + ) # crop resample/average multilooking of raster - check_call('gdalwarp {co} -te\t{extents}\t-tr\t{res}\t-r\taverage \t{p}\t{l}'.format( - co=COMMON_OPTIONS, extents=extents, res=res, p=corrected_p, l=l), shell=True) + check_call( + f'gdalwarp {opts} -te\t{extents}\t-tr\t{res}\t-r\taverage \t{corrected_p}\t{sampled}', + shell=True + ) - check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} --overwrite -A {nan_frac_avg}\t-B {l}\t' + check_call(f'{GDAL_CALC} {COMMON_OPTIONS2} --overwrite -A {nan_frac_avg}\t-B {sampled}\t' f'--calc=\"B*(A < {thresh}) -99999*(A >= {thresh})\"\t' - f'--outfile={l}\t' + f'--outfile={sampled}\t' f'--NoDataValue=nan', shell=True) - __update_meta_data(p_unset.as_posix(), c, l, params) + __update_meta_data(p_unset.as_posix(), coh, sampled, params) # clean up nan_frac_avg.unlink() @@ -247,7 +264,7 @@ def __prepifg_system(exts, gtiff, params, res): p_unset.unlink() -def __update_meta_data(p_unset, c, l, params): +def __update_meta_data(p_unset, coh, sampled, params): # update metadata ds = gdal.Open(p_unset) md = ds.GetMetadata() @@ -257,26 +274,32 @@ def __update_meta_data(p_unset, c, l, params): md[ifc.IFG_LKSY] = str(params[C.IFG_LKSY]) md[ifc.IFG_CROP] = str(params[C.IFG_CROP_OPT]) # update data type - if c is not None: # it's a interferogram when COH_MASK=1 - md_str = '-mo {k}={v}'.format(k=ifc.DATA_TYPE, v=ifc.MLOOKED_COH_MASKED_IFG) + if coh is not None: # it's a interferogram when COH_MASK=1 + md_str = f'-mo {ifc.DATA_TYPE}={ifc.MLOOKED_COH_MASKED_IFG}' md[ifc.COH_THRESH] = str(params[C.COH_THRESH]) else: if v == ifc.DEM: # it's a dem - md_str = '-mo {k}={v}'.format(k=ifc.DATA_TYPE, v=ifc.MLOOKED_DEM) + md_str = f'-mo {ifc.DATA_TYPE}={ifc.MLOOKED_DEM}' elif v == ifc.COH: - md_str = '-mo {k}={v}'.format(k=ifc.DATA_TYPE, v=ifc.MULTILOOKED_COH) + md_str = f'-mo {ifc.DATA_TYPE}={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(): - md_str += ' -mo {k}={v}'.format(k=k.replace(' ', '_'), v=v.replace(' ', '_')) - check_call('gdal_edit.py -unsetmd {md} {f}'.format(md=md_str, f=l), shell=True) + md_str = f'-mo {ifc.DATA_TYPE}={ifc.MULTILOOKED}' + for key, val in md.items(): + key = key.replace(' ', '_') + val = val.replace(' ', '_') + md_str += ' -mo {key}={val}' + check_call(f'gdal_edit.py -unsetmd {md_str} {sampled}', shell=True) ds = None # make prepifg output readonly - Path(l).chmod(0o444) # readonly output + Path(sampled).chmod(0o444) # readonly output -def _prepifg_multiprocessing(m_path: MultiplePaths, exts: Tuple[float, float, float, float], params: dict): +def _prepifg_multiprocessing( + m_path: MultiplePaths, + exts: Tuple[float, float, float, float], + params: dict +): """ Multiprocessing wrapper for prepifg """ @@ -298,20 +321,30 @@ def _prepifg_multiprocessing(m_path: MultiplePaths, exts: Tuple[float, float, fl if params[C.LARGE_TIFS]: return m_path.converted_path, coherence_path, m_path.sampled_path - else: - prepifg_helper.prepare_ifg(m_path.converted_path, xlooks, ylooks, exts, thresh, crop, - out_path=m_path.sampled_path, header=hdr, coherence_path=coherence_path, - coherence_thresh=coherence_thresh) - Path(m_path.sampled_path).chmod(0o444) # readonly output + + prepifg_helper.prepare_ifg( + m_path.converted_path, xlooks, ylooks, exts, thresh, crop, + out_path=m_path.sampled_path, + header=hdr, + coherence_path=coherence_path, + coherence_thresh=coherence_thresh + ) + Path(m_path.sampled_path).chmod(0o444) # readonly output + + # FIXME: I'm not sure how this has ever worked? (guessing params[C.LARGE_TIFS] is always true?) + # - we're supposed to return stuff here..? def find_header(path: MultiplePaths, params: dict) -> Dict[str, str]: + """ + Extract the header information out of an interferogram path, + attempting to automatically handle the differences between data file types + """ processor = params[C.PROCESSOR] # roipac, gamma or geotif tif_path = path.converted_path - if (processor == GAMMA) or (processor == GEOTIF): + if processor in (GAMMA, GEOTIF): header = gamma.gamma_header(tif_path, params) elif processor == ROIPAC: - import warnings warnings.warn("Warning: ROI_PAC support will be deprecated in a future PyRate release", category=DeprecationWarning) header = roipac.roipac_header(tif_path, params) @@ -348,22 +381,23 @@ def __write_geometry_files(params: dict, exts: Tuple[float, float, float, float] # using metadata of the first image in the stack # get pixel values of crop (needed to crop lookup table file) # pixel extent of cropped area (original IFG input) - xmin, ymax = convert_geographic_coordinate_to_pixel_value(exts[0], exts[1], transform) - xmax, ymin = convert_geographic_coordinate_to_pixel_value(exts[2], exts[3], transform) + xmin, ymax = convert_lat_lon_to_pixel_coord(exts[0], exts[1], transform) + xmax, ymin = convert_lat_lon_to_pixel_coord(exts[2], exts[3], transform) # xmin, xmax: columns of crop # ymin, ymax: rows of crop # calculate per-pixel radar coordinates - az, rg = geometry.calc_radar_coords(ifg, params, xmin, xmax, ymin, ymax) + az_loc, rg_loc = geometry.calc_radar_coords(ifg, params, xmin, xmax, ymin, ymax) # Read height data from DEM dem_file = params[C.DEM_FILE_PATH].sampled_path dem = DEM(dem_file) # calculate per-pixel look angle (also calculates and saves incidence and azimuth angles) - lk_ang, inc_ang, az_ang, rg_dist = geometry.calc_pixel_geometry(ifg, rg, lon.data, lat.data, dem.data) + geom = geometry.calc_pixel_geometry(ifg, rg_loc, lon.data, lat.data, dem.data) + lk_ang, inc_ang, az_ang, rg_dist = geom # save radar coordinates and angles to geotiff files - combinations = zip([az, rg, lk_ang, inc_ang, az_ang, rg_dist], C.GEOMETRY_OUTPUT_TYPES) + combinations = zip([az_loc, rg_loc, lk_ang, inc_ang, az_ang, rg_dist], C.GEOMETRY_OUTPUT_TYPES) shared.iterable_split(__parallelly_write, combinations, params, ifg_path) @@ -392,7 +426,7 @@ def __save_geom_files(ifg_path, dest, array, out_type): """ Convenience function to save geometry geotiff files """ - log.debug('Saving PyRate outputs {}'.format(out_type)) + log.debug(f'Saving PyRate outputs {out_type}') gt, md, wkt = shared.get_geotiff_header_info(ifg_path) md[ifc.DATA_TYPE] = out_type_md_dict[out_type] shared.remove_file_if_exists(dest) diff --git a/tests/test_geometry.py b/tests/test_geometry.py index fa5f5b35f..635caa2cb 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -7,7 +7,7 @@ import pyrate.constants as C from pyrate.core import ifgconstants as ifc from pyrate.core.geometry import get_lonlat_coords, get_sat_positions, vincinv -from pyrate.core.refpixel import convert_pixel_value_to_geographic_coordinate +from pyrate.core.refpixel import convert_pixel_coord_to_lat_lon from tests import common from pyrate.configuration import Configuration from subprocess import run, PIPE @@ -31,7 +31,7 @@ def get_lonlat_coords_slow(ifg: Ifg) -> Tuple[np.ndarray, np.ndarray]: lat = np.zeros((nrows, ncols)) # pre-allocate 2D numpy array for i in range(0, nrows): # rows are y-direction for j in range(0, ncols): # cols are x-direction - lon[i, j], lat[i, j] = convert_pixel_value_to_geographic_coordinate(j, i, transform) + lon[i, j], lat[i, j] = convert_pixel_coord_to_lat_lon(j, i, transform) return lon, lat diff --git a/tests/test_refpixel.py b/tests/test_refpixel.py index 432be278f..88f00b65a 100644 --- a/tests/test_refpixel.py +++ b/tests/test_refpixel.py @@ -30,7 +30,7 @@ import pyrate.constants as C import pyrate.core.refpixel from pyrate.core.refpixel import ref_pixel, _step, RefPixelError, ref_pixel_calc_wrapper, \ - convert_geographic_coordinate_to_pixel_value, convert_pixel_value_to_geographic_coordinate + convert_lat_lon_to_pixel_coord, convert_pixel_coord_to_lat_lon from pyrate.core import shared, ifgconstants as ifc from pyrate import correct, conv2tif, prepifg from pyrate.configuration import Configuration, ConfigException @@ -477,7 +477,7 @@ def x_y_pixel(): def test_convert_pixel_value_to_geographic_coordinate(x_y_pixel): transform = dem_transform() for x, y in x_y_pixel: - lon, lat = convert_pixel_value_to_geographic_coordinate(x, y, transform) + lon, lat = convert_pixel_coord_to_lat_lon(x, y, transform) out = run(f"gdallocationinfo -geoloc {SML_TEST_DEM_TIF} {lon} {lat}", shell=True, universal_newlines=True, stdout=PIPE).stdout xs = (x, x+1, x-1) @@ -495,6 +495,6 @@ def dem_transform(): def test_convert_geographic_coordinate_to_pixel_value(x_y_pixel): transform = dem_transform() for x, y in x_y_pixel: - lon, lat = convert_pixel_value_to_geographic_coordinate(x, y, transform) - xp, yp = convert_geographic_coordinate_to_pixel_value(lon, lat, transform) + lon, lat = convert_pixel_coord_to_lat_lon(x, y, transform) + xp, yp = convert_lat_lon_to_pixel_coord(lon, lat, transform) assert (xp == x) & (yp == y) diff --git a/tests/test_system.py b/tests/test_system.py index e7facb30a..2ab676d9e 100644 --- a/tests/test_system.py +++ b/tests/test_system.py @@ -50,8 +50,8 @@ def test_workflow(system_conf): shutil.rmtree(params[C.OUT_DIR]) -def test_single_workflow(gamma_or_mexicoa_conf): - if MPI_INSTALLED: +def test_single_workflow(pytestconfig, gamma_or_mexicoa_conf): + if MPI_INSTALLED and "not mpi" not in pytestconfig.getoption("-m"): check_call(f"mpirun -n 4 pyrate workflow -f {gamma_or_mexicoa_conf}", shell=True) else: check_call(f"pyrate workflow -f {gamma_or_mexicoa_conf}", shell=True) From 5e900bf4324ee08729509d9500dab026de1a2bce Mon Sep 17 00:00:00 2001 From: Mitchell Wheeler Date: Tue, 31 May 2022 17:00:42 +1000 Subject: [PATCH 2/5] Removed a lot of disabled pylint message types (which had no really valid reason for being disabled), kept a few disabled with justification for why... and tweaked some other parameters to be more suitable for algorithmic/math heavy code in PyRate --- .pylintrc | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/.pylintrc b/.pylintrc index 51ef01791..866a9de7b 100644 --- a/.pylintrc +++ b/.pylintrc @@ -65,7 +65,14 @@ confidence= # --enable=similarities". If you want to run only the classes checker, but have # no Warning level messages displayed, use"--disable=all --enable=classes # --disable=W" -disable=logging-format-interpolation,too-few-public-methods,import-star-module-level,old-octal-literal,oct-method,print-statement,unpacking-in-except,parameter-unpacking,backtick,old-raise-syntax,old-ne-operator,long-suffix,dict-view-method,dict-iter-method,metaclass-assignment,next-method-called,raising-string,indexing-exception,raw_input-builtin,long-builtin,file-builtin,execfile-builtin,coerce-builtin,cmp-builtin,buffer-builtin,basestring-builtin,apply-builtin,filter-builtin-not-iterating,using-cmp-argument,useless-suppression,range-builtin-not-iterating,suppressed-message,no-absolute-import,old-division,cmp-method,reload-builtin,zip-builtin-not-iterating,intern-builtin,unichr-builtin,reduce-builtin,standarderror-builtin,unicode-builtin,xrange-builtin,coerce-method,delslice-method,getslice-method,setslice-method,input-builtin,round-builtin,hex-method,nonzero-method,map-builtin-not-iterating,fixme +disable=missing-module-docstring,logging-fstring-interpolation,import-star-module-level,too-few-public-methods,backtick,fixme + +# JUSTIFICATION: +# missing-module-docstring -> many of our modules are self-explanitory, they hold a single class or set of functions which are well-documented already, the few modules which are more like 'scripts' are indeed documented in multiple ways as well. +# logging-fstring-interpolation -> this is a relic of the past, f-string everywhere is more consisten and more performant. +# too-few-public-methods -> there's a few data classes and the likes that exist currently which aren't promoted to real data classes (needs to be fixed eventually) +# backtick -> this is used in places it genuinely looks cleaner to split code over multiple lines. +# fixme -> these are used as reminders w/ links to GH issues for future work that needs to happen. [REPORTS] @@ -226,7 +233,7 @@ expected-line-ending-format= [BASIC] # Good variable names which should always be accepted, separated by a comma -good-names=ts,i,j,k,ex,Run,_,x,w,f,n,v,t,r,d,c,cd,g,r,q,V,m,A,b,log,y,id,ds,gt,md +good-names=ts,i,j,k,ex,Run,_,x,w,f,n,v,t,r,d,c,cd,g,r,q,V,m,A,b,log,y,id,ds,gt,md,ax,im # Bad variable names which should always be refused, separated by a comma bad-names=foo,bar,baz,toto,tutu,tata @@ -243,10 +250,10 @@ include-naming-hint=no property-classes=abc.abstractproperty # Regular expression matching correct function names -function-rgx=[a-z_][a-z0-9_]{2,30}$ +function-rgx=[a-z_][a-z0-9_]{2,40}$ # Naming hint for function names -function-name-hint=[a-z_][a-z0-9_]{2,30}$ +function-name-hint=[a-z_][a-z0-9_]{2,40}$ # Regular expression matching correct variable names variable-rgx=[a-z_][a-z0-9_]{2,30}$ @@ -320,14 +327,16 @@ max-nested-blocks=5 [DESIGN] # Maximum number of arguments for function / method -max-args=5 +# - PyRate allows up to 8 (some functions have a few settings that don't make sense to make data-classes for) +max-args=8 # Argument names that match this expression will be ignored. Default to name # with leading underscore ignored-argument-names=_.* # Maximum number of locals for function / method body -max-locals=15 +# - PyRate has some numerically complicated functions that don't make sense to split up +max-locals=20 # Maximum number of return / yield for function / method body max-returns=6 From 779702ba64292d8372c09f27e4c18ad7c3fc2c47 Mon Sep 17 00:00:00 2001 From: Mitchell Wheeler Date: Wed, 1 Jun 2022 13:31:37 +1000 Subject: [PATCH 3/5] Fixed a small logging bug I introduced in a refactor (this is not the cause of the unit test failure, which I still haven't found the source of) --- pyrate/core/shared.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyrate/core/shared.py b/pyrate/core/shared.py index f36ff9cef..d31f57941 100644 --- a/pyrate/core/shared.py +++ b/pyrate/core/shared.py @@ -1376,11 +1376,11 @@ def close_all(ifgs): if not all(flags) and any(flags): log.debug('Detected mix of corrected and uncorrected interferograms') - for flag in flags: + for i, flag in zip(ifgs, flags): if flag: - msg = '{i.data_path}: correction detected' + msg = f'{i.data_path}: correction detected' else: - msg = '{i.data_path}: correction NOT detected' + msg = f'{i.data_path}: correction NOT detected' log.debug(msg) close_all(ifgs) raise CorrectionStatusError(msg) From 8fe921e123955d41532f4ef6ea9b5a9f177daaaf Mon Sep 17 00:00:00 2001 From: Mitchell Wheeler Date: Wed, 1 Jun 2022 14:31:18 +1000 Subject: [PATCH 4/5] Testing a possible cause of MPI failures (for lack of anything else I can find that could impact how anything runs) --- pyrate/configuration.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyrate/configuration.py b/pyrate/configuration.py index 3db878210..1aa9e919e 100644 --- a/pyrate/configuration.py +++ b/pyrate/configuration.py @@ -226,11 +226,11 @@ def __init__(self, config_file_path: Union[str, Path]): config_file_path = Path(config_file_path) # Setup default values (in case they're not in the config file) - self.parallel = False - self.processes = 1 - self.cohfilelist = None - self.basefilelist = None - self.demfile = None + #self.parallel = False + #self.processes = 1 + #self.cohfilelist = None + #self.basefilelist = None + #self.demfile = None # Load config file parser = ConfigParser() From 6252e3d17ab28124fb8229e631690dab872233d3 Mon Sep 17 00:00:00 2001 From: Mitchell Wheeler Date: Wed, 1 Jun 2022 15:54:08 +1000 Subject: [PATCH 5/5] Un-did last checking which incorrectly attempted to solve an issue in Configuration, and fixed the issue that has been identified a missing f-string in my prepifg.py refactoring --- pyrate/configuration.py | 10 +++++----- pyrate/prepifg.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/pyrate/configuration.py b/pyrate/configuration.py index 1aa9e919e..3db878210 100644 --- a/pyrate/configuration.py +++ b/pyrate/configuration.py @@ -226,11 +226,11 @@ def __init__(self, config_file_path: Union[str, Path]): config_file_path = Path(config_file_path) # Setup default values (in case they're not in the config file) - #self.parallel = False - #self.processes = 1 - #self.cohfilelist = None - #self.basefilelist = None - #self.demfile = None + self.parallel = False + self.processes = 1 + self.cohfilelist = None + self.basefilelist = None + self.demfile = None # Load config file parser = ConfigParser() diff --git a/pyrate/prepifg.py b/pyrate/prepifg.py index 603f0e86f..9e660214e 100644 --- a/pyrate/prepifg.py +++ b/pyrate/prepifg.py @@ -287,7 +287,7 @@ def __update_meta_data(p_unset, coh, sampled, params): for key, val in md.items(): key = key.replace(' ', '_') val = val.replace(' ', '_') - md_str += ' -mo {key}={val}' + md_str += f' -mo {key}={val}' check_call(f'gdal_edit.py -unsetmd {md_str} {sampled}', shell=True) ds = None