utils0.py

"""Miscellaneous utilities - independent."""
############################################################
# Program is part of MintPy                                #
# Copyright (c) 2013, Zhang Yunjun, Heresh Fattahi         #
# Author: Zhang Yunjun, Heresh Fattahi, 2013               #
############################################################
# Contents
#   InSAR
#   File Operation
#   Coordinate
#   Orbit
#   Geometry
#   Image Processing
#   User Interaction
#   Math / Statistics
# Recommend import:
#   from mintpy.utils import utils as ut


import math
import os

import h5py
import numpy as np

# global variables
SPEED_OF_LIGHT = 299792458  # m/s
EARTH_RADIUS = 6371.0088e3  # Earth radius in meters
K = 40.31                   # m^3/s^2, constant


#################################### InSAR ##########################################
def misregistration2coherence(mu):
    """Calculate the resulting coherence due to mis-registration (coregistration error).

    Reference:
      Equation (30) in Just and Bamler (1994); Equation (4.4.27) in Hanssen (2001).

    Parameters: mu  - float / np.ndarray, mis-registration in the unit of resolution cell
                      NOTE: the unit is resolution cell, NOT pixel size/spacing.
                      Applicable to both SAR range and azimuth directions.
    Returns:    coh - float / np.ndarray, coherence.
    """

    # https://numpy.org/doc/stable/reference/generated/numpy.sinc.html
    coh = np.sinc(mu)

    # for coregistration errors >1, set coherence to zero
    if isinstance(mu, np.ndarray):
        coh[mu > 1] = 0
    else:
        if mu > 1:
            coh = 0

    return coh


def range_distance(atr, dimension=2, print_msg=True):
    """Calculate slant range distance from input attribute dict
    Parameters: atr : dict, including the following ROI_PAC attributes:
                    STARTING_RANGE
                    RANGE_PIXEL_SIZE
                    LENGTH
                    WIDTH
                dimension : int, choices = [0,1,2]
                    2 for 2d matrix, vary in range direction, constant in az direction,
                        for radar coord only
                    1 for 1d matrix, in range direction, for radar coord file
                    0 for center value
    Returns:    np.array (0, 1 or 2 D) : range distance between antenna and ground target in meters
    """
    # return center value for geocoded input file
    if 'Y_FIRST' in atr.keys() and dimension > 0:
        dimension = 0
        if print_msg:
            print('input file is geocoded, return center range distance for the whole area')

    range_n, dR = float(atr['STARTING_RANGE']), float(atr['RANGE_PIXEL_SIZE'])
    length, width = int(atr['LENGTH']), int(atr['WIDTH'])

    range_f = range_n + dR*(width-1)
    range_c = (range_f + range_n)/2.0
    if print_msg:
        print('near   range : %.2f m' % (range_n))
        print('center range : %.2f m' % (range_c))
        print('far    range : %.2f m' % (range_f))

    if dimension == 0:
        return np.array(range_c, np.float32)

    range_x = np.linspace(range_n, range_f, num=width)
    if dimension == 1:
        return np.array(range_x, np.float32)
    else:
        range_xy = np.tile(range_x, (length, 1))
        return np.array(range_xy, np.float32)


def incidence_angle(atr, dem=None, dimension=2, print_msg=True):
    """Calculate 2D matrix of incidence angle from ROI_PAC attributes, very accurate.
    Parameters: atr : dict - ROI_PAC attributes including the following items:
                     STARTING_RANGE
                     RANGE_PIXEL_SIZE
                     EARTH_RADIUS
                     HEIGHT
                     WIDTH
                     LENGTH                  #for dimension=2
                     CENTER_INCIDENCE_ANGLE  #for dimension=0
                dem : 2D array for height to calculate local incidence angle
                dimension : int,
                            2 for 2d matrix
                            1 for 1d array
                            0 for one center value
                print_msg : bool
    Returns:    inc_angle : 2D np.array, incidence angle in degree for each pixel
    Example:    dem = readfile.read('hgt.rdr')[0]
                atr = readfile.read_attribute('filt_fine.unw')
                inc_angle = ut.incidence_angle(atr, dem=dem)
    """
    vprint = print if print_msg else lambda *args, **kwargs: None

    # Return center value for geocoded input file
    if 'Y_FIRST' in atr.keys() and dimension > 0:
        dimension = 0
        vprint('input file is geocoded, return center incident angle only')

    # Check if the center inc angle already exist in the metadata
    # Notes on Mar 2024 by Alex Handwerger & Talib Oliver-Cabrera:
    # Proposing these changes after encountering a range_n value smaller than the platform height
    # for UAVSAR dataset swatch_00540, thus, the calc equation w/o considering topography won't work.
    if dimension == 0 and 'CENTER_INCIDENCE_ANGLE' in atr.keys():
        inc_angle = float(atr['CENTER_INCIDENCE_ANGLE'])
        vprint(f'center incidence angle : {inc_angle:.4f} degree (grabbed from metadata directly)')
        return inc_angle

    # Read Attributes
    range_n = float(atr['STARTING_RANGE'])
    dR = float(atr['RANGE_PIXEL_SIZE'])
    r = float(atr.get('EARTH_RADIUS', EARTH_RADIUS))
    H = float(atr['HEIGHT'])
    width = int(atr['WIDTH'])

    # Calculation
    range_f = range_n + dR * width
    inc_angle_n = (np.pi - np.arccos((r**2 + range_n**2 - (r+H)**2)/(2*r*range_n))) * 180.0/np.pi
    inc_angle_f = (np.pi - np.arccos((r**2 + range_f**2 - (r+H)**2)/(2*r*range_f))) * 180.0/np.pi
    inc_angle_c = (inc_angle_n + inc_angle_f) / 2.0
    vprint(f'near   incidence angle : {inc_angle_n:.4f} degree')
    vprint(f'center incidence angle : {inc_angle_c:.4f} degree')
    vprint(f'far    incidence angle : {inc_angle_f:.4f} degree')

    if dimension == 0:
        inc_angle = inc_angle_c

    elif dimension == 1:
        inc_angle = np.linspace(inc_angle_n, inc_angle_f, num=width,
                                endpoint='FALSE', dtype=np.float32)

    elif dimension == 2:
        length = int(atr['LENGTH'])

        # consider the local variable due to topography
        if dem is not None:
            range_dist = range_distance(atr, dimension=2, print_msg=False)
            inc_angle = (np.pi - np.arccos(((r+dem)**2 + range_dist**2 - (r+H)**2) /
                                           (2*(r+dem)*range_dist))) * 180.0/np.pi
        else:
            inc_angle = np.tile(np.linspace(inc_angle_n, inc_angle_f, num=width,
                                            endpoint='FALSE', dtype=np.float32), (length, 1))
    else:
        raise Exception(f'un-supported dimension input: {dimension}')
    return inc_angle


