import re
import glob
import os
import time
from datetime import datetime

import utm
import scipy.io
import numpy as num

from kite import util

__all__ = ['Gamma', 'Matlab', 'ISCE', 'GMTSAR', 'ROI_PAC', 'SARscape', 'SNAP_Gamma']

try:
    from osgeo import gdal
    __all__.append('LiCSAR')
    __all__.append('ARIA')
except ImportError:
    pass


d2r = num.pi/180.
km = 1e3
op = os.path

LAMBDA_SENTINEL = 0.055465763


def check_required(required, params):
    for r in required:
        if r not in params:
            return False
    return True


def safe_cast(val, to_type, default=None):
    try:
        return to_type(val)
    except (ValueError, TypeError):
        return default


class HeaderError(Exception):
    pass


class AttribDict(dict):

    def __getattr__(self, item):
        return self[item]

    def __setattr__(self, item, value):
        self[item] = value


class SceneIO(object):
    """ Prototype class for SARIO objects. """
    def __init__(self, scene=None):
        if scene is not None:
            self._log = scene._log.getChild('IO/%s' % self.__class__.__name__)
        else:
            import logging
            self._log = logging.getLogger('SceneIO/%s'
                                          % self.__class__.__name__)

        self.container = AttribDict(
            phi=0.,    # Look orientation counter-clockwise angle from east
            theta=0.,  # Look elevation angle (up from horizontal in degree)
                       # 90 deg North
            displacement=None,  # Displacement towards LOS
            frame=AttribDict(
                llLon=None,  # Lower left corner latitude
                llLat=None,  # Lower left corner londgitude
                dN=None,   # Pixel delta in north, meter or degree
                dE=None,   # Pixel delta in east, meter or degree
                spacing='meter',  # Pixel spacing unit
            ),
            # Meta information
            meta=AttribDict(
                title=None,
                orbital_node=None,
                satellite_name=None,
                wavelength=None,
                time_master=None,
                time_slave=None
            ),
            # All extra information
            extra={}
        )

    def read(self, filename, **kwargs):
        """ Read function of the file format

        :param filename: file to read
        :type filename: string
        :param kwargs: Keyword arguments
        :type kwargs: {dict}
        """
        raise NotImplementedError('read not implemented')

    def write(self, filename, **kwargs):
        """ Write method for IO

        :param filename: file to write to
        :type filename: string
        :param **kwargs: Keyword arguments
        :type **kwargs: {dict}
        """
        raise NotImplementedError('write not implemented')

    def validate(self, filename, **kwargs):
        """ Validate file format

        :param filename: file to validate
        :type filename: string
        :returns: Validation
        :rtype: {bool}
        """
        pass
        raise NotImplementedError('validate not implemented')


class Matlab(SceneIO):
    """
    Variable naming conventions for Matlab :file:`.mat` container:

        ================== ==================== ===================== =====
        Property           Matlab ``.mat`` name type                  unit
        ================== ==================== ===================== =====
        Scene.displacement ``ig_``              n x m array           [m]
        Scene.phi          ``phi``              float or n x m array  [rad]
        Scene.theta        ``theta``            float or n x m array  [rad]
        Scene.frame.x      ``xx``               n x 1 vector          [m]
        Scene.frame.y      ``yy``               m x 1 vector          [m]
        Scene.utm_zone     ``utm_zone``         str ('33T')
        ================== ==================== ===================== =====

    Displacement is expected to be in meters. Note that the displacement maps
    could also be pixel offset maps rather than unwrapped SAR interferograms.
    For SAR azimuth pixel offset maps calculate ``phi`` from the heading
    direction and set ``theta=0.``. For SAR range pixel offsets use the same
    LOS angles as for InSAR.
    """
    def validate(self, filename, **kwargs):
        if filename[-4:] == '.mat':
            return True
        else:
            return False
        try:
            variables = self.io.whosmat(filename)
            if len(variables) > 50:
                return False
            return True
        except ValueError:
            return False

    def read(self, filename, **kwargs):
        c = self.container

        mat = scipy.io.loadmat(filename)
        utm_e = None
        utm_n = None
        utm_zone = None
        utm_zone_letter = None
        phi0 = None
        theta0 = None

        for mat_k, v in mat.items():
            for io_k in c.keys():
                if io_k in mat_k:
                    c[io_k] = num.rot90(mat[mat_k])
                elif 'ig_' in mat_k:
                    c.displacement = num.rot90(mat[mat_k])
                elif 'xx' in mat_k:
                    utm_e = mat[mat_k].flatten()
                elif 'yy' in mat_k:
                    utm_n = mat[mat_k].flatten()
                elif 'utm_zone' in mat_k:
                    utm_zone = int(mat['utm_zone'][0][:-1])
                    utm_zone_letter = str(mat['utm_zone'][0][-1])
                elif 'phi' in mat_k:
                    phi0 = mat[mat_k].flatten()
                elif 'theta' in mat_k:
                    theta0 = mat[mat_k].flatten()

        if len(theta0) == 1:
            c.theta = num.ones(num.shape(c.displacement)) * theta0

        if len(theta0) == 1:
            c.phi = num.ones(num.shape(c.displacement)) * phi0

        if utm_zone is None:
            utm_zone = 33
            utm_zone_letter = 'N'
            self._log.warning(
                'Variable utm_zone not defined. Defaulting to UTM Zone %d%s!'
                % (utm_zone, utm_zone_letter))

        if not (num.all(utm_e) or num.all(utm_n)):
            self._log.warning(
                'Could not find referencing UTM vectors in .mat file!')
            utm_e = num.linspace(100000, 110000, c.displacement.shape[0])
            utm_n = num.linspace(1100000, 1110000, c.displacement.shape[1])

        if utm_e.min() < 1e4 or utm_n.min() < 1e4:
            utm_e *= km
            utm_n *= km

        c.frame.dE = num.abs(utm_e[1] - utm_e[0])
        c.frame.dN = num.abs(utm_n[1] - utm_n[0])
        try:
            c.frame.llLat, c.frame.llLon =\
                utm.to_latlon(utm_e.min(), utm_n.min(),
                              utm_zone, utm_zone_letter)

        except utm.error.OutOfRangeError:
            self._log.warning(
                'Could not interpret spatial vectors,'
                ' referencing to 0, 0 (lat, lon)')
            c.frame.llLat, c.frame.llLon = (0., 0.)
        return c


