shared.py

#   This Python module is part of the PyRate software package.
#
#   Copyright 2020 Geoscience Australia
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.
# coding: utf-8
"""
This Python module contains utilities and classes shared by
all other PyRate modules
"""
# pylint: disable=too-many-lines
import re
from typing import List, Union

import errno
import math
from joblib import Parallel, delayed
from math import floor
import os
from os.path import basename, join
from pathlib import Path
import struct
from datetime import date
from itertools import product
from enum import Enum
import numpy as np
from numpy import where, nan, isnan, sum as nsum, isclose
import pyproj
import pkg_resources
try:
    from osgeo import osr, gdal
    from osgeo.gdalconst import GA_Update, GA_ReadOnly
except ImportError:
    import gdal

from pyrate.core import ifgconstants as ifc, mpiops, config as cf
from pyrate.core.logger import pyratelogger as log


gdal.UseExceptions()

VERBOSE = True

# Constants
PHASE_BAND = 1
RADIANS = 'RADIANS'
MILLIMETRES = 'MILLIMETRES'
GAMMA = 'GAMMA'
ROIPAC = 'ROIPAC'

# GDAL projection list
GDAL_X_CELLSIZE = 1
GDAL_Y_CELLSIZE = 5
GDAL_X_FIRST = 0
GDAL_Y_FIRST = 3


class InputTypes(Enum):
    IFG = 'ifg'
    COH = 'coh'
    DEM = 'dem'
    HEADER = 'header'


def joblib_log_level(level: str) -> int:
    """
    Convert python log level to joblib int verbosity.
    """ 
    if level == 'INFO':
        return 0
    else:
        return 60

def mkdir_p(path):
    """
    Make new directory and create parent directories as necessary.
    Copied from stackoverflow:
    http://stackoverflow.com/questions/600268/mkdir-p-functionality-in-python

    :param str path: Path name for new directory
    """
    try:
        os.makedirs(path)
    except OSError as exc:  # Python >2.5
        if exc.errno == errno.EEXIST and os.path.isdir(path):
            pass
        else:
            raise