def incidence_angle2slant_range_distance(atr, inc_angle):
    """Calculate the corresponding slant range distance given an incidence angle

    Law of sines:
               r + H                   r               range_dist
       --------------------- = ----------------- = ------------------ = 2R
        sin(pi - inc_angle)     sin(look_angle)     sin(range_angle)

    where range_angle = inc_angle - look_angle
          R is the radius of the circumcircle.

    link: http://www.ambrsoft.com/TrigoCalc/Triangle/BasicLaw/BasicTriangle.htm

    Parameters: atr         - dict, metadata including the following items:
                                  EARTH_RADIUS
                                  HEIGHT
                inc_angle   - float / np.ndarray, incidence angle in degree
    Returns:    slant_range - float, slant range distance
    """
    if isinstance(inc_angle, str):
        inc_angle = float(inc_angle)
    inc_angle = np.array(inc_angle, dtype=np.float32) / 180 * np.pi
    r = float(atr['EARTH_RADIUS'])
    H = float(atr['HEIGHT'])

    # calculate 2R based on the law of sines
    R2 = (r + H) / np.sin(np.pi - inc_angle)

    look_angle = np.arcsin( r / R2 )
    range_angle = inc_angle - look_angle
    range_dist = R2 * np.sin(range_angle)

    return range_dist


def range_ground_resolution(atr, print_msg=False):
    """Get range resolution on the ground in meters,
        from ROI_PAC attributes, for file in radar coord
    """
    if 'X_FIRST' in atr.keys():
        print('Input file is in geo coord, no range resolution info.')
        return
    inc_angle = incidence_angle(atr, dimension=0, print_msg=print_msg)
    rg_step = float(atr['RANGE_PIXEL_SIZE'])/np.sin(inc_angle/180.0*np.pi)
    return rg_step


def azimuth_ground_resolution(atr):
    """Get azimuth resolution on the ground in meters,
        from ROI_PAC attributes, for file in radar coord
    """
    if 'X_FIRST' in atr.keys():
        print('Input file is in geo coord, no azimuth resolution info.')
        return

    proc = atr.get('PROCESSOR', 'isce')
    if proc in ['roipac', 'isce']:
        Re = float(atr['EARTH_RADIUS'])
        height = float(atr['HEIGHT'])
        az_step = float(atr['AZIMUTH_PIXEL_SIZE']) * Re / (Re + height)
    elif proc == 'gamma':
        az_step = float(atr['AZIMUTH_PIXEL_SIZE'])

    return az_step



#################################### File Operation ##########################################
def touch(fname_list, times=None):
    """python equivalent function to Unix utility - touch
    It sets the modification and access times of files to the current time of day.
    If the file doesn't exist, it is created with default permissions.
    Inputs/Output:
        fname_list - string / list of string
    """
    if not fname_list:
        return None

    if isinstance(fname_list, str):
        fname_list = [fname_list]

    fname_list = [x for x in fname_list if x is not None]
    for fname in fname_list:
        if os.path.isfile(fname):
            with open(fname, 'a'):
                os.utime(fname, times)
                print('touch '+fname)

    if len(fname_list) == 1:
        fname_list = fname_list[0]
    return fname_list



################################## Coordinate ##########################################
def utm_zone2epsg_code(utm_zone):
    """Convert UTM Zone string to EPSG code.

    Reference:
        https://docs.up42.com/data/reference/utm#utm-wgs84
        https://pyproj4.github.io/pyproj/stable/examples.html#initializing-crs

    Parameters: utm_zone  - str, atr['UTM_ZONE'], comprises a zone number
                            and a hemisphere, e.g. 11N, 36S, etc.
    Returns:    epsg_code - str, EPSG code
    Examples:   epsg_code = utm_zone2epsg_code('11N')
    """
    from pyproj import CRS
    crs = CRS.from_dict({
        'proj': 'utm',
        'zone': int(utm_zone[:-1]),
        'south': utm_zone[-1].upper() == 'S',
    })
    epsg_code = crs.to_authority()[1]
    return epsg_code


def epsg_code2utm_zone(epsg_code):
    """Convert EPSG code to UTM Zone string.

    Reference:
        https://docs.up42.com/data/reference/utm#utm-wgs84
        https://pyproj4.github.io/pyproj/stable/examples.html#initializing-crs

    Parameters: epsg_code - str / int, EPSG code
    Returns:    utm_zone  - str, atr['UTM_ZONE'], comprises a zone number
                            and a hemisphere, e.g. 11N, 36S, etc. None for
                            a EPSG code not in a UTM coordnate system
    Examples:   utm_zone = epsg_code2utm_zone('32736')
    """
    from pyproj import CRS
    crs = CRS.from_epsg(epsg_code)
    utm_zone = crs.utm_zone
    if not utm_zone:
        print(f'WARNING: input EPSG code ({epsg_code}) is NOT a UTM zone, return None and continue.')
    return utm_zone


def to_latlon(infile, x, y):
    """Convert x, y in the projection coordinates of the file to lat/lon in degree.

    Similar functionality also exists in utm.to_latlon() at:
        https://github.com/Turbo87/utm#utm-to-latitudelongitude

    Parameters: infile - str, GDAL supported file path
                x/y    - scalar or 1/2D np.ndarray, coordinates in x and y direction
    Returns:    y/x    - scalar or 1/2D np.ndarray, coordinates in latitutde and longitude
    """
    from osgeo import gdal
    from pyproj import Proj, Transformer

    # read projection info using gdal
    ds = gdal.Open(infile)
    srs = ds.GetSpatialRef()

    # if input file is already in lat/lon, do nothing and return
    if (not srs.IsProjected()) and (srs.GetAttrValue('unit') == 'degree'):
        return y, x

    # convert coordinates using pyproj
    # note that Transform.from_proj(x, y, always_xy=True) convert the x, y to lon, lat
    p_in = Proj(ds.GetProjection())
    p_out = Proj('epsg:4326')
    transformer = Transformer.from_proj(p_in, p_out)
    y, x = transformer.transform(x, y)
    return y, x