class Gamma(SceneIO):
    """

    Reading geocoded displacement maps (unwrapped igs) originating
        from GAMMA software.

    .. note :: Expects:

        * [:file:`*`] Binary file from Gamma with displacement in radians
        * [:file:`*.slc.par`] If you want to translate radians to
          meters using the `radar_frequency`.
        * [:file:`*par`] Parameter file, describing ``corner_lat, corner_lon,
          nlines, width, post_lat, post_lon`` or ``post_north, post_east,
          corner_east, corner_north, nlines, width``.
        * [:file:`*theta*`, :file:`*phi*`] Two look vector files,
          generated by GAMMA command ``look_vector``.

    .. warning ::

        * Data has to be georeferenced to latitude/longitude or UTM!
        * Look vector files - expected to have a particular name
    """
    @staticmethod
    def _parseParameterFile(filename):
        params = {}
        rc = re.compile(r'^(\w*):\s*([a-zA-Z0-9+-.*]*\s[a-zA-Z0-9_]*).*')

        with open(filename, mode='r') as par:
            for line in par:
                parsed = rc.match(line)
                if parsed is None:
                    continue

                groups = parsed.groups()
                params[groups[0]] = safe_cast(groups[1], float,
                                              default=groups[1].strip())
        return params

    def _getParameters(self, path, log=False):
        required_utm = ['post_north', 'post_east', 'corner_east',
                        'corner_north', 'nlines', 'width']
        required_lat_lon = ['corner_lat', 'corner_lon', 'nlines',
                            'width', 'post_lat', 'post_lon']

        path = op.dirname(op.realpath(path))
        par_files = glob.glob('%s/*par' % path)

        for file in par_files:
            params = self._parseParameterFile(file)

            if check_required(required_utm, params)\
               or check_required(required_lat_lon, params):
                if not log:
                    self._log.info('Found parameter file %s' % file)
                return params

        raise ImportError(
                    'Parameter file does not hold required parameters')

    def _getSLCParameters(self, path):
        required_params = ['radar_frequency']
        path = op.dirname(op.realpath(path))
        par_files = glob.glob('%s/*slc.par' % path)

        for file in par_files:
            params = self._parseParameterFile(file)

            if check_required(required_params, params):
                self._log.info('Found SLC parameter file %s' % file)
                return params

        raise ImportError('Could not find SLC parameter file *.slc.par'
                          ' with parameters %s' % required_params)

    def validate(self, filename, **kwargs):
        try:
            par_file = kwargs.pop('par_file', filename)
            self._getParameters(par_file)
            return True
        except ImportError:
            return False

    def _getLOSAngles(self, filename, pattern):
        path = op.dirname(op.realpath(filename))
        phi_files = glob.glob('%s/%s' % (path, pattern))
        if len(phi_files) == 0:
            self._log.warning('Could not find LOS file %s, '
                              'defaulting to angle to 0.' % pattern)
            return 0.
        elif len(phi_files) > 1:
            self._log.warning('Found multiple LOS files %s, '
                              'defaulting to angle 0.' % pattern)
            return 0.

        filename = phi_files[0]
        self._log.info('Loading LOS %s from %s' % (pattern, filename))
        return num.memmap(filename, mode='r', dtype='>f4')

    def read(self, filename, **kwargs):
        """
        :param filename: Gamma software parameter file
        :type filename: str
        :param par_file: Corresponding parameter (:file:`*par`) file.
                         (optional)
        :type par_file: str
        :returns: Import dictionary
        :rtype: dict
        :raises: ImportError
        """
        par_file = kwargs.pop('par_file', filename)

        params = self._getParameters(par_file, log=True)

        try:
            params_slc = self._getSLCParameters(par_file)
        except ImportError as e:
            raise e

        fill = None

        ncols = int(params['width'])
        nlines = int(params['nlines'])
        radar_frequency = float(params_slc.get('radar_frequency', None))

        displ = num.fromfile(filename, dtype='>f4')
        # Resize array if last line is not scanned completely
        if (displ.size % ncols) != 0:
            fill = num.empty(ncols - displ.size % ncols)
            fill.fill(num.nan)
            displ = num.append(displ, fill)

        displ = displ.reshape(nlines, ncols)
        displ[displ == -0.] = num.nan
        displ = num.flipud(displ)

        if radar_frequency is not None:
            radar_frequency = float(radar_frequency)
            self._log.info('Scaling displacement by radar_frequency %f GHz'
                           % (radar_frequency/1e9))
            wavelength = util.C / radar_frequency
            displ /= -4*num.pi
            displ *= wavelength

        else:
            wavelength = 'None'
            self._log.warning(
                'Could not determine radar_frequency from *.slc.par file!'
                ' Leaving displacement to radians.')

        phi = self._getLOSAngles(filename, '*phi*')
        theta = self._getLOSAngles(filename, '*theta*')
        theta = theta

        if isinstance(phi, num.ndarray):
            phi = phi.reshape(nlines, ncols)
            phi = num.flipud(phi)
        if isinstance(theta, num.ndarray):
            theta = theta.reshape(nlines, ncols)
            theta = num.flipud(theta)

        if fill is not None:
            theta = num.append(theta, fill)
            phi = num.append(phi, fill)

        c = self.container

        c.displacement = displ
        c.theta = theta
        c.phi = phi

        c.meta.wavelength = wavelength
        c.meta.title = params.get('title', 'None')

        c.bin_file = filename
        c.par_file = par_file

        if params['DEM_projection'] == 'UTM':
            utm_zone = params['projection_zone']
            try:
                utm_zone_letter = utm.latitude_to_zone_letter(
                    params['center_latitude'])
            except ValueError:
                self._log.warning('Could not parse UTM Zone letter,'
                                  ' defaulting to N!')
                utm_zone_letter = 'N'

            self._log.info('Using UTM reference: Zone %d%s'
                           % (utm_zone, utm_zone_letter))

            dN = params['post_north']
            dE = params['post_east']

            utm_corn_e = params['corner_east']
            utm_corn_n = params['corner_north']

            utm_corn_eo = utm_corn_e + dE * displ.shape[1]
            utm_corn_no = utm_corn_n + dN * displ.shape[0]

            utm_e = num.linspace(utm_corn_e, utm_corn_eo, displ.shape[1])
            utm_n = num.linspace(utm_corn_n, utm_corn_no, displ.shape[0])

            llLat, llLon = utm.to_latlon(utm_e.min(), utm_n.min(),
                                         utm_zone, utm_zone_letter)

            c.frame.llLat = llLat
            c.frame.llLon = llLon

            c.frame.dE = abs(dE)
            c.frame.dN = abs(dN)

        else:
            self._log.info('Using Lat/Lon reference')
            c.frame.spacing = 'degree'
            c.frame.llLat = params['corner_lat'] \
                + params['post_lat'] * nlines
            c.frame.llLon = params['corner_lon']
            c.frame.dE = abs(params['post_lon'])
            c.frame.dN = abs(params['post_lat'])

        return c