class RasterBase(object):
    """
    Base class for PyRate GeoTIFF based raster datasets.
    """
    # pylint: disable=missing-docstring
    # pylint: disable=too-many-instance-attributes
    def __init__(self, path: Union[gdal.Dataset, str]):
        if isinstance(path, gdal.Dataset):
            self.dataset = path  # path will be Dataset in this case
            self.data_path = self.dataset  # data_path dummy
            self.add_geographic_data()
        else:
            self.data_path = path
            self.dataset = None  # for GDAL dataset obj
            self._readonly = not os.access(path, os.R_OK | os.W_OK)

            if self._readonly is None:
                raise NotImplementedError  # os.access() has failed?

    def __str__(self):
        name = self.__class__.__name__
        return "%s('%s')" % (name, self.data_path)

    def __repr__(self):
        name = self.__class__.__name__
        return "%s('%s')" % (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
            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))

        # unless read only, by default open files as writeable
        if readonly not in [True, False, None]:
            raise ValueError("readonly must be True, False or None")

        if readonly is False and self._readonly is True:
            raise IOError("Cannot open write protected file for writing")

        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)

        self.add_geographic_data()

    def add_geographic_data(self):
        """
        Determine and add geographic data to object
        """
        # add some geographic data
        self.x_centre = int(self.ncols / 2)
        self.y_centre = int(self.nrows / 2)
        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)

    @property
    def ncols(self):
        """
        Number of raster columns
        """
        return self.dataset.RasterXSize

    @property
    def nrows(self):
        """
        Number of raster rows
        """
        return self.dataset.RasterYSize

    @property
    def x_step(self):
        """
        Raster pixel size in X (easting) dimension
        """
        return float(self.dataset.GetGeoTransform()[GDAL_X_CELLSIZE])

    @property
    def y_step(self):
        """
        Raster pixel size in Y (northing) dimension
        """
        return float(self.dataset.GetGeoTransform()[GDAL_Y_CELLSIZE])

    @property
    def x_first(self):
        """
        Raster western bounding coordinate
        """
        return float(self.dataset.GetGeoTransform()[GDAL_X_FIRST])

    @property
    def x_last(self):
        """
        Raster eastern bounding coordinate
        """
        return self.x_first + (self.x_step * self.ncols)

    @property
    def y_first(self):
        """
        Raster northern bounding coordinate
        """
        return float(self.dataset.GetGeoTransform()[GDAL_Y_FIRST])

    @property
    def y_last(self):
        """
        Raster southern bounding coordinate
        """
        return self.y_first + (self.y_step * self.nrows)

    @property
    def shape(self):
        """
        Returns tuple of (Y,X) shape of the raster (as per numpy.shape).
        """
        return self.dataset.RasterYSize, self.dataset.RasterXSize

    @property
    def num_cells(self):
        """
        Total number of pixels in raster dataset
        """
        if self.is_open:
            return self.dataset.RasterXSize * self.dataset.RasterYSize
        else:
            raise RasterException('Dataset not open')

    @property
    def is_open(self):
        """
        Returns True if the underlying dataset has been opened by GDAL.
        """
        return self.dataset is not None

    def close(self):
        """
        Explicitly closes file opened by gdal.Open()
        This is required in windows, otherwise opened files can not be removed,
        because windows locks open files.
        """
        if self.is_open:
            self.dataset = None

    @property
    def is_read_only(self):
        """
        Determines file permissions
        """
        return self._readonly

    def _get_band(self, band):
        """
        Wrapper (with error checking) for GDAL's Band.GetRasterBand() method.

        :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)


class Ifg(RasterBase):
    """
    Interferogram (Ifg) class objects; double as a container for
    interferometric phase raster band data and related data.
    """
    # pylint: disable=too-many-instance-attributes
    def __init__(self, path: Union[str, Path]):
        """
        Interferogram constructor, for 2-band Ifg raster datasets.

        :param str path: Path to interferogram file
        """
        path = path.as_posix() if isinstance(path, Path) else path
        RasterBase.__init__(self, path)
        self._phase_band = None
        self._phase_data = None
        self.first = None
        self.second = None
        self.nan_converted = False
        self.mm_converted = False
        self.meta_data = None
        self.wavelength = None
        self._nodata_value = None
        self.time_span = None

    def open(self, readonly=None):
        """
        Open raster file.

        :param bool readonly: True/False, or None to open as underlying file setting
        """
        RasterBase.open(self, readonly)
        self.initialize()

    def initialize(self):
        """
        Read basic interferogram properties upon opening interferogram.
        """
        self._init_dates()
        md = self.dataset.GetMetadata()
        self.wavelength = float(md[ifc.PYRATE_WAVELENGTH_METRES])
        self.meta_data = md
        self.nan_converted = False # This flag set True after NaN conversion

    def _init_dates(self):
        """
        Determine first and second image dates, and interferogram timespan
        """
        def _to_date(datestr):
            year, month, day = [int(i) for i in datestr.split('-')]
            return date(year, month, day)

        md = self.dataset.GetMetadata()
        datestrs = [md[k] for k in [ifc.FIRST_DATE, ifc.SECOND_DATE]]

        if all(datestrs):
            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
            raise IfgException(msg)

    def convert_to_nans(self):
        """
        Convert phase data of given value to NaN
        """
        if (self._nodata_value is None) \
                or (self.dataset is None):  # pragma: no cover
            msg = 'nodata value needs to be set for nan conversion.' \
                  '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)
            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.write_modified_phase(self.phase_data)

    @property
    def phase_band(self):
        """
        Returns a GDAL Band object for the phase band.
        """
        if self._phase_band is None:
            self._phase_band = self._get_band(PHASE_BAND)
        return self._phase_band

    @property
    def nodata_value(self):
        """
        Determine the no-data value in phase band.
        """
        return self._nodata_value

    @nodata_value.setter
    def nodata_value(self, val):
        """
        Set the no-data value for phase band.
        """
        self._nodata_value = val

    @property
    def phase_data(self):
        """
        Returns phase band as an array.
        """
        if self._phase_data is None:
            self._phase_data = self.phase_band.ReadAsArray()
        return self._phase_data

    def convert_to_mm(self):
        """
        Convert phase data units from radians to millimetres.
        """
        self.mm_converted = True
        if self.dataset.GetMetadataItem(ifc.DATA_UNITS) == MILLIMETRES:
            msg = '{}: ignored as previous phase unit conversion ' \
                  'already applied'.format(self.data_path)
            log.debug(msg)
            self.phase_data = self.phase_data
            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
            # self.write_modified_phase()
            # 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)
            log.debug(msg)
        else:  # pragma: no cover
            msg = 'Phase units are not millimetres or radians'
            raise IfgException(msg)

        # self.write_modified_phase(self.phase_data)

    @phase_data.setter
    def phase_data(self, data):
        """
        Set phase data value
        """
        self._phase_data = data

    @property
    def phase_rows(self):
        """
        Generator returning each row of the phase data.
        """
        for y in range(self.nrows):
            r = self.phase_band.ReadAsArray(yoff=y,
                                            win_xsize=self.ncols, win_ysize=1)
            yield r[0] # squeezes row from (1, WIDTH) to 1D array

    @property
    def nan_count(self):
        """
        Returns total number of NaN cells in the phase data.
        """
        return nsum(isnan(self.phase_data))

    @property
    def nan_fraction(self):
        """
        Returns decimal fraction of NaN cells in the phase band.
        """
        if (self._nodata_value is None) or (self.dataset is None):
            msg = 'nodata_value needs to be set for nan fraction calc.' \
                  'Use ifg.nondata = NoDataValue to set nodata'
            raise RasterException(msg)
        # don't cache nan_count as client code may modify phase data
        nan_count = self.nan_count
        # handle datasets with no 0 -> NaN replacement
        if not self.nan_converted and (nan_count == 0):
            nan_count = nsum(np.isclose(self.phase_data,
                                        self._nodata_value, atol=1e-6))
        return nan_count / float(self.num_cells)

    def write_modified_phase(self, data=None):
        """
        Write updated phase data to file on disk.
        """
        if self.is_read_only:
            raise IOError("Cannot write to read only Ifg")

        # keep this block
        # if new_data_path is None:
        #     self.dataset = gdal.Open(self.data_path, GA_Update)
        # else:
        #     self.dataset = gdal.Open(new_data_path, GA_Update)
        # self._phase_band = None

        if data is not None:
            assert isinstance(data, np.ndarray)
            data_r, data_c = data.shape
            assert data_r == self.nrows and data_c == self.ncols
            self.phase_data = data
        self.phase_band.WriteArray(self.phase_data)
        for k, v in self.meta_data.items():
            self.dataset.SetMetadataItem(k, v)
        self.dataset.FlushCache()

    def add_metadata(self, **kwargs):
        if (not self.is_open) or self.is_read_only:
            raise IOError("Ifg not open or readonly. Cannot write!")

        for k, v in kwargs.items():
            self.dataset.SetMetadataItem(k, v)
        self.dataset.FlushCache()  # write to disc


class IfgPart(object):
    """
    Create a tile (subset) of an Ifg data object
    """
    def __init__(self, ifg_or_path, tile, ifg_dict=None, params=None):
        """
        Interferogram tile constructor.

        :param str path: Path to interferogram file
        """
        self.tile = tile
        self.r_start = self.tile.top_left_y
        self.r_end = self.tile.bottom_right_y
        self.c_start = self.tile.top_left_x
        self.c_end = self.tile.bottom_right_x
        if ifg_dict is not None:  # should be used with MPI
            ifg = ifg_dict[ifg_or_path]
            self.nan_fraction = ifg.nan_fraction
            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)
            self.phase_data = np.load(join(params[cf.TMPDIR], phase_file))
        else:
            # check if Ifg was sent.
            if isinstance(ifg_or_path, Ifg):
                ifg = ifg_or_path
            else:
                self.data_path = ifg_or_path
                ifg = Ifg(ifg_or_path)
            self.phase_data = None
            self.nan_fraction = None
            self.first = None
            self.second = None
            self.time_span = None
        if isinstance(ifg, Ifg):
            self.read_required(ifg)

    def read_required(self, ifg):
        """
        Read interferogram file if not already open.
        """
        if not ifg.is_open:
            ifg.open(readonly=True)
        ifg.nodata_value = 0
        self.phase_data = ifg.phase_data[self.r_start:self.r_end,
                                         self.c_start:self.c_end]
        self.nan_fraction = ifg.nan_fraction
        self.first = ifg.first
        self.second = ifg.second
        self.time_span = ifg.time_span
        ifg.phase_data = None
        ifg.close()  # close base ifg

    @property
    def nrows(self):
        """
        Determine number of rows in tile.
        """
        return self.r_end - self.r_start

    @property
    def ncols(self):
        """
        Determine number of columns in tile.
        """
        return self.c_end - self.c_start


class Incidence(RasterBase):   # pragma: no cover
    """
    Class for storing viewing geometry data.
    e.g. incidence and azimuth raster values
    """

    def __init__(self, path):
        """
        Incidence object constructor.
        """
        RasterBase.__init__(self, path)
        self._incidence_band = None
        self._azimuth_band = None
        self._incidence_data = None
        self._azimuth_data = None

    @property
    def incidence_band(self):
        """
        Returns the GDALBand for the incidence angle layer.
        """

        if self._incidence_band is None:
            self._incidence_band = self._get_band(1)
        return self._incidence_band

    @property
    def incidence_data(self):
        """
        Returns the incidence band as an array.
        """
        if self._incidence_data is None:
            self._incidence_data = self.incidence_band.ReadAsArray()
        return self._incidence_data

    @property
    def azimuth_band(self):
        """
        Returns the GDALBand for the azimuth layer.
        """
        if self._azimuth_band is None:
            self._azimuth_band = self._get_band(2)
        return self._azimuth_band

    @property
    def azimuth_data(self):
        """
        Returns the azimuth band as an array.
        """
        if self._azimuth_data is None:
            self._azimuth_data = self.azimuth_band.ReadAsArray()
        return self._azimuth_data


class DEM(RasterBase):
    """
    Generic raster class for single band DEM files.
    """

    def __init__(self, path):
        """
        DEM constructor.
        """
        RasterBase.__init__(self, path)
        self._band = None

    @property
    def height_band(self):
        """
        Returns the GDALBand for the elevation layer.
        """
        if self._band is None:
            self._band = self._get_band(1)
        return self._band


class IfgException(Exception):
    """
    Generic exception class for interferogram errors.
    """


class RasterException(Exception):
    """
    Generic exception for raster errors.
    """


class EpochList(object):
    """
    Metadata container for epoch related information.
    """

    def __init__(self, dates=None, repeat=None, spans=None):
        """
        Construct epochlist object
        """
        self.dates = dates # list of unique dates from all the ifgs
        self.repeat = repeat
        self.spans = spans  # time span from earliest ifg

    def __str__(self):
        return "EpochList: %s" % str(self.dates)

    def __repr__(self):
        return "EpochList: %s" % repr(self.dates)


def convert_radians_to_mm(data, wavelength):
    """
    Function to translates phase in units of radians to units in millimetres.

    :param ndarray data: Interferogram phase data array
    :param float wavelength: Radar wavelength in metres

    :return: data: converted phase data
    :rtype: ndarray
    """
    return data * ifc.MM_PER_METRE * (wavelength / (4 * math.pi))


def nanmedian(x):
    """
    Determine the median of values excluding nan values.
    Use different numpy algorithm dependent on numpy version.

    :param ndarray x: array of numeric data.

    :return: y: median value
    :rtype: float
    """
    # pylint: disable=no-member
    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)])


def _is_interferogram(hdr):
    """
    Convenience function to determine if file is interferogram
    """
    return (ifc.PYRATE_WAVELENGTH_METRES in hdr) and \
           (hdr[ifc.INPUT_TYPE] == InputTypes.IFG if ifc.INPUT_TYPE in hdr else True)


def _is_coherence(hdr):
    """
    Convenience function to determine if file is interferogram
    """
    return (ifc.PYRATE_WAVELENGTH_METRES in hdr) and \
           (hdr[ifc.INPUT_TYPE] == InputTypes.COH if ifc.INPUT_TYPE in hdr else False)


def _is_incidence(hdr):
    """
    Convenience function to determine if incidence file
    """
    return 'FILE_TYPE' in hdr


def write_fullres_geotiff(header, data_path, dest, nodata):
    # pylint: disable=too-many-statements
    """
    Creates a copy of input image data (interferograms, DEM, incidence maps
    etc) in GeoTIFF format with PyRate metadata.

    :param dict header: Interferogram metadata dictionary
    :param str data_path: Input file
    :param str dest: Output destination file
    :param float nodata: No-data value

    :return: None, file saved to disk
    """
    # pylint: disable=too-many-branches
    # pylint: disable=too-many-locals
    ifg_proc = header[ifc.PYRATE_INSAR_PROCESSOR]
    ncols = header[ifc.PYRATE_NCOLS]
    nrows = header[ifc.PYRATE_NROWS]
    bytes_per_col, fmtstr = _data_format(ifg_proc, _is_interferogram(header), ncols)
    if _is_interferogram(header) and ifg_proc == ROIPAC:
        # roipac ifg has 2 bands
        _check_raw_data(bytes_per_col*2, data_path, ncols, nrows)
    else:
        _check_raw_data(bytes_per_col, data_path, ncols, nrows)

    _check_pixel_res_mismatch(header)

    # 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]]
    srs = osr.SpatialReference()
    res = srs.SetWellKnownGeogCS(header[ifc.PYRATE_DATUM])
    if res:
        msg = 'Unrecognised projection: %s' % header[ifc.PYRATE_DATUM]
        raise GeotiffException(msg)

    wkt = srs.ExportToWkt()
    dtype = 'float32' if (_is_interferogram(header) or _is_incidence(header) or _is_coherence(header)) else 'int16'

    # get subset of metadata relevant to PyRate
    md = collate_metadata(header)

    # create GDAL object
    ds = gdal_dataset(
        dest, ncols, nrows, driver="GTiff", bands=1, dtype=dtype, metadata=md, crs=wkt, geotransform=gt,
        creation_opts=["compress=packbits"]
    )

    # copy data from the binary file
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(nodata)

    row_bytes = ncols * bytes_per_col

    with open(data_path, 'rb') as f:
        for y in range(nrows):
            if ifg_proc == ROIPAC:
                if _is_interferogram(header):
                    f.seek(row_bytes, 1)  # skip interleaved band 1

            data = struct.unpack(fmtstr, f.read(row_bytes))
            band.WriteArray(np.array(data).reshape(1, ncols), yoff=y)

    ds = None  # manual close
    del ds


def gdal_dataset(out_fname, columns, rows, driver="GTiff", bands=1,
                 dtype='float32', metadata=None, crs=None,
                 geotransform=None, creation_opts=None):
    """
    Initialises a py-GDAL dataset object for writing image data.
    """
    if dtype == 'float32':
        gdal_dtype = gdal.GDT_Float32
    elif dtype == 'int16':
        gdal_dtype = gdal.GDT_Int16
    else:
        # assume gdal.GDT val is passed to function
        gdal_dtype = dtype

    # create output dataset
    driver = gdal.GetDriverByName(driver)
    outds = driver.Create(out_fname, columns, rows, bands, gdal_dtype, options=creation_opts)

    # geospatial info
    outds.SetGeoTransform(geotransform)
    outds.SetProjection(crs)

    # add metadata
    if metadata is not None:
        for k, v in metadata.items():
            outds.SetMetadataItem(k, str(v))

    return outds


def collate_metadata(header):
    """
    Grab metadata relevant to PyRate from input metadata

    :param dict header: Input file metadata dictionary

    :return: dict of relevant metadata for PyRate
    """
    md = dict()

    def __common_ifg_coh_update(header, md):
        for k in [ifc.PYRATE_WAVELENGTH_METRES, ifc.PYRATE_TIME_SPAN,
                  ifc.PYRATE_INSAR_PROCESSOR,
                  ifc.FIRST_DATE, ifc.SECOND_DATE,
                  ifc.DATA_UNITS]:
            md.update({k: str(header[k])})
        if header[ifc.PYRATE_INSAR_PROCESSOR] == GAMMA:
            for k in [ifc.FIRST_TIME, ifc.SECOND_TIME, ifc.PYRATE_INCIDENCE_DEGREES]:
                md.update({k: str(header[k])})

    if _is_coherence(header):
        __common_ifg_coh_update(header, md)
        md.update({ifc.DATA_TYPE: ifc.COH})
    elif _is_interferogram(header):
        __common_ifg_coh_update(header, md)
        md.update({ifc.DATA_TYPE: ifc.ORIG})
    elif _is_incidence(header):
        md.update({ifc.DATA_TYPE: ifc.INCIDENCE})
    else:  # must be dem
        md.update({ifc.DATA_TYPE: ifc.DEM})

    return md


def _data_format(ifg_proc, is_ifg, ncols):
    """
    Convenience function to determine the bytesize and format of input files
    """
    if ifg_proc == GAMMA:
        fmtstr = '!' + ('f' * ncols)  # data format is big endian float32s
        bytes_per_col = 4
    elif ifg_proc == ROIPAC:
        if is_ifg:
            fmtstr = '<' + ('f' * ncols)  # roipac ifgs are little endian float32s
            bytes_per_col = 4
        else:
            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
        raise GeotiffException(msg)
    return bytes_per_col, fmtstr


def _check_raw_data(bytes_per_col, data_path, ncols, nrows):
    """
    Convenience function to check the file size is as expected
    """
    size = ncols * nrows * bytes_per_col
    act_size = os.stat(data_path).st_size
    if act_size != size:
        msg = '%s should have size %s, not %s. Is the correct file being used?'
        raise GeotiffException(msg % (data_path, size, act_size))


def _check_pixel_res_mismatch(header):
    """
    Convenience function to check equality of pixel resolution in X and Y dimensions
    """
    # pylint: disable=invalid-name
    xs, ys = [abs(i) for i in [header[ifc.PYRATE_X_STEP], header[ifc.PYRATE_Y_STEP]]]

    if xs != ys:
        msg = 'X and Y cell sizes do not match: %s & %s'
        raise GeotiffException(msg % (xs, ys))


def write_unw_from_data_or_geotiff(geotif_or_data, dest_unw, ifg_proc):
    """
    Function to write numpy array data or a geotiff to a GAMMA-format
    big-endian float32 .unw file

    :param str geotif_or_data: path name of geotiff file to convert
        OR
    :param ndarray geotif_or_data: numpy array of data to convert
    :param str dest_unw: destination unw file
    :param int ifg_proc: processor type, GAMMA=1, ROIPAC=0

    :return: None, file saved to disk
    """
    if ifg_proc != 1:
        raise NotImplementedError('only supports GAMMA format for now')
    if isinstance(geotif_or_data, str):
        assert os.path.exists(geotif_or_data), 'make sure geotif exists'
        ds = gdal.Open(geotif_or_data)
        data = ds.ReadAsArray()
        ds = None
    else:
        data = geotif_or_data

    nrows, ncols = data.shape
    fmtstr = '!' + ('f' * ncols)  # data format is big endian float32s

    with open(dest_unw, 'wb') as f:
        for y in range(nrows):
            col_data = struct.pack(fmtstr, *data[y, :])
            f.write(col_data)


# This function may be able to be deprecated
def write_output_geotiff(md, gt, wkt, data, dest, nodata):
    # pylint: disable=too-many-arguments
    """
    Writes PyRate output data to a GeoTIFF file.

    :param dict md: Dictionary containing PyRate metadata
    :param list gt: GDAL geotransform for the data
    :param list wkt: GDAL projection information for the data
    :param ndarray data: Output data array to save
    :param str dest: Destination file name
    :param float nodata: No data value of data

    :return None, file saved to disk
    """

    driver = gdal.GetDriverByName("GTiff")
    nrows, ncols = data.shape
    ds = driver.Create(dest, ncols, nrows, 1, gdal.GDT_Float32, options=['compress=packbits'])
    # set spatial reference for geotiff
    ds.SetGeoTransform(gt)
    ds.SetProjection(wkt)
    ds.SetMetadataItem(ifc.EPOCH_DATE, str(md[ifc.EPOCH_DATE]))

    # set other metadata
    ds.SetMetadataItem('DATA_TYPE', str(md['DATA_TYPE']))

    # sequence position for time series products

    for k in [ifc.SEQUENCE_POSITION, ifc.PYRATE_REFPIX_X, ifc.PYRATE_REFPIX_Y, ifc.PYRATE_REFPIX_LAT,
              ifc.PYRATE_REFPIX_LON, ifc.PYRATE_MEAN_REF_AREA, ifc.PYRATE_STDDEV_REF_AREA]:
        if k in md:
            ds.SetMetadataItem(k, str(md[k]))

    # write data to geotiff
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(nodata)
    band.WriteArray(data, 0, 0)


def write_geotiff(data, outds, nodata):
    # pylint: disable=too-many-arguments
    """
    A generic routine for writing a NumPy array to a geotiff.

    :param ndarray data: Output data array to save
    :param obj outds: GDAL destination object
    :param float nodata: No data value of data

    :return None, file saved to disk
    """
    # only support "2 <= dims <= 3"
    if data.ndim == 3:
        count, height, width = data.shape
    elif data.ndim == 2:
        height, width = data.shape
    else:
        msg = "Only support dimensions of '2 <= dims <= 3'."
        raise GeotiffException(msg)

    # write data to geotiff
    band = outds.GetRasterBand(1)
    band.SetNoDataValue(nodata)
    band.WriteArray(data, 0, 0)

    outds = None
    band = None
    del outds
    del band


class GeotiffException(Exception):
    """
    Geotiff exception class
    """

def create_tiles(shape, nrows=2, ncols=2):
    """
    Return a list of tiles containing nrows x ncols with each tile preserving
    the physical layout of original array. The number of rows can be changed
    (increased) such that the resulting tiles with float32's do not exceed
    500MB in memory. When the array shape (rows, columns) are not divisible
    by (nrows, ncols) then some of the array dimensions can change according
    to numpy.array_split.

    :param tuple shape: Shape tuple (2-element) of interferogram.
    :param int nrows: Number of rows of tiles
    :param int ncols: Number of columns of tiles

    :return: List of Tile class instances.
    :rtype: list
    """

    if len(shape) != 2:
        raise ValueError('shape must be a length 2 tuple')

    no_y, no_x = shape

    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))]


class Tile():
    """
    Tile class for containing a sub-part of an interferogram
    """
    def __init__(self, index, top_left, bottom_right):
        """
        Parameters
        ----------
        index: int
            identifying index of a tile
        top_left: tuple
            ifg index of top left of tile
        bottom_right: tuple
            ifg index of bottom right of tile
        """

        self.index = index
        self.top_left = top_left
        self.bottom_right = bottom_right
        self.top_left_y, self.top_left_x = top_left
        self.bottom_right_y, self.bottom_right_x = bottom_right

    def __str__(self):
        return "Convenience Tile class containing tile co-ordinates"


def get_tiles(ifg_path, rows, cols) -> List[Tile]:
    """
    Break up the interferograms into smaller tiles based on user supplied
    rows and columns.

    :param list ifg_path: List of destination geotiff file names
    :param int rows: Number of rows to break each interferogram into
    :param int cols: Number of columns to break each interferogram into

    :return: tiles: List of shared.Tile instances
    :rtype: list
    """
    ifg = Ifg(ifg_path)
    ifg.open(readonly=True)
    tiles = create_tiles(ifg.shape, nrows=rows, ncols=cols)
    ifg.close()
    return tiles


def nan_and_mm_convert(ifg, params):
    """
    Perform millimetre and nan conversion on interferogram data

    :param Ifg instance ifg: Interferogram class instance
    :param dict params: Dictionary of parameters

    :return: None, data modified internally
    """
    nan_conversion = params[cf.NAN_CONVERSION]
    if nan_conversion:  # nan conversion happens here in networkx mst
        # if not ifg.nan_converted:
        ifg.nodata_value = params[cf.NO_DATA_VALUE]
        ifg.convert_to_nans()
    if not ifg.mm_converted:
        ifg.convert_to_mm()


def cell_size(lat, lon, x_step, y_step):
    # pylint: disable=invalid-name
    """
    Converts X|Y_STEP in degrees to X & Y cell size in metres.
    This function depends on PyProj/PROJ4 to implement the function

    :param float lat: Latitude in degrees
    :param float lon: Longitude in degrees
    :param float x_step: Horizontal step size in degrees
    :param float y_step: Vertical step size in degrees

    :return: tuple of X and Y cell size floats
    :rtype: tuple
    """
    if lat > 84.0 or lat < -80:
        msg = "No UTM zone for polar region: > 84 degrees N or < 80 degrees S. " \
              "Provided values are lat: "+str(lat) +" long: " +str(lon)
        raise ValueError(msg)

    zone = _utm_zone(lon)
    p0 = pyproj.Proj(proj='latlong', ellps='WGS84')
    p1 = pyproj.Proj(proj='utm', zone=zone, ellps='WGS84')

    x0, y0 = pyproj.transform(p0, p1, lon, lat, errcheck=True)
    x1, y1 = pyproj.transform(p0, p1, lon + x_step, lat + y_step, errcheck=True)

    return tuple(abs(e) for e in (x1 - x0, y1 - y0))


def _utm_zone(longitude):
    """
    Returns basic UTM zone for given longitude in degrees.
    Currently does NOT handle the sub-zoning around Scandanavian countries.
    See http://www.dmap.co.uk/utmworld.htm
    """
    if longitude == 180:
        return 60.0
    return floor((longitude + 180) / 6.0) + 1


class PrereadIfg:
    """
    Convenience class for handling pre-calculated ifg params
    """
    # pylint: disable=too-many-arguments
    # pylint: disable=too-many-instance-attributes
    def __init__(self, path, tmp_path, nan_fraction, first, second, time_span,
                 nrows, ncols, metadata):
        self.path = path
        self.tmp_path = tmp_path
        self.nan_fraction = nan_fraction
        self.first = first
        self.second = second
        self.time_span = time_span
        self.nrows = nrows
        self.ncols = ncols
        self.shape = (nrows, ncols)
        self.metadata = metadata


def _prep_ifg(ifg_path, params):
    """
    Wrapper for reading an interferogram file and creating an Ifg object

    :param str ifg_path: Interferogram file path
    :param dict params: Dictionary of configuration parameters

    :return: ifg: Interferogram class instance
    :rtype: xxxx (eg flaot)
    """
    # Only used in pyrate.scripts.run_pyrate?
    ifg = Ifg(ifg_path)
    ifg.open()
    nan_and_mm_convert(ifg, params)
    return ifg


def save_numpy_phase(ifg_paths, params):
    """
    Save interferogram phase data as numpy array file on disk.

    :param list ifg_paths: List of strings for interferogram paths
    :param list tiles: List of pyrate.shared.Tile instances
    :param dict params: Dictionary of configuration parameters

    :return: None, file saved to disk
    """
    tiles = params['tiles']
    outdir = params[cf.TMPDIR]
    if not os.path.exists(outdir):
        mkdir_p(outdir)
    for ifg_path in mpiops.array_split(ifg_paths):
        ifg = Ifg(ifg_path)
        ifg.open()
        phase_data = ifg.phase_data
        bname = basename(ifg_path).split('.')[0]
        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)
            np.save(file=join(outdir, phase_file),
                    arr=p_data)
        ifg.close()
    mpiops.comm.barrier()


def get_geotiff_header_info(ifg_path):
    """
    Return information from a geotiff interferogram header using GDAL methods.

    :param str ifg_path: path to interferogram geotiff file

    :return: md: PyRate metadata
    :rtype: list
    :return: gt: GDAL geotransform for the data
    :rtype: list
    :return: wkt: GDAL projection information for the data
    :rtype: list
    """
    ds = gdal.Open(ifg_path)
    md = ds.GetMetadata()  # get metadata for writing on output tifs
    gt = ds.GetGeoTransform()  # get geographical bounds of data
    wkt = ds.GetProjection()  # get projection of data
    ds = None  # close dataset
    return gt, md, wkt


def warp_required(xlooks, ylooks, crop):
    """
    Check if a crop or multi-look operation is required.

    :param int xlooks: Resampling/multi-looking in x dir
    :param int ylooks: Resampling/multilooking in y dir
    :param int crop: Interferogram crop option

    :return: True if params show rasters need to be cropped and/or resized
    :rtype: bool
    """
    if xlooks > 1 or ylooks > 1:
        return True
    if crop is None:
        return False
    return True


def check_correction_status(ifgs, meta):  # pragma: no cover
    """
    Generic function for checking if a correction has already been performed
    in a previous run by interrogating PyRate meta data entries

    :param preread_ifgs: Dictionary of pre-read interferogram information
    :param str meta: Meta data flag to check for

    :return: True if correction has been performed, otherwise False
    :rtype: bool
    """
    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]

    for ifg in ifgs:
        if not ifg.is_open:
            ifg.open()

    flags = [meta in ifg.meta_data for ifg in ifgs]
    if all(flags):
        log.info('Skipped: interferograms already corrected')
        return True
    elif not all(flags) and any(flags):
        log.debug('Detected mix of corrected and uncorrected interferograms')
        for i, flag in zip(ifgs, flags):
            if flag:
                msg = '{}: correction detected'.format(i.data_path)
            else:
                msg = '{}: correction NOT detected'.format(i.data_path)
            log.debug(msg)
            close_all(ifgs)
            raise CorrectionStatusError(msg)
    else:
        log.debug('Calculating corrections')
        close_all(ifgs)
        return False


class CorrectionStatusError(Exception):
    """
    Generic class for correction status errors.
    """


def extract_epochs_from_filename(filename_with_epochs: str) -> List[str]:
    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))
    return src_epochs


def mpi_vs_multiprocess_logging(step, params):
    if mpiops.size > 1:  # Over-ride input options if this is an MPI job
        log.info(f"Running '{step}' step with MPI using {mpiops.size} processes")
        log.warning("Disabling joblib parallel processing (setting parallel = 0)")
        params[cf.PARALLEL] = 0
    else:
        if params[cf.PARALLEL] == 1:
            log.info(f"Running '{step}' step in parallel using {params[cf.PROCESSES]} processes")
        else:
            log.info(f"Running '{step}' step in serial")


def dem_or_ifg(data_path):
    """
    Returns an Ifg or DEM class object from input geotiff file.

    :param str data_path: file path name

    :return: Interferogram or DEM object from input file
    :rtype: Ifg or DEM class object
    """
    ds = gdal.Open(data_path)
    md = ds.GetMetadata()
    if ifc.FIRST_DATE in md:  # ifg
        return Ifg(data_path)
    else:
        return DEM(data_path)


def join_dicts(dicts: List[dict]) -> dict:
    """
    Function to concatenate dictionaries
    """
    if dicts is None:  # pragma: no cover
        return {}
    assembled_dict = {k: v for D in dicts for k, v in D.items()}
    return assembled_dict


def tiles_split(func, params, *args, **kwargs):
    tiles = params[cf.TILES]
    process_tiles = mpiops.array_split(tiles)
    if params[cf.PARALLEL]:
        Parallel(n_jobs=params[cf.PROCESSES], verbose=joblib_log_level(cf.LOG_LEVEL))(
            delayed(func)(t, params, *args, **kwargs) for t in process_tiles)
    else:
        for t in process_tiles:
            func(t, params, *args, **kwargs)
    mpiops.comm.barrier()


def output_tiff_filename(inpath, outpath):
    """
    Output geotiff filename for a given input filename.

    :param str inpath: path of input file location
    :param str outpath: path of output file location

    :return: Geotiff filename for the given file.
    :rtype: str
    """
    fname, ext = os.path.basename(inpath).split('.')
    outpath = os.path.dirname(inpath) if outpath is None else outpath
    if ext == 'tif':
        name = os.path.join(outpath, fname + '.tif')
    else:
        name = os.path.join(outpath, fname + '_' + ext + '.tif')
    return name