def utm2latlon(meta, easting, northing):
    """Convert UTM easting/northing in meters to lat/lon in degrees.

    Parameters: meta     - dict, mintpy attributes that includes:
                           UTM_ZONE
                easting  - scalar/list/tuple/1-2D np.ndarray, UTM    coordinates in x direction
                northing - scalar/list/tuple/1-2D np.ndarray, UTM    coordinates in y direction
    Returns:    lat      - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in y direction
                lon      - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in x direction
    """
    import utm
    zone_num = int(meta['UTM_ZONE'][:-1])
    northern = meta['UTM_ZONE'][-1].upper() == 'N'
    # set 'strict=False' to allow coordinates outside the range of a typical single UTM zone,
    # which can be common for large area analysis, e.g. the Norwegian mapping authority
    # publishes a height data in UTM zone 33 coordinates for the whole country, even though
    # most of it is technically outside zone 33.
    lat, lon = utm.to_latlon(np.array(easting), np.array(northing), zone_num,
                             northern=northern, strict=False)

    # output format
    if any(isinstance(x, (list, tuple)) for x in [easting, northing]):
        lat = lat.tolist()
        lon = lon.tolist()

    return lat, lon


def latlon2utm(meta, lat, lon):
    """Convert latitude/longitude in degrees to UTM easting/northing in meters.

    Parameters: meta     - dict, mintpy attributes that includes:
                           UTM_ZONE
                lat      - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in y direction
                lon      - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in x direction
    Returns:    easting  - scalar/list/tuple/1-2D np.ndarray, UTM    coordinates in x direction
                northing - scalar/list/tuple/1-2D np.ndarray, UTM    coordinates in y direction
    """
    import utm

    # invoke zone_num to ensure all coordinates are converted into the same single UTM zone,
    # even if they cross a UTM boundary.
    zone_num = int(meta['UTM_ZONE'][:-1])
    easting, northing = utm.from_latlon(np.array(lat), np.array(lon), force_zone_number=zone_num)[:2]

    # output format
    if any(isinstance(x, (list, tuple)) for x in [lat, lon]):
        easting = easting.tolist()
        northing = northing.tolist()

    return northing, easting


def snwe_to_wkt_polygon(snwe):
    """Convert the input bounding box in SNWE into WKT format POLYGON.

    Parameters: snwe    - list of 4 float, south, north, west and east in degrees/meters
    Returns:    polygon - str, WKT format POLYGON
    """
    S, N, W, E = snwe
    lats = [N, N, S, S, N]
    lons = [W, E, E, W, W]
    polygon = "POLYGON((" + ",".join([f"{lon} {lat}" for lon, lat in zip(lons, lats)])  + "))"
    return polygon


def get_lat_lon(meta, geom_file=None, box=None, dimension=2, ystep=1, xstep=1):
    """Extract precise pixel-wise lat/lon.

    For meta dict in geo-coordinates OR geom_file with latitude/longitude dataset

    Returned lat/lon are corresponds to the pixel center

    Parameters: meta      - dict, including LENGTH, WIDTH and Y/X_FIRST/STEP
                box       - 4-tuple of int for (x0, y0, x1, y1)
                dimension - int, output lat/lon matrix dimension, 1 or 2
                y/xstep   - int, number of pixels to skip for each output pixel
    Returns:    lats      - 1/2D np.array for latitude  in size of (length, _width_)
                lons      - 1/2D np.array for longitude in size of (_length_, width)
    """
    length, width = int(meta['LENGTH']), int(meta['WIDTH'])
    if box is None:
        box = (0, 0, width, length)

    ds_list = []
    if geom_file is not None:
        with h5py.File(geom_file, 'r') as f:
            ds_list = list(f.keys())

    if 'latitude' in ds_list:
        # read 2D matrices from geometry file
        with h5py.File(geom_file, 'r') as f:
            lats = f['latitude'][box[1]:box[3], box[0]:box[2]]
            lons = f['longitude'][box[1]:box[3], box[0]:box[2]]

    elif 'Y_FIRST' in meta.keys():
        # get lat/lon0/1
        lat_step = float(meta['Y_STEP'])
        lon_step = float(meta['X_STEP'])
        lat0 = float(meta['Y_FIRST']) + lat_step * (box[1] + 0.5)
        lon0 = float(meta['X_FIRST']) + lon_step * (box[0] + 0.5)

        lat_num = box[3] - box[1]
        lon_num = box[2] - box[0]

        lat1 = lat0 + lat_step * (lat_num - 1)
        lon1 = lon0 + lon_step * (lon_num - 1)

        # get matrix of lat/lon
        if dimension == 2:
            lats, lons = np.mgrid[lat0:lat1:lat_num*1j,
                                  lon0:lon1:lon_num*1j]

        elif dimension == 1:
            lats = np.linspace(lat0, lat1, num=lat_num, endpoint=True)
            lons = np.linspace(lon0, lon1, num=lon_num, endpoint=True)

        else:
            raise ValueError(f'un-supported dimension = {dimension}')

        # UTM to lat/lon
        if not meta['Y_UNIT'].startswith('deg') and 'UTM_ZONE' in meta.keys():
            print('UTM coordinates detected, convert UTM into lat/lon')
            lats, lons = utm2latlon(meta, easting=lons, northing=lats)

    else:
        msg = 'Can not get pixel-wise lat/lon!'
        msg += '\nmeta dict is not geocoded and/or geometry file does not contains latitude/longitude dataset.'
        raise ValueError(msg)

    lats = np.array(lats, dtype=np.float32)
    lons = np.array(lons, dtype=np.float32)

    # y/xstep
    if ystep * xstep > 1:
        if lats.ndim == 1:
            lats = lats[::ystep]
            lons = lons[::xstep]
        elif lats.ndim == 2:
            lats = lats[::ystep, ::xstep]
            lons = lons[::ystep, ::xstep]

    return lats, lons