class ROI_PAC(SceneIO):
    """
    .. note:: Expects:

        * Binary file from ROI_PAC (:file:`*`)
        * Parameter file (:file:`<binary_file>.rsc`),
          describing ``WIDTH, FILE_LENGTH, X_FIRST, Y_FIRST, X_STEP,
          Y_STEP, WAVELENGTH``
        * If the georeferencing is in UTM coordinates, further needed
          entries in parameter file are 'X_UNIT' and 'Y_UNIT' that give
          'meters' and 'LAT_REF3' as well as 'LON_REF3'.

        The unwrapped displacement is expected in radians and will be scaled
        to meters by ``WAVELENGTH`` parsed from the :file:`*.rsc` file.

    """

    def validate(self, filename, **kwargs):
        try:
            par_file = kwargs.pop('par_file',
                                  self._getParameterFile(filename))
            self._parseParameterFile(par_file)
            return True
        except ImportError:
            return False

    def _getParameterFile(self, bin_file):
        par_file = op.realpath(bin_file) + '.rsc'
        try:
            self._parseParameterFile(par_file)
            self._log.info('Found parameter file %s' % par_file)
            return par_file
        except (ImportError, IOError):
            raise ImportError('Could not find ROI_PAC parameter file (%s)'
                              % par_file)

    @staticmethod
    def _parseParameterFile(par_file):
        params = {}
        required_L0 = ['WIDTH', 'FILE_LENGTH', 'X_FIRST', 'Y_FIRST', 'X_STEP',
                       'Y_STEP', 'WAVELENGTH']
        required_utm = ['X_UNIT', 'LAT_REF1', 'LON_REF1']

        rc = re.compile(r'([\w]*)\s*([\w.+-]*)')
        with open(par_file, 'r') as par:
            for line in par:
                parsed = rc.match(line)
                if parsed is None:
                    continue
                groups = parsed.groups()
                params[groups[0]] = safe_cast(groups[1], float,
                                              default=groups[1].strip())

        if check_required(required_L0, params):
            if check_required(required_utm, params):
                geo_ref = 'all'
                return params, geo_ref
            else:
                geo_ref = 'latlon'
                return params, geo_ref

        raise ImportError(
            'Parameter file %s does not hold the basic \
             required parameters' % par_file)

    def read(self, filename, **kwargs):
        """
        :param filename: ROI_PAC binary file
        :type filename: str
        :param par_file: Corresponding parameter (:file:`*rsc`) file.
                         (optional)
        :type par_file: str
        :returns: Import dictionary
        :rtype: dict
        :raises: ImportError
        """
        par_file = kwargs.pop('par_file', self._getParameterFile(filename))

        par, geo_ref = self._parseParameterFile(par_file)
        nlines = int(par['FILE_LENGTH'])
        ncols = int(par['WIDTH'])
        wavelength = par['WAVELENGTH']
        heading = par['HEADING_DEG']
        if geo_ref == 'latlon':
            lat_ref = par['Y_FIRST']
            lon_ref = par['X_FIRST']
        elif geo_ref == 'all':
            lat_ref = par['LAT_REF3']
            lon_ref = par['LON_REF3']

        look_ref1 = par['LOOK_REF1']
        look_ref2 = par['LOOK_REF2']
        look_ref3 = par['LOOK_REF3']
        look_ref4 = par['LOOK_REF4']

        utm_zone_letter = utm.latitude_to_zone_letter(lat_ref)
        utm_zone = utm.latlon_to_zone_number(lat_ref, lon_ref)

        look = num.mean(
            num.array([look_ref1, look_ref2, look_ref3, look_ref4]))

        data = num.memmap(filename, dtype='<f4')
        data = data.reshape(nlines, ncols*2)

        displ = data[:, ncols:]
        displ = num.flipud(displ)
        displ[displ == -0.] = num.nan
        displ = displ / (4.*num.pi) * wavelength

        z_scale = par.get('Z_SCALE', 1.)
        z_offset = par.get('Z_OFFSET', 0.)
        displ += z_offset
        displ *= z_scale

        c = self.container

        c.displacement = displ
        c.theta = num.deg2rad(90. - look)
        c.phi = num.deg2rad(-heading + 180.)

        c.meta.title = par.get('TITLE', 'None')
        c.meta.wavelength = par['WAVELENGTH']
        c.bin_file = filename
        c.par_file = par_file

        if geo_ref == 'all':
            if par['X_UNIT'] == 'meters':
                c.frame.spacing = 'meter'
                c.frame.dE = par['X_STEP']
                c.frame.dN = -par['Y_STEP']
                geo_ref = 'utm'

            elif par['X_UNIT'] == 'degree':
                c.frame.spacing = 'degree'
                geo_ref = 'latlon'

        elif geo_ref == 'latlon':
            self._log.info('Georeferencing is in Lat-Lon [degrees].')
            c.frame.spacing = 'degree'
            c.frame.llLat = par['Y_FIRST'] + par['Y_STEP'] * nlines
            c.frame.llLon = par['X_FIRST']

            # c_utm_0 = utm.from_latlon(lat_ref, lon_ref)
            # c_utm_1 = utm.from_latlon(lat_ref + par['Y_STEP'],
            #                           lon_ref + par['X_STEP'])

            # c.frame.dE = c_utm_1[0] - c_utm_0[0]
            # c.frame.dN = abs(c_utm_1[1] - c_utm_0[1])
            c.frame.dE = par['X_STEP']
            c.frame.dN = -par['Y_STEP']

        elif geo_ref == 'utm':
            self._log.info('Georeferencing is in UTM (zone %d%s)',
                           utm_zone, utm_zone_letter)
            y_ll = par['Y_FIRST'] + par['Y_STEP'] * nlines
            c.frame.llLat, c.frame.llLon = utm.to_latlon(
                par['X_FIRST'], y_ll, utm_zone,
                zone_letter=utm_zone_letter)

        return self.container