def get_lat_lon_rdc(meta):
    """Get 2D array of lat and lon for metadata dict in radar-coord.

    WARNING: This is a rough lat/lon value, NOT accurate!

    Parameters: meta : dict, including LENGTH, WIDTH and LAT/LON_REF1/2/3/4
    Returns:    lats : 2D np.array for latitude  in size of (length, width)
                lons : 2D np.array for longitude in size of (length, width)
    """
    if 'Y_FIRST' in meta.keys():
        raise Exception('Input file is in geo-coordinates, use more accurate get_lat_lon() instead.')

    length, width = int(meta['LENGTH']), int(meta['WIDTH'])
    lats = [float(meta[f'LAT_REF{i}']) for i in [1,2,3,4]]
    lons = [float(meta[f'LON_REF{i}']) for i in [1,2,3,4]]

    lat = np.zeros((length,width),dtype = np.float32)
    lon = np.zeros((length,width),dtype = np.float32)

    for i in range(length):
        for j in range(width):
            lat[i,j] = lats[0] + j*(lats[1] - lats[0])/width + i*(lats[2] - lats[0])/length
            lon[i,j] = lons[0] + j*(lons[1] - lons[0])/width + i*(lons[2] - lons[0])/length
    return lat, lon


def four_corners(atr):
    """Get the 4 corners coordinates from metadata dict in geo-coordinates.
    Parameters: atr - dict
    Returns:    south, north, west, east - float, in degrees or meters
    Examples:   S, N, W, E = ut.four_corners(atr)
                SNWE = ut.four_corners(atr)
    """
    width  = int(atr['WIDTH'])
    length = int(atr['LENGTH'])
    lon_step = float(atr['X_STEP'])
    lat_step = float(atr['Y_STEP'])
    west  = float(atr['X_FIRST'])
    north = float(atr['Y_FIRST'])
    south = north + lat_step * length
    east  = west  + lon_step * width
    return south, north, west, east


def get_lalo_digit4display(meta, coord_unit='degree'):
    """Get the digit of the decimal place for the lat/lon info for display (e.g., at the status bar).

    Parameters: meta       - dict, metadata for the following attributes:
                             X_STEP
                             Y_STEP
                             RANGE_PIXEL_SIZE
                             AZIMUTH_PIXEL_SIZE
                coord_unit - str, coordinate unit, degree or meter
    Returns:    digit      - int, the digit for the decimal places of lat/lon
    """

    if coord_unit.startswith('meter'):
        digit = 2

    else:
        geo_step_keys = ['X_STEP', 'Y_STEP']
        rdr_step_keys = ['RANGE_PIXEL_SIZE', 'AZIMUTH_PIXEL_SIZE']

        # get step size in degree
        if all(x in meta.keys() for x in geo_step_keys):
            min_step = min([abs(float(meta[x])) for x in geo_step_keys])

        elif all(x in meta.keys() for x in rdr_step_keys):
            # default scaling for spaceborne system to ground range / azimuth
            rg_pix_size = float(meta['RANGE_PIXEL_SIZE']) / np.cos(np.deg2rad(30))
            az_pix_size = float(meta['AZIMUTH_PIXEL_SIZE']) * 0.9
            min_step = min([rg_pix_size, az_pix_size]) / 108e3
        else:
            # default pixel size 30 m
            min_step = 30 / 108e3

        # set the decimal place one order   (precision for step-range)
        if   min_step >= 2e-2:  digit = 3  # 110  m for >2.2 km
        elif min_step >= 2e-3:  digit = 4  # 11   m for >220 m
        elif min_step >= 2e-4:  digit = 5  # 1    m for >22  m
        elif min_step >= 2e-5:  digit = 6  # 0.1  m for >2.2 m
        else:                   digit = 7  # 0.01 m for <2.2 m

    return digit


###################################### Orbit ###########################################
def xyz_to_local_radius(xyz):
    """Calculate satellite height and ellipsoid local radius from orbital state vector.

    This is a simplified version of the following functions from ISCE-2:
    + isce.isceobj.Planet.xyz_to_llh()
    + isce.isceobj.Ellipsoid.localRadius()

    Parameters: xyz    - tuple of 3 float, orbital state vector
    Reference:  height - float, satellite altitude in m
                radius - float, Earth radius in m
    """

    # parameters from isce.isceobj.Planet.AstronomicalHandbook
    a = 6378137.000       # WGS84 semimajor
    e2 = 0.0066943799901  # WGS84 eccentricity squared

    # xyz --> llh
    a2 = a**2
    e4 = e2**2
    r_llh =  [0]*3
    d_llh = [[0]*3 for i in range(len(xyz))]
    xy2 = xyz[0]**2 + xyz[1]**2
    p = (xy2) / a2
    q = (1. - e2) * xyz[2]**2 / a2
    r = (p + q - e4) / 6.
    s = e4 * p * q / (4. * r**3)
    t = (1. + s + math.sqrt(s * (2. + s))) **(1. / 3.)
    u = r * (1. + t + 1./t)
    v = math.sqrt(u**2 + e4 * q)
    w = e2 * (u + v - q) / (2. * v)
    k = math.sqrt(u + v + w**2) - w
    D = k * math.sqrt(xy2) / (k + e2)
    r_llh[0] = math.atan2(xyz[2], D)
    r_llh[1] = math.atan2(xyz[1], xyz[0])
    r_llh[2] = (k + e2 - 1.) * math.sqrt(D**2 + xyz[2]**2) / k
    d_llh[0] = math.degrees(r_llh[0])
    d_llh[1] = math.degrees(r_llh[1])
    d_llh[2] = r_llh[2]
    height = r_llh[2]

    # calculate local radius
    b = a * (1.0 - e2)**0.5
    latg = math.atan( math.tan(math.radians(d_llh[0])) * a**2 / b**2 )
    arg = (math.cos(latg)**2 / a**2) + (math.sin(latg)**2 / b**2)
    radius = 1.0 / math.sqrt(arg)

    return height, radius



#################################### Geometry ##########################################
# Definition of angles:
# (los_)inc_angle - the incidence angle of the LOS vector (from the ground to the SAR platform)
#                   measured from vertical. Used in isce2.
# (los_)az_angle  - the azimuth   angle of the LOS vector (from the ground to the SAR platform)
#                   measured from the north, with anti-clockwise as positive. Used in isce2.
# orb_az_angle    - the azimuth   angle of the SAR platform's orbit (along-track direction)
#                   measured from the north, with anti-clockwise as positive
# head_angle      - the azimuth   angle of the SAR platform's orbit (along-track direction)
#                   measured from the north, with      clockwise as positive. Used in roipac.
#
# Typical values in deg for satellites with near-polar orbit:
#     AlosA: los_inc_angle = 34, los_az_angle =  102, orb_az_angle =  12, head_angle =  -12
#     AlosD: los_inc_angle = 34, los_az_angle = -102, orb_az_angle = 168, head_angle = -168
#     SenA : los_inc_angle = 40, los_az_angle =  102, orb_az_angle =  12, head_angle =  -12
#     SenD : los_inc_angle = 40, los_az_angle = -102, orb_az_angle = 168, head_angle = -168
#     NiA  : los_inc_angle = 42, los_az_angle =  -78, orb_az_angle =  12, head_angle =  -12
#     NiD  : los_inc_angle = 42, los_az_angle =   78, orb_az_angle = 168, head_angle = -168

def los2orbit_azimuth_angle(los_az_angle, look_direction='right'):
    """Convert the azimuth angle of the LOS vector to the one of the orbit flight vector.
    Parameters: los_az_angle - np.ndarray or float, azimuth angle of the LOS vector from the ground to the SAR platform
                               measured from the north with anti-clockwise direction as positive, in the unit of degrees
    Returns:    orb_az_angle - np.ndarray or float, azimuth angle of the SAR platform along track/orbit direction
                               measured from the north with anti-clockwise direction as positive, in the unit of degrees
    """
    if look_direction == 'right':
        orb_az_angle = los_az_angle - 90
    else:
        orb_az_angle = los_az_angle + 90
    orb_az_angle -= np.round(orb_az_angle / 360.) * 360.
    return orb_az_angle


def azimuth2heading_angle(az_angle, look_direction='right'):
    """Convert azimuth angle from ISCE los.rdr band2 into satellite orbit heading angle

    ISCE-2 los.* file band2 is azimuth angle of LOS vector from ground target to the satellite
        measured from the north in anti-clockwise as positive

    Below are typical values in deg for satellites with near-polar orbit:
        ascending  orbit: heading angle of -12  and azimuth angle of 102
        descending orbit: heading angle of -168 and azimuth angle of -102
    """
    if look_direction == 'right':
        head_angle = (az_angle - 90) * -1
    else:
        head_angle = (az_angle + 90) * -1
    head_angle -= np.round(head_angle / 360.) * 360.
    return head_angle


def heading2azimuth_angle(head_angle, look_direction='right'):
    """Convert satellite orbit heading angle into azimuth angle as defined in ISCE-2."""
    if look_direction == 'right':
        az_angle = (head_angle - 90) * -1
    else:
        az_angle = (head_angle + 90) * -1
    az_angle -= np.round(az_angle / 360.) * 360.
    return az_angle


def enu2los(v_e, v_n, v_u, inc_angle, head_angle=None, az_angle=None):
    """Project east/north/up motion into the line-of-sight (LOS) direction defined by incidence/azimuth angle.
    Parameters: v_e        - np.ndarray or float, displacement in east-west   direction, east  as positive
                v_n        - np.ndarray or float, displacement in north-south direction, north as positive
                v_u        - np.ndarray or float, displacement in vertical    direction, up    as positive
                inc_angle  - np.ndarray or float, incidence angle from vertical, in the unit of degrees
                head_angle - np.ndarray or float, azimuth angle of the SAR platform along track direction
                             measured from the north with clockwise direction as positive, in the unit of degrees
                az_angle   - np.ndarray or float, azimuth angle of the LOS vector from the ground to the SAR platform
                             measured from the north with anti-clockwise direction as positive, in the unit of degrees
                             head_angle = 90 - az_angle
    Returns:    v_los      - np.ndarray or float, displacement in LOS direction, motion toward satellite as positive
    """
    # unite (los_)head/az_angle into (los_)az_angle
    if az_angle is None:
        if head_angle is not None:
            az_angle = heading2azimuth_angle(head_angle)
        else:
            raise ValueError(f'invalid az_angle: {az_angle}!')

    # project ENU onto LOS
    v_los = (  v_e * np.sin(np.deg2rad(inc_angle)) * np.sin(np.deg2rad(az_angle)) * -1
             + v_n * np.sin(np.deg2rad(inc_angle)) * np.cos(np.deg2rad(az_angle))
             + v_u * np.cos(np.deg2rad(inc_angle)))

    return v_los


def en2az(v_e, v_n, orb_az_angle):
    """Project east/north motion into the radar azimuth direction.
    Parameters: v_e          - np.ndarray or float, displacement in east-west   direction, east  as positive
                v_n          - np.ndarray or float, displacement in north-south direction, north as positive
                orb_az_angle - np.ndarray or float, azimuth angle of the SAR platform along track/orbit direction
                               measured from the north with anti-clockwise direction as positive, in the unit of degrees
                               orb_az_angle = los_az_angle - 90 for right-looking radar.
    Returns:    v_az         - np.ndarray or float, displacement in azimuth direction,
                               motion along flight direction as positive
    """
    # project EN onto azimuth
    v_az = (  v_e * np.sin(np.deg2rad(orb_az_angle)) * -1
            + v_n * np.cos(np.deg2rad(orb_az_angle)))
    return v_az


def calc_azimuth_from_east_north_obs(east, north):
    """Calculate the azimuth angle of the given horizontal observation (in east and north)

    Parameters: east     - float,  eastward motion
                north    - float, northward motion
    Returns:    az_angle - float, azimuth angle in degree
                           measured from the north with anti-clockwise as positive
    """
    az_angle = -1 * np.rad2deg(np.arctan2(east, north)) % 360
    return az_angle