class ISCEXMLParser(object):
    def __init__(self, filename):
        import xml.etree.ElementTree as ET
        self.root = ET.parse(filename).getroot()

    @staticmethod
    def type_convert(value):
        for t in (float, int, str):
            try:
                return t(value)
            except ValueError:
                continue
        raise ValueError('Could not convert value')

    def getProperty(self, name):
        name = name.lower()

        for child in self.root.iter():
            child_name = child.get('name')
            if isinstance(child_name, str):
                child_name = child_name.lower()
            if child_name == name.lower():
                if child.tag == 'property':
                    return self.type_convert(child.find('value').text)
                elif child.tag == 'component':
                    values = {}
                    for prop in child.iter('property'):
                        values[prop.get('name')] =\
                            self.type_convert(prop.find('value').text)
                    return values
        return None


class ISCE(SceneIO):
    """
    Reading geocoded, unwraped displacement maps
    processed with ISCE software (https://winsar.unavco.org/isce.html).

    .. note :: Expects:

        * Unwrapped displacement binary (:file:`*.unw.geo`)
        * Metadata XML (:file:`*.unw.geo.xml`)
        * LOS binary data (:file:`*.rdr.geo`)

    .. note ::

        When using ``gdal_translate`` to crop the scene, use the argument
        ``-co SCHEME=BIL`` to make the output

    .. note ::

        Data are in radians but no transformation to
        meters yet, as ``wavelength`` or at least sensor name is not
        provided in the XML file.
    """
    def validate(self, filename, **kwargs):
        try:
            self._getDisplacementFile(filename)
            self._getLOSFile(filename)
            return True
        except ImportError:
            return False

    def _getLOSFile(self, path):
        if not op.isdir(path):
            path = op.dirname(path)
        rdr_files = glob.glob(op.join(path, '*.rdr.geo'))

        if len(rdr_files) == 0:
            raise ImportError('Could not find LOS file (*.rdr.geo)')

        rdr_file = rdr_files[0]
        self._log.info('Found LOS file: %s', rdr_file)
        return rdr_file

    def _getDisplacementFile(self, path):
        if op.isfile(path):
            disp_file = path
        else:
            files = glob.glob(op.join(path, '*.unw.geo'))
            if len(files) == 0:
                raise ImportError('Could not find displacement file '
                                  '(.unw.geo) at %s', path)
            disp_file = files[0]

        if not op.isfile('%s.xml' % disp_file):
            raise ImportError('Could not find displacement XML file '
                              '(%s.unw.geo.xml)' % op.basename(disp_file))
        self._log.info('Found Displacement file: %s', disp_file)
        return disp_file

    def read(self, path, **kwargs):
        path = op.abspath(path)
        c = self.container

        xml_file = self._getDisplacementFile(path) + '.xml'
        self._log.info('Parsing ISCE XML file %s' % xml_file)
        isce_xml = ISCEXMLParser(xml_file)

        coord_lon = isce_xml.getProperty('coordinate1')
        coord_lat = isce_xml.getProperty('coordinate2')
        c.frame.dN = num.abs(coord_lat['delta'])
        c.frame.dE = num.abs(coord_lon['delta'])

        nlon = int(coord_lon['size'])
        nlat = int(coord_lat['size'])

        c.frame.spacing = 'degree'
        c.frame.llLat = coord_lat['startingvalue'] +\
            (nlat * coord_lat['delta'])
        c.frame.llLon = coord_lon['startingvalue']

        displ = num.memmap(self._getDisplacementFile(path),
                           dtype='<f4')\
            .reshape(nlat, nlon*2)[:, nlon:]

        displ = num.flipud(displ)
        displ[displ == 0.] = num.nan
        c.displacement = displ

        los_file = self._getLOSFile(path)
        los_data = num.fromfile(los_file, dtype='<f4')\
            .reshape(nlat*2, nlon)

        theta = num.flipud(los_data[0::2, :])
        phi = num.flipud(los_data[1::2, :])

        def los_is_degree():
            return num.abs(theta).max() > num.pi or num.abs(phi).max() > num.pi

        if not los_is_degree():
            raise ImportError(
                'The LOS file (%s) seems to be in radians! '
                'Change it to degree!' % op.basename(los_file))

        phi[phi == 0.] = num.nan
        theta[theta == 0.] = num.nan

        phi *= d2r
        theta *= d2r

        phi = num.pi/2 + phi
        theta = num.pi/2 - theta

        c.phi = phi
        c.theta = theta

        return c