def get_unit_vector4component_of_interest(los_inc_angle, los_az_angle, comp='enu2los', horz_az_angle=None):
    """Get the unit vector for the component of interest.
    Parameters: los_inc_angle - np.ndarray or float, incidence angle from vertical, in the unit of degrees
                los_az_angle  - np.ndarray or float, azimuth angle of the LOS vector from the ground to the SAR platform
                                measured from the north with anti-clockwise direction as positive, in the unit of degrees
                comp          - str, component of interest, choose among the following values:
                                enu2los, en2los, hz2los, u2los, up2los, orb(it)_az, vert, horz
                horz_az_angle - np.ndarray or float, azimuth angle of the horizontal direction of interest
                                measured from the north with anti-clockwise direction as positive, in the unit of degrees
    Returns:    unit_vec      - list(np.ndarray/float), unit vector of the ENU component for the component of interest
    """
    # check input arguments
    comps = [
        'enu2los', 'en2los', 'hz2los', 'horz2los', 'u2los', 'vert2los',   # radar LOS / cross-track
        'en2az', 'hz2az', 'orb_az', 'orbit_az',                           # radar azimuth / along-track
        'vert', 'vertical', 'horz', 'horizontal',                         # vertical / arbitrary horizontal
    ]

    if comp not in comps:
        raise ValueError(f'un-recognized comp input: {comp}.\nchoose from: {comps}')

    if comp == 'horz' and horz_az_angle is None:
        raise ValueError('comp=horz requires horz_az_angle input!')

    # initiate output
    unit_vec = None

    if comp in ['enu2los']:
        unit_vec = [
            np.sin(np.deg2rad(los_inc_angle)) * np.sin(np.deg2rad(los_az_angle)) * -1,
            np.sin(np.deg2rad(los_inc_angle)) * np.cos(np.deg2rad(los_az_angle)),
            np.cos(np.deg2rad(los_inc_angle)),
        ]

    elif comp in ['en2los', 'hz2los', 'horz2los']:
        unit_vec = [
            np.sin(np.deg2rad(los_inc_angle)) * np.sin(np.deg2rad(los_az_angle)) * -1,
            np.sin(np.deg2rad(los_inc_angle)) * np.cos(np.deg2rad(los_az_angle)),
            np.zeros_like(los_inc_angle),
        ]

    elif comp in ['u2los', 'vert2los']:
        unit_vec = [
            np.zeros_like(los_inc_angle),
            np.zeros_like(los_inc_angle),
            np.cos(np.deg2rad(los_inc_angle)),
        ]

    elif comp in ['en2az', 'hz2az', 'orb_az', 'orbit_az']:
        orb_az_angle = los2orbit_azimuth_angle(los_az_angle)
        unit_vec = [
            np.sin(np.deg2rad(orb_az_angle)) * -1,
            np.cos(np.deg2rad(orb_az_angle)),
            np.zeros_like(orb_az_angle),
        ]

    elif comp in ['vert', 'vertical']:
        unit_vec = [0, 0, 1]

    elif comp in ['horz', 'horizontal']:
        unit_vec = [
            np.sin(np.deg2rad(horz_az_angle)) * -1,
            np.cos(np.deg2rad(horz_az_angle)),
            np.zeros_like(horz_az_angle),
        ]

    return unit_vec



#################################### Image Processing ##########################################
def wrap(data_in, wrap_range=[-1.*np.pi, np.pi]):
    """Wrap data into a range.
    Parameters: data_in    : np.array, array to be wrapped
                wrap_range : list of 2 float, range to be wrapped into
    Returns:    data       : np.array, data after wrapping
    """
    w0, w1 = wrap_range
    data = np.array(data_in)
    data = w0 + np.mod(data - w0, w1 - w0)
    return data


def diff_wrapped_phase(pha1, pha2):
    """Calculate the difference between two input wrapped phase.

    Parameters: pha1     - np.ndarray, (un)wrapped phase
                pha2     - np.ndarray / float, (un)wrapped phase
    Returns:    pha_diff - np.ndarray, wrapped phase difference for (pha1 - pha2)
    """
    pha_diff = np.angle(np.exp(-1j * pha1) * np.conj(np.exp(-1j * pha2)))
    return pha_diff


def get_all_conn_components(mask_in, min_num_pixel=1e4):
    """Get all connected component with number of pixels larger than threshold
    Parameters: mask_in  : 2D np.array with zero as background and non-zero as foreground
                min_num_pixel : int/float, minimum number of pixels to be identified and output
    Returns:    mask_out : list of 2D np.array in np.bool_ format
    """
    mask_in = np.array(mask_in)
    mask_out = []  # list of 2D np.array in bool
    mask_cc = get_largest_conn_component(mask_in, min_num_pixel=1e4)
    while not np.all(~mask_cc):
        mask_out.append(mask_cc)
        mask_in ^= mask_cc
        mask_cc = get_largest_conn_component(mask_in, min_num_pixel=1e4)
    return mask_out


def get_largest_conn_component(mask_in, min_num_pixel=1e4, display=False):
    """Extract the largest connected component from an 2D array
       with zero as background value
    Parameters: mask_in  : 2D np.array with zero as background and non-zero as foreground
                min_num_pixel : int/float, minimum number of pixels to be identified and output
                display : bool, display the result or not.
    Returns:    mask_out : 2D np.array in np.bool_ format
    """
    import matplotlib.pyplot as plt
    from scipy import ndimage

    mask_out = np.zeros(mask_in.shape, np.bool_)
    labels = ndimage.label(mask_in)[0]
    num_pixel = np.max(np.bincount(labels.flatten())[1:])
    if num_pixel < min_num_pixel:
        return mask_out

    max_label = np.argmax(np.bincount(labels.flatten())[1:]) + 1
    mask_out = labels == max_label
    if display:
        _, ax = plt.subplots(nrows=1, ncols=3, figsize=[15, 5])
        ax[0].imshow(mask_in)
        ax[1].imshow(mask_out)
        ax[2].imshow(mask_in ^ mask_out)
        plt.show()
    return mask_out


def min_region_distance(mask1, mask2, display=False):
    """Calculate the min distance between two regions of pixels marked by mask1 and mask2
    Parameters: mask1/2 : 2D np.array in size of (length, width) in np.bool_ format
    Returns:    pts1 : tuple of 2 int, bridge point in mask1, in (x, y)
                pts2 : tuple of 2 int, bridge point in mask2, in (x, y)
                min_dist : float, min euclidean distance
    """
    import matplotlib.pyplot as plt
    from scipy.spatial import cKDTree

    y, x = np.where(mask1 != 0)
    pts1 = np.hstack((x.reshape(-1, 1), y.reshape(-1, 1)))
    tree = cKDTree(pts1)

    y, x = np.where(mask2 != 0)
    pts2 = np.hstack((x.reshape(-1, 1), y.reshape(-1, 1)))
    dist, idx = tree.query(pts2)

    idx_min = np.argmin(dist)
    xy2 = pts2[idx_min]
    xy1 = pts1[idx[idx_min]]
    min_dist = dist[idx_min]

    if display:
        plt.figure()
        plt.imshow(mask1 * 1 + mask2 * 2)
        plt.plot([xy1[0], xy2[0]], [xy1[1], xy2[1]], '-o')
        plt.show()

    return xy1, xy2, min_dist


def interpolate_data(in_data, out_shape, interp_method='linear'):
    """Interpolate input 2D matrix into different shape.

    Used to get full resolution perp baseline from ISCE coarse grid baseline file.

    Parameters: in_data       : 2D np.ndarray
                out_shape     : tuple of 2 int in (length, width)
                interp_method : string, choose in [nearest, linear, cubic]
    Returns:    out_data      : 2D np.ndarray in out_shape
    """
    from scipy.interpolate import RegularGridInterpolator

    # prepare interpolation function
    in_shape = in_data.shape
    in_pts = (np.arange(in_shape[0]), np.arange(in_shape[1]))
    interp_func = RegularGridInterpolator(
        in_pts,
        in_data,
        method=interp_method,
        bounds_error=False,
    )

    # prepare output coordinates
    xx, yy = np.meshgrid(np.linspace(0, in_shape[1]-1, out_shape[1], endpoint=False),
                         np.linspace(0, in_shape[0]-1, out_shape[0], endpoint=False))
    out_pts = np.hstack((yy.reshape(-1, 1), xx.reshape(-1, 1)))

    # run interpolation
    out_data = interp_func(out_pts).reshape(out_shape)

    return out_data


def polygon2mask(polygon, shape):
    """Create a 2D mask (numpy array in binary) from a polygon.

    Link: https://stackoverflow.com/questions/3654289/scipy-create-2d-polygon-mask
    Parameters: polygon - list of tuples of 2 int, e.g. [(x1, y1), (x2, y2), ...]
                shape   - list/tuple of 2 int, for length and width
    Returns:    mask    - 2D np.ndarray in bool in size of (length, width)
    """
    from PIL import Image, ImageDraw
    length, width = shape

    img = Image.new('L', (width, length), 0)
    ImageDraw.Draw(img).polygon(polygon, outline=1, fill=1)
    mask = np.array(img, dtype=np.bool_)

    return mask


def get_circular_mask(x, y, radius, shape):
    """Get mask of pixels within circle defined by (x, y, r)"""
    length, width = shape
    yy, xx = np.ogrid[-y:length-y,
                      -x:width-x]
    cmask = (xx**2 + yy**2 <= radius**2)
    return cmask


def circle_index(atr, circle_par):
    """Return Index of Elements within a Circle centered at input pixel
    Parameters: atr : dictionary
                    containing the following attributes:
                    WIDT
                    LENGTH
                circle_par : string in the format of 'y,x,radius'
                    i.e. '200,300,20'          for radar coord
                         '31.0214,130.5699,20' for geo   coord
    Returns:    idx : 2D np.array in bool type
                    mask matrix for those pixel falling into the circle
                    defined by circle_par
    Examples:   idx_mat = ut.circle_index(atr, '200,300,20')
                idx_mat = ut.circle_index(atr, '31.0214,130.5699,20')
    """

    width = int(atr['WIDTH'])
    length = int(atr['LENGTH'])

    if isinstance(circle_par, tuple):
        cir_par = circle_par
    elif isinstance(circle_par, list):
        cir_par = circle_par
    else:
        cir_par = circle_par.replace(',', ' ').split()
    cir_par = [str(i) for i in cir_par]

    try:
        c_y = int(cir_par[0])
        c_x = int(cir_par[1])
        radius = int(float(cir_par[2]))
        print('Input circle index in y/x coord: %d, %d, %d' % (c_y, c_x, radius))
    except:
        try:
            c_lat = float(cir_par[0])
            c_lon = float(cir_par[1])
            radius = int(float(cir_par[2]))
            c_y = np.rint((c_lat-float(atr['Y_FIRST'])) / float(atr['Y_STEP']))
            c_x = np.rint((c_lon-float(atr['X_FIRST'])) / float(atr['X_STEP']))
            print('Input circle index in lat/lon coord: '
                  '{:.4f}, {:.4f}, {}'.format(c_lat, c_lon, radius))
        except:
            print('\nERROR: Unrecognized circle index format: '+circle_par)
            print('Supported format:')
            print('--circle 200,300,20            for radar coord input')
            print('--circle 31.0214,130.5699,20   for geo   coord input\n')
            return 0

    y, x = np.ogrid[-c_y:length-c_y, -c_x:width-c_x]
    idx = x**2 + y**2 <= radius**2

    return idx



#################################### User Interaction #####################################
def yes_or_no(question):
    """garrettdreyfus on Github: https://gist.github.com/garrettdreyfus/8153571"""
    reply = str(input(question+' (y/n): ')).lower().strip()
    if reply[0] == 'y':
        return True
    elif reply[0] == 'n':
        return False
    else:
        return yes_or_no("Uhhhh... please enter ")


def update_attribute_or_not(atr_new, atr_orig):
    """Compare new attributes with existing ones"""
    for key in atr_new.keys():
        value = str(atr_new[key])
        if not ((key in atr_orig.keys() and value == str(atr_orig[key]) and value != 'None')
                or (key not in atr_orig.keys() and value == 'None')):
            return True
    return False


def which(program):
    """Test if executable exists"""
    def is_exe(fpath):
        return os.path.isfile(fpath) and os.access(fpath, os.X_OK)

    if os.path.split(program)[0]:
        if is_exe(program):
            return program
    else:
        for path in os.environ["PATH"].split(os.pathsep):
            path = path.strip('"')
            exe_file = os.path.join(path, program)
            if is_exe(exe_file):
                return exe_file
    return None