class GMTSAR(SceneIO):
    """
    Reading GMTSAR grid files.

    .. note ::

        Expects:

        * Displacement grid (NetCDF, :file:`*los_ll.grd`) in meter
          (in case use "gmt grdmath los_cm_ll.grd 0.01 MUL = los_m_ll.grd')
        * LOS binary data (see instruction, :file:`*.los.enu`)

    Calculate the corresponding unit look vectors with GMT5SAR ``SAT_look``:

    .. code-block:: sh

        gmt grd2xyz los_ll.grd | gmt grdtrack -Gdem.grd | \\
        awk {'print $1, $2, $4'} | \\
        SAT_look 20050731.PRM -bos > 20050731.los.enu
    """
    def validate(self, filename, **kwargs):
        try:
            if self._getDisplacementFile(filename)[-4:] == '.grd':
                return True
        except ImportError:
            return False
        return False

    def _getLOSFile(self, path):
        if not op.isdir(path):
            path = op.dirname(path)
        los_files = glob.glob(op.join(path, '*.los.*'))
        if len(los_files) == 0:
            self._log.warning(GMTSAR.__doc__)
            raise ImportError('Could not find LOS file (*.los.*)')
        los_file = los_files[0]
        self._log.debug('Found LOS file: %s', los_file)
        return los_file

    def _getDisplacementFile(self, path):
        if op.isfile(path):
            return path
        else:
            files = glob.glob(op.join(path, '*.grd'))
            if len(files) == 0:
                raise ImportError('Could not find displacement file '
                                  '(*.grd) at %s', path)
            disp_file = files[0]
        self._log.debug('Found Displacement file: %s', disp_file)
        return disp_file

    def read(self, path, **kwargs):
        from scipy.io import netcdf
        path = op.abspath(path)
        c = self.container

        grd = netcdf.netcdf_file(self._getDisplacementFile(path),
                                 mode='r', version=2)
        displ = grd.variables['z'][:].copy()
        c.displacement = displ
        shape = c.displacement.shape
        # LatLon
        c.frame.spacing = 'degree'
        c.frame.llLat = grd.variables['lat'][:].min()
        c.frame.llLon = grd.variables['lon'][:].min()

        c.frame.dN = (grd.variables['lat'][:].max() -
                      c.frame.llLat) / shape[0]
        c.frame.dE = (grd.variables['lon'][:].max() -
                      c.frame.llLon) / shape[1]

        # Theta and Phi
        try:
            los = num.memmap(self._getLOSFile(path), dtype='<f4')
            e = los[3::6].copy().reshape(shape)
            n = los[4::6].copy().reshape(shape)
            u = los[5::6].copy().reshape(shape)

            phi = num.arctan(n/e)
            theta = num.arcsin(u)
            # phi[n < 0] += num.pi

            c.phi = phi
            c.theta = theta
        except ImportError:
            self._log.warning(self.__doc__)
            self._log.warning('Defaulting theta to pi/2 and phi to 0.')
            c.theta = num.pi/2
            c.phi = 0.
        return c