def check_parallel(file_num=1, print_msg=True, maxParallelNum=8):
    """Check parallel option based file num and installed module
    Link: https://joblib.readthedocs.io/en/latest/parallel.html
    Examples:
        num_cores, inps.parallel, Parallel, delayed = ut.check_parallel(len(file_list))
        Parallel(n_jobs=num_cores)(delayed(subset_file)(file, vars(inps)) for file in file_list)
    """
    enable_parallel = True

    # Disable parallel option for one input file
    if file_num <= 1:
        enable_parallel = False
        if print_msg:
            print('parallel processing is disabled for one input file')
        return 1, enable_parallel, None, None

    # Check required python module
    try:
        from joblib import Parallel, delayed
    except:
        print('Can not import joblib')
        print('parallel is disabled.')
        enable_parallel = False
        return 1, enable_parallel, None, None

    # Find proper number of cores for parallel processing
    num_cores = min(os.cpu_count(), file_num, maxParallelNum)
    if num_cores <= 1:
        enable_parallel = False
        print('parallel processing is disabled because min of the following two numbers <= 1:')
        print(f'available cpu number of the computer: {os.cpu_count()}')
    elif print_msg:
        print('parallel processing using %d cores ...' % (num_cores))

    try:
        return num_cores, enable_parallel, Parallel, delayed
    except:
        return num_cores, enable_parallel, None, None


def print_command_line(script_name, args):
    """print the command line with "" for arguments containing *.

    Parameters: script_name - str, e.g. prep_isce.py
                args        - list(str), list of input arguments
    """
    cmd = script_name
    for arg in args:
        # for option values containing *, add "" in the print out msg
        # so that it can be copied and pasted to run directly.
        if not arg.startswith('-') and '*' in arg:
            cmd += f' "{arg}"'
        else:
            cmd += f' {arg}'

    print(cmd)
    return


#################################### Math / Statistics ###################################
def median_abs_deviation(data, center=None, scale=0.67449):
    """Compute the median absolute deviation (MAD) of the data along the LAST axis.

    This function is also available in the following packages:
    + statsmodels.robust.mad()
      https://www.statsmodels.org/dev/generated/statsmodels.robust.scale.mad.html
    + scipy.stats.median_abs_deviation() since v1.5.0
      https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.median_abs_deviation.html

    The implementation here is preferred because we would like to:
    1. omit the NaN value in the data for the median and MAD calculation
    2. scale the returned value to be comparable with standard deviation (STD)
       for easy interpretation with 1/2/3-sigma rule.
    While statsmodels could not handle the NaN value in the data; and scipy does not hebave
    as desired in the default setting (propagate nan and do not scale).

    The following two are equivalent for a 2D np.ndarray X:
    + mintpy.utils.utils0.median_abs_deviation(X)
    + scipy.stats.median_abs_deviation(X, axis=-1, center=np.nanmedian, scale=0.67449, nan_policy='omit')

    Parameters: data   - 1/2D np.ndarray, input array
                center - 0/1D np.ndarray or None
                scale  - float, the normalization constant
    Returns:    mad    - 0/1D np.ndarray
    """
    # check input: compute MAD over 1/2D matrix only
    data = np.array(data)
    if data.ndim > 2:
        ntime = data.shape[0]
        data = data.reshape(ntime, -1)

    # calculate the center value: median along the last axis
    if center is None:
        center = np.nanmedian(data, axis=-1)

        # replicate center matrix for 2D matrix calculation later on
        if data.ndim == 2:
            center = np.tile(center.reshape(-1, 1), (1, data.shape[1]))

    # calculate median absolute deviation
    mad = np.nanmedian(np.abs(data - center), axis=-1) / scale

    return mad


def median_abs_deviation_threshold(data, center=None, cutoff=3.):
    """calculate rms_threshold based on the standardized residual

    Outlier detection with median absolute deviation.
    """
    if center is None:
        center = np.nanmedian(data)
    mad = median_abs_deviation(data, center=center)
    threshold = center + cutoff * mad
    return threshold


def root_mean_sq_error(x, y=None):
    """Calculate the root-mean-square error between x and y."""
    # make a copy & ensure 1D numpy array format
    x = np.array(x).flatten()
    if y is not None:
        y = np.array(y).flatten()
        if x.size != y.size:
            raise ValueError(f'Input x & y have different size: {x.size} vs. {y.size}!')
        x -= y

    # omit NaN values
    x = x[~np.isnan(x)]
    rmse = np.sqrt(np.sum(x**2) / (x.size - 1))
    return rmse


def ceil_to_1(x):
    """Return the most significant digit of input number and ceiling it"""
    digit = int(np.floor(np.log10(abs(x))))
    x_round = round(x, -digit)
    # round to ceil
    if x_round >= x:
        x_ceil = x_round
    else:
        x_ceil = x_round + 10**digit
    return x_ceil


def round_to_1(x):
    """Return the most significant digit of input number"""
    digit = int(np.floor(np.log10(abs(x))))
    return round(x, -digit)


def round_up_to_odd(x):
    """Round a float up to the next odd integer .
    Link: https://stackoverflow.com/questions/31648729
    """
    y = np.ceil(x) // 2 * 2 + 1
    return y.astype(np.int16)


def highest_power_of_2(x):
    """Given a number x, find the highest power of 2 that <= x"""
    res = np.power(2, np.floor(np.log2(x)))
    res = np.int16(res)
    return res


def most_common(L, k=1):
    """Return the k most common item in the list L.
    Examples:
        5, 8 = most_common([4,5,5,5,5,8,8,8,9], k=2)
        'duck' = most_common(['goose','duck','duck','dog'])
        'goose' = most_common(['goose','duck','duck','goose'])
    """
    from collections import Counter
    cnt = Counter(L)
    item_mm = [i[0] for i in cnt.most_common(k)]
    if k == 1:
        item_mm = item_mm[0]
    return item_mm


def is_number(string):
    """Check string is a number.
    Not using str.isnumeric() because it can not handle floating point nor negative sign.
    Link: https://elearning.wsldp.com/python3/python-check-string-is-a-number/
    Parameters: string - str, a string
    Returns:    True/False
    """
    try:
        float(string)
        return True
    except ValueError:
        return False