class SARscape(SceneIO):
    """
    Reading SARscape :file:`*_disp` files.

    .. note ::

        Expects:

        * Header file in :file:`*_disp.hdr`
        * Displacement data in cm in :file:`*_disp`
        * LOS data in :file:`*disp_ILOS` and :file:`*disp_ALOS` files.
    """
    def read(self, filename, **kwargs):
        header = self.parseHeaderFile(filename)

        def load_data(filename):
            self._log.debug('Loading %s' % filename)
            return num.flipud(
                num.fromfile(filename, dtype=num.float32)
                .reshape((header.lines, header.samples)))

        displacement = load_data(filename)
        theta_file, phi_file = self.getLOSFiles(filename)

        if not theta_file:
            theta = num.full_like(displacement, 0.)
        else:
            theta = load_data(theta_file)
            theta = num.deg2rad(theta)

        if not phi_file:
            phi = num.full_like(displacement, num.pi/2)
        else:
            phi = load_data(phi_file)
            phi = num.pi/2 - num.rad2deg(phi)

        c = self.container
        c.displacement = displacement
        c.phi = phi
        c.theta = theta

        map_info = header.map_info
        c.frame.dE = float(map_info[5])
        c.frame.dN = dN = float(map_info[6])
        c.frame.spacing = 'meter'

        c.frame.llLat, c.frame.llLon = utm.to_latlon(
            float(map_info[3]) - header.lines * dN,
            float(map_info[4]),
            zone_number=int(map_info[7]),
            northern=True if map_info[8] == 'Northern' else False)

        return c

    def parseHeaderFile(self, filename):
        hdr_file = self._getHDRFile(filename)
        conf = re.compile(r'^(.+)\s+=\s+(.+)\n', re.MULTILINE)

        header = AttribDict()
        with open(hdr_file) as f:
            s = f.read()

            linebreaks = re.compile(r'{(.+)\n?(.+)}')
            s = linebreaks.sub(r'{ \g<1> \g<2> }', s)

            for match in conf.finditer(s):
                groups = match.groups()
                key = groups[0].strip().replace(' ', '_')
                value = groups[1].strip()
                try:
                    value = int(value)
                except ValueError:
                    pass

                header[key] = value

            header.map_info = header.map_info.strip('{} ').split(', ')
            if not len(header.map_info) == 11:
                raise HeaderError('`map info` header is not consistent!')
            if header.map_info[0] != 'UTM':
                raise HeaderError('`map info` is not UTM!')

        return header

    def getLOSFiles(self, filename):
        ilos_file = op.abspath(filename + '_ILOS')
        if not op.exists(ilos_file):
            self._log.warning('Could not find ILOS file! (%s)' % ilos_file)
            ilos_file = False

        alos_file = op.abspath(filename + '_ALOS')
        if not op.exists(alos_file):
            self._log.warning('Could not find ALOS file! (%s)' % alos_file)
            alos_file = False
        return ilos_file, alos_file

    def _getHDRFile(self, filename):
        hdr_file = op.abspath(op.splitext(filename)[0] + '.hdr')
        if not op.exists(hdr_file):
            raise OSError('SARscape .hdr file not found (%s)' % hdr_file)
        return hdr_file

    def validate(self, filename, **kwargs):
        val = re.compile(r'SARscape|ENVI Standard', re.MULTILINE)
        try:
            hdr_file = self._getHDRFile(filename)
        except OSError:
            return False

        with open(hdr_file) as f:
            res = val.search(f.read())
            if res is not None:
                return True
            return False


class LiCSAR(SceneIO):
    '''
    Import unwrapped Geotiffs from the
    `COMET LiCSAR Portal <https://comet.nerc.ac.uk/COMET-LiCS-portal/>`_.

    .. note ::

        Requires the Python package
        `gdal/osgeo <https://pypi.org/project/GDAL/>`_! Or through

        Expects:

        * Unwrapped geotiff in :file:`*.unw.tif`
        * LOS data in :file:`*.geo.[NEU].tif` files

    See also the download script in :mod:`kite.clients`.
    '''

    def _getLOS(self, filename, component):
        path = op.dirname(filename)
        fn = glob.glob(op.join(path, component))
        if len(fn) != 1:
            raise ImportError('Cannot find LOS vector file %s!' % component)

        dataset = gdal.Open(fn[0], gdal.GA_ReadOnly)
        return self._readBandData(dataset)

    @staticmethod
    def _readBandData(dataset, band=1):
        band = dataset.GetRasterBand(band)
        array = band.ReadAsArray()
        array[array == band.GetNoDataValue()] = num.nan

        return num.flipud(array)

    def read(self, filename, **kwargs):
        dataset = gdal.Open(filename, gdal.GA_ReadOnly)
        georef = dataset.GetGeoTransform()

        llLon = georef[0]
        llLat = georef[3] + dataset.RasterYSize * georef[5]

        c = self.container

        c.frame.spacing = 'degree'
        c.frame.llLat = llLat
        c.frame.llLon = llLon
        c.frame.dE = georef[1]
        c.frame.dN = abs(georef[5])

        displacement = self._readBandData(dataset)
        c.displacement = -displacement / (4*num.pi) * LAMBDA_SENTINEL

        try:
            los_n = self._getLOS(filename, '*.geo.N.tif')
            los_e = self._getLOS(filename, '*.geo.E.tif')
            los_u = self._getLOS(filename, '*.geo.U.tif')
        except ImportError:
            self._log.warning(
                'Cannot find LOS angle files *.geo.[NEU].tif,'
                ' using static Sentinel-1 descending angles.')

            heading = 83.
            incident = 50.

            un = num.sin(d2r*incident) * num.cos(d2r*heading)
            ue = num.sin(d2r*incident) * num.sin(d2r*heading)
            uz = num.cos(d2r*incident)

            los_n = num.full_like(c.displacement, un)
            los_e = num.full_like(c.displacement, ue)
            los_u = num.full_like(c.displacement, uz)

        c.phi = num.arctan2(los_n, los_e)
        c.theta = num.arcsin(los_u)

        c.meta.title = dataset.GetDescription()

        return c

    def validate(self, filename, **kwargs):
        if gdal.IdentifyDriver(filename) is None:
            return False
        return True


class ARIA(SceneIO):
    '''
    Import unwrapped InSAR scenes from the
    `NASA/JPL ARIA <https://aria.jpl.nasa.gov/>`_ GUNW data products.

    .. note ::

        Requires the Python package
        `gdal/osgeo <https://pypi.org/project/GDAL/>`_! Or through

        Expects:

        * Extracted layers: unwrappedPhase, lookAngle, incidenceAngle,
          connectedComponents

    Use ``ariaExtract.py`` to extract the layers:

    .. code-block:: sh

        ariaExtract.py -w ascending -f aria-data.nc -d download \\
        -l unwrappedPhase,incidenceAngle,lookAngle

    '''

    @staticmethod
    def _readBandData(dataset, band=1):
        band = dataset.GetRasterBand(band)
        array = band.ReadAsArray()
        if array.dtype != num.int16 and array.dtype != num.int:
            array[array == band.GetNoDataValue()] = num.nan

        return num.flipud(array)

    @staticmethod
    def _dataset_from_dir(folder):
        files = set(f for f in os.scandir(folder) if f.is_file())
        for f in files:
            if op.splitext(f.name)[-1] == '':
                break
        else:
            raise ImportError('could not load dataset from %s' % folder)

        return gdal.Open(f.path, gdal.GA_ReadOnly)

    def read(self, folder, **kwargs):
        unw_phase = self._dataset_from_dir(op.join(folder, 'unwrappedPhase'))
        georef = unw_phase.GetGeoTransform()

        llLon = georef[0]
        llLat = georef[3] + unw_phase.RasterYSize * georef[5]

        c = self.container

        c.frame.spacing = 'degree'
        c.frame.llLat = llLat
        c.frame.llLon = llLon
        c.frame.dE = georef[1]
        c.frame.dN = abs(georef[5])

        conn_comp = self._dataset_from_dir(op.join(
            folder, 'connectedComponents'))

        displacement = self._readBandData(unw_phase)
        conn_mask = self._readBandData(conn_comp)  # Mask from snaphu
        displacement *= num.where(conn_mask, 1., num.nan)

        c.displacement = displacement / (4*num.pi) * LAMBDA_SENTINEL

        inc_angle = self._dataset_from_dir(op.join(folder, 'incidenceAngle'))
        azi_angle = self._dataset_from_dir(op.join(folder, 'azimuthAngle'))

        c.theta = num.pi/2 - self._readBandData(inc_angle) * d2r
        c.phi = self._readBandData(azi_angle) * d2r

        c.meta.scene_id = op.basename(unw_phase.GetDescription())
        c.meta.scene_title = c.meta.scene_id

        t_slave, t_master = c.meta.scene_id.split('_')
        c.meta.time_master = datetime(*time.strptime(t_master, '%Y%m%d')[:6]) \
            .timestamp()
        c.meta.time_slave = datetime(*time.strptime(t_slave, '%Y%m%d')[:6]) \
            .timestamp()

        c.meta.satellite_name = 'undefined (ARIA)'

        return c

    def validate(self, folder, **kwargs):
        expected_dirs = set(
            ['unwrappedPhase', 'incidenceAngle', 'lookAngle',
             'connectedComponents'])
        if not op.isdir(folder):
            return False

        dirs = set(d.name for d in os.scandir(folder) if d.is_dir())
        if not expected_dirs - dirs:
            return True

        return False


class SNAP_Gamma(SceneIO):
    """SNAP import

    Reading geocoded displacement maps (unwrapped igs) originating
    from SNAP software using the export option.

    When georeferencing the scene, export the incidence angle from ellpsoid in
    the tab **Processing Parameters**.
    Export the unwrapped **displacement** band and
    **incidenceAngleFromEllipsoid** to Gamma format:

    1. File -> Export -> SAR Formats -> Gamma

    Select **Metadata -> Abstracted_Metadata** and
    2. File -> Export -> Other -> Product Metadata.

    .. note :: Expects:

        * [:file:`*`] Binary file from SNAP Gamma Export with displacement in
            radians
        * [:file:`*.Abstracted_Metadata.txt`] Metadata (parameter) file
            If you want to translate radians to
            meters using the `radar_frequency`.
        * [:file:`*par`] Parameter file, describing ``first_near_lat,
        last_near_long, num_output_lines, num_samples_per_line,
                lat_pixel_res, lon_pixel_res, radar_frequency and heading.``
        * [:file:`incidenceAngleFromEllipsoid.rslc`] Incidence angle file.

    .. warning ::

        * Data has to be georeferenced to latitude/longitude or UTM!
        * Look vector files - expected to have a particular name
    """
    @staticmethod
    def _parseParameterFile(filename):
        params = {}
        rc = re.compile(r'([\w]*)\s*([\w.+-]*)')
        with open(filename, mode='r') as par:
            for line in par:
                if line[:29] == 'metadata.Abstracted_Metadata.':
                    parsed = rc.match(line[29:])
                    if parsed is None:
                        continue

                    groups = parsed.groups()
                    params[groups[0]] = safe_cast(
                        groups[1], float,
                        default=groups[1].strip())
        return params

    def _getParameters(self, path, log=False):
        required_utm = (
            'post_north', 'post_east', 'corner_east',
            'corner_north', 'num_output_lines',
            'num_samples_per_line')
        required_lat_lon = (
            'first_near_lat', 'last_near_long',
            'num_output_lines', 'num_samples_per_line',
            'lat_pixel_res', 'lon_pixel_res')

        par_file = glob.glob(op.join(path, '*Abstracted_Metadata.txt'))
        for file in par_file:
            params = self._parseParameterFile(file)
            if check_required(required_utm, params)\
               or check_required(required_lat_lon, params):
                if not log:
                    self._log.info('Found parameter file %s' % file)
                return op.basename(file), params

        raise ImportError(
            'Parameter file %s does not hold required parameters' %
            par_file)

    def validate(self, filename, **kwargs):
        try:
            folder = op.dirname(op.abspath(filename))
            self._getParameters(folder)
            return True
        except ImportError:
            return False

    def _getLOSAngles(self, filename, pattern):
        path = op.dirname(op.realpath(filename))
        phi_files = glob.glob(op.join(path, pattern))
        if len(phi_files) == 0:
            self._log.warning('Could not find LOS file %s, '
                              'defaulting to angle to 0.', pattern)
            return 0.
        elif len(phi_files) > 1:
            self._log.warning('Found multiple LOS files %s, '
                              'defaulting to angle 0.', pattern)
            return 0.

        filename = phi_files[0]
        self._log.info('Loading LOS %s from %s' % (pattern, filename))
        return num.memmap(filename, mode='r', dtype='>f4')

    def read(self, filename, **kwargs):
        """
        :param filename: Gamma software parameter file
        :type filename: str
        :param par_file: Corresponding parameter (:file:`*par`) file.
                         (optional)
        :type par_file: str
        :returns: Import dictionary
        :rtype: dict
        :raises: ImportError
        """
        par_file, params = self._getParameters(
            op.dirname(op.abspath(filename)), log=True)

        ncols = int(params['num_samples_per_line'])
        nlines = int(params['num_output_lines'])
        radar_frequency = params.get('radar_frequency', None)
        heading_par = float(params.get('centre_heading', None))
        displ = num.fromfile(filename, dtype='>f4')

        # Resize array if last line is not scanned completely
        fill = 0
        if (displ.size % ncols) != 0:
            fill = num.empty(ncols - displ.size % ncols)
            fill.fill(num.nan)
            displ = num.append(displ, fill)

        displ = displ.reshape(nlines, ncols)
        displ[displ == -0.] = num.nan
        displ = num.flipud(displ)

        if radar_frequency and '_dsp_' not in par_file:
            radar_frequency = float(radar_frequency)
            radar_frequency *= 1e6  # SNAP gives MHz
            self._log.info('Scaling displacement by radar_frequency %f GHz',
                           radar_frequency/1e9)
            wavelength = util.C / radar_frequency
            displ /= -4*num.pi
            displ *= wavelength

        elif not radar_frequency and '_dsp_' not in par_file:
            self._log.warning('Could not determine radar_frequency!')
            wavelength = None

        else:
            wavelength = None

        inc_angle = self._getLOSAngles(
            filename, 'incidenceAngleFromEllipsoid.rslc').copy()
        if fill:
            inc_angle = num.append(inc_angle, fill)
        inc_angle[inc_angle == 0.] = num.nan

        phi = num.full_like(displ, (180. - heading_par))
        theta = 90. - inc_angle.reshape(displ.shape)
        theta = num.flipud(theta)

        c = self.container

        c.displacement = displ
        c.theta = theta * d2r
        c.phi = phi * d2r

        c.meta.wavelength = wavelength
        c.meta.title = params.get('PRODUCT', 'SNAP Import')
        c.meta.satellite_name = params.get('SPH_DESCRIPTOR', 'None')

        orb = params.get('PASS', None)
        c.meta.orbital_node = orb.title() if orb else None

        c.bin_file = filename
        c.par_file = par_file

        if params['map_projection'] == 'UTM':
            utm_zone = params['projection_zone']
            try:
                utm_zone_letter = utm.latitude_to_zone_letter(
                    params['center_latitude'])
            except ValueError:
                self._log.warning('Could not parse UTM Zone letter,'
                                  ' defaulting to N!')
                utm_zone_letter = 'N'

            self._log.info('Using UTM reference: Zone %d%s',
                           utm_zone, utm_zone_letter)
            c.frame.spacing = 'meter'

            dN = abs(params['post_north'])
            dE = abs(params['post_east'])

            utm_corn_e = params['corner_east']
            utm_corn_n = params['corner_north']

            utm_corn_eo = utm_corn_e + dE * displ.shape[1]
            utm_corn_no = utm_corn_n + dN * displ.shape[0]

            utm_e = num.linspace(utm_corn_e, utm_corn_eo, displ.shape[1])
            utm_n = num.linspace(utm_corn_n, utm_corn_no, displ.shape[0])

            llLat, llLon = utm.to_latlon(utm_e.min(), utm_n.min(),
                                         utm_zone, utm_zone_letter)

        else:
            self._log.info('Using Lat/Lon reference')
            c.frame.spacing = 'degree'

            if orb.lower() == 'ascending':
                llLat = params['last_near_lat']
                llLon = params['last_near_long']

            elif orb.lower() == 'descending':
                llLat = params['last_far_lat']
                llLon = params['first_near_long']

            else:
                raise AttributeError('cannot determine orbit')

            dE = abs(params['lon_pixel_res'])
            dN = abs(params['lat_pixel_res'])

        c.frame.llLat = llLat
        c.frame.llLon = llLon

        c.frame.dE = dE
        c.frame.dN = dN

        return c
