ISCE_INSAR/components/isceobj/Location/Offset.py

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Copyright 2010 California Institute of Technology. ALL RIGHTS RESERVED.
# 
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# 
# http://www.apache.org/licenses/LICENSE-2.0
# 
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# 
# United States Government Sponsorship acknowledged. This software is subject to
# U.S. export control laws and regulations and has been classified as 'EAR99 NLR'
# (No [Export] License Required except when exporting to an embargoed country,
# end user, or in support of a prohibited end use). By downloading this software,
# the user agrees to comply with all applicable U.S. export laws and regulations.
# The user has the responsibility to obtain export licenses, or other export
# authority as may be required before exporting this software to any 'EAR99'
# embargoed foreign country or citizen of those countries.
#
# Author: Walter Szeliga
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~




import math
from isce import logging

from isceobj.Util.decorators import type_check, force, pickled, logged
import numpy as np

from iscesys.Component.Component import Component

X = Component.Parameter('x',
        public_name='x',
        default=None,
        type=float,
        mandatory=True,
        doc = 'Range location')

DX = Component.Parameter('dx',
        public_name='dx',
        default=None,
        type=float,
        mandatory=True,
        doc = 'Range offset')

Y = Component.Parameter('y',
        public_name='y',
        default=None,
        type=float,
        mandatory=True,
        doc = 'Azimuth location')

DY = Component.Parameter('dy',
        public_name='dy',
        default=None,
        type=float,
        mandatory=True,
        doc = 'Azimuth offset')

SIGMA_X = Component.Parameter('sigmax',
        public_name='sigmax',
        default=0,
        type=float,
        mandatory=True,
        doc = 'Range covariance')

SIGMA_Y = Component.Parameter('sigmay',
        public_name='sigmay',
        default=0,
        type=float,
        mandatory=True,
        doc = 'Azimuth covariance')

SIGMA_XY = Component.Parameter('sigmaxy',
        public_name='sigmaxy',
        default=0,
        type=float,
        mandatory=True,
        doc = 'Cross covariance')

SNR = Component.Parameter('snr',
        public_name='snr',
        default=0,
        type=float,
        mandatory=True,
        doc = 'Signal to Noise Ratio')

@pickled
class Offset(Component):
    """A class to represent the two-dimensional offset of a particular
    location"""

    logging_name = "isceobj.Location.Offset"

    family = 'offset'
    parameter_list = (
                      X,
                      Y,
                      DX,
                      DY,
                      SIGMA_X,
                      SIGMA_Y,
                      SIGMA_XY,
                      SNR
                     )
    @logged
    def __init__(self, x=None, y=None, dx=None, dy=None, snr=0.0,
                 sigmax=0.0, sigmay=0.0, sigmaxy=0.0,
                 family = None, name = None):
        super(Offset, self).__init__(family=family if family else  self.__class__.family, name=name)
        self.x = x
        self.dx = dx
        self.y = y
        self.dy = dy
        self.setSignalToNoise(snr)
        self.setCovariance(sigmax, sigmay, sigmaxy)
        return None

    def setCoordinate(self, x, y):
        self.x = x
        self.y = y

    def setOffset(self, dx, dy):
        self.dx = dx
        self.dy = dy
        pass
    def setCovariance(self, covx, covy, covxy):
        self.sigmax = covx
        self.sigmay = covy
        self.sigmaxy = covxy

    @force(float)
    def setSignalToNoise(self, snr):
        self.snr = snr if not math.isnan(snr) else 0.0

    def getCoordinate(self):
        return self.x,self.y

    def getOffset(self):
        return self.dx,self.dy

    def getSignalToNoise(self):
        return self.snr

    def getCovariance(self):
        return self.sigmax, self.sigmay, self.sigmaxy


    def __str__(self):
        retstr = "%s %s %s %s %s %s %s %s" % (self.x,self.dx,self.y,self.dy,self.snr, self.sigmax, self.sigmay, self.sigmaxy)
        return retstr


OFFSETS = Component.Parameter('_offsets',
        public_name='offsets',
        default=[],
        type=float,
        mandatory=False,
        intent='output',
        doc = 'List of offsets')

@pickled
class OffsetField(Component):
    """A class to represent a collection of offsets defining an offset field"""
    logging_name = "isceobj.Location.OffsetField"

    family = 'offsetfield'
    parameter_list = (
                      OFFSETS,
                     )
    @logged
    def __init__(self,family=None,name=None):

        super(OffsetField, self).__init__(
            family=family if family else  self.__class__.family, name=name)
        self._last = 0
        self._cpOffsets = None
        return None

    #extend dump method. Convert Offset object to list before dumping it
    def dump(self,filename):

        self.adaptToRender()
        super(OffsetField,self).dump(filename)
        #restore to list of Offset
        self.restoreAfterRendering()

    def load(self,filename):
        import copy

        super(OffsetField,self).load(filename)
        #make a copy
        cpOffsets = copy.deepcopy(self._offsets)
        self.packOffsetswithCovariance(cpOffsets)

    def adaptToRender(self):
        import copy
        #make a copy before dumping
        self._cpOffsets = copy.deepcopy(self._offsets)
        #change the offsets to a list on numbers instead of Offset
        self._offsets = self.unpackOffsetswithCovariance()

    def restoreAfterRendering(self):
        self._offsets = self._cpOffsets

    def initProperties(self,catalog):
        if 'offsets' in catalog:
            offsets = catalog['offsets']
            import numpy as np
            offsets = np.array(offsets)
            self.packOffsetswithCovariance(offsets.T)
            catalog.pop('offsets')
        super().initProperties(catalog)

    def getLocationRanges(self):
        xdxydysnr = self.unpackOffsets()
        numEl = len(xdxydysnr)
        x = np.zeros(numEl)
        y = np.zeros(numEl)
        for i in range(numEl):
            x[i] =  xdxydysnr[i][0]
            y[i] =  xdxydysnr[i][2]
        xr = sorted(x)
        yr = sorted(y)
        return [xr[0],xr[-1],yr[0],yr[-1]]

    def plot(self,type,xmin = None, xmax = None, ymin = None, ymax = None):
        try:
            import numpy as np
            from  scipy.interpolate import griddata
            import matplotlib.pyplot as plt
            from pylab import quiver,quiverkey
        except ImportError:
            self.logger.error('This method requires scipy, numpy and matplotlib to be installed.')
        xdxydysnr = self.unpackOffsets()
        numEl = len(xdxydysnr)
        x = np.zeros(numEl)
        y = np.zeros(numEl)
        dx = np.zeros(numEl)
        dy = np.zeros(numEl)
        for i in range(numEl):
            x[i] =  xdxydysnr[i][0]
            dx[i] =  xdxydysnr[i][1]
            y[i] =  xdxydysnr[i][2]
            dy[i] =  xdxydysnr[i][3]
        if xmin is None: xmin = np.min(x)
        if xmax is None: xmax = np.max(x)
        if ymin is None: ymin = np.min(y)
        if ymax is None: ymax = np.max(y)
        legendL = np.floor(max(np.max(dx),np.max(dy)))
        #normally the width in range is much smaller that the length in azimuth, so normalize so that we have the same number os sample for each axis
        step = min(np.min(int(np.ceil(((ymax - ymin)/(xmax - xmin))))),5)
        X , Y = np.mgrid[xmin:xmax,ymin:ymax:step]
        skip = int(np.ceil(xmax - xmin)/100)*5
        if type == 'field':
            U = griddata(np.array([x,y]).T,dx, (X,Y), method='linear')
            V = griddata(np.array([x,y]).T,dy, (X,Y), method='linear')
            Q = quiver(X[::skip,::skip], Y[::skip,::skip],
                       U[::skip,::skip], V[::skip,::skip],
                       pivot='mid', color='g')
            arrow = str(legendL) + ' pixles'
            qk = quiverkey(Q, 0.8, 0.95, legendL, arrow,
                           labelpos='E',
                           coordinates='figure',
                           fontproperties={'weight':'bold'})
            ax = Q.axes
            ax.set_xlabel('range location')
            ax.set_ylabel('azimuth location')
        elif(type == 'pcolor'):
            M = griddata(np.array([x,y]).T,
                         np.sqrt(dx**2 + dy**2),
                         (X,Y),
                         method='linear')
            P = griddata(np.array([x,y]).T,
                         np.arctan2(dy, dx),
                         (X,Y)
                         ,method='linear')
            plt.subplot(2, 1, 1)
            plt.imshow(M,aspect='auto', extent=[xmin, xmax, ymin, ymax])
            plt.colorbar()
            ax1 = plt.gca()
            ax1.set_ylabel('azimuth location')
            ax1.set_title('offset magnitude')
            plt.subplot(2, 1, 2)
            plt.imshow(P, aspect='auto', extent=[xmin,xmax,ymin,ymax])
            plt.colorbar()
            ax2 = plt.gca()
            ax2.set_xlabel('range location')
            ax2.set_ylabel('azimuth location')
            ax2.set_title('offset phase')
        plt.show()
        return plt

    @type_check(Offset)
    def addOffset(self, offset):
        self._offsets.append(offset)
        pass

    def __next__(self):
        if self._last < len(self._offsets):
            next = self._offsets[self._last]
            self._last += 1
            return next
        else:
            self._last = 0 # This is so that we can restart iteration
            raise StopIteration()

    def packOffsets(self, offsets):#create an offset field from a list of offets
        self._offset = []
        for i in range(len(offsets[0])):
            #note that different ordering (x,y,dx,dy,snr) instead of (x,dx,y,dy,snr)
            self.addOffset(
                Offset(x=offsets[0][i],
                       y=offsets[2][i],
                       dx=offsets[1][i],
                       dy=offsets[3][i],
                       snr=offsets[4][i])
                )

    def packOffsetswithCovariance(self, offsets):
        self._offset = []
        for i in range(len(offsets[0])):
            self.addOffset(
                    Offset(x=offsets[0][i],
                           y=offsets[2][i],
                           dx=offsets[1][i],
                           dy=offsets[3][i],
                           snr=offsets[4][i],
                           sigmax=offsets[5][i],
                           sigmay=offsets[6][i],
                           sigmaxy=offsets[7][i])
                    )

    def unpackOffsets(self):
        """A convenience method for converting our iterator into a flat
        list for use in Fortran and C code"""
        offsetArray = []
        for offset in self.offsets:
            x, y = offset.getCoordinate()
            dx, dy = offset.getOffset()
            snr = offset.getSignalToNoise()
            offsetArray.append([x,dx,y,dy,snr])
            pass
        return offsetArray

    def unpackOffsetswithCovariance(self):
        offsetArray = []
        for offset in self.offsets:
            x,y = offset.getCoordinate()
            dx,dy = offset.getOffset()
            snr = offset.getSignalToNoise()
            sx, sy, sxy = offset.getCovariance()
            offsetArray.append([x,dx,y,dy,snr,sx,sy,sxy])
            pass
        return offsetArray

    def cull(self, snr=0.0):
        """Cull outliers based on their signal-to-noise ratio.

        @param snr: the signal-to-noise ratio to use in the culling.  Values with greater signal-to-noise will be kept.
        """
        culledOffsetField = OffsetField()
        i = 0
        for offset in self.offsets:
            if (offset.getSignalToNoise() < snr):
                i += 1
            else:
                culledOffsetField.addOffset(offset)

        self.logger.info("%s offsets culled" % (i))
        return culledOffsetField

    def __iter__(self):
        return self

    def __str__(self):
        return '\n'.join(map(str, self.offsets))+'\n' #2013-06-03 Kosal: wrong use of map

    @property
    def offsets(self):
        return self._offsets

    pass

    def getFitPolynomials(self,rangeOrder=2,azimuthOrder=2,maxOrder=True, usenumpy=False):
        from stdproc.stdproc.offsetpoly.Offsetpoly import Offsetpoly
        from isceobj.Util import Poly2D

        numCoeff = 0
        ####Try and use Howard's polynomial fit code whenever possible
        if (rangeOrder == azimuthOrder) and (rangeOrder <= 3):
            if (rangeOrder == 1):
                if maxOrder:
                    numCoeff = 3
                else:
                    numCoeff = 4
            elif (rangeOrder == 2):
                if maxOrder:
                    numCoeff = 6
            elif (rangeOrder == 3):
                if maxOrder:
                    numCoeff = 10


        inArr = np.array(self.unpackOffsets())
        azmin = np.min(inArr[:,2])
        inArr[:,2] -= azmin

        ####Use Howard's code
        if (numCoeff != 0) and not usenumpy:
            x = list(inArr[:,0])
            y = list(inArr[:,2])
            dx = list(inArr[:,1])
            dy = list(inArr[:,3])
            sig = list(inArr[:,4])

            ####Range Offsets
            obj = Offsetpoly()
            obj.setLocationAcross(x)
            obj.setLocationDown(y)
            obj.setSNR(sig)
            obj.setOffset(dx)
            obj.numberFitCoefficients = numCoeff
            obj.offsetpoly()

            val = obj.offsetPoly

            #####Unpack into 2D array
            if numCoeff == 3:
                coeffs = [[val[0], val[1]],
                          [val[2], 0.0]]

            elif numCoeff == 4:
                coeffs = [[val[0], val[1]],
                          [val[2], val[3]]]

            elif numCoeff == 6:
                coeffs = [[val[0], val[1], val[4]],
                          [val[2], val[3], 0.0],
                          [val[5], 0.0, 0.0]]

            elif numCoeff == 10:
                ####Unpacking needs to be checked.
                coeffs = [[val[0], val[1], val[4], val[8]],
                          [val[2], val[3], val[6], 0.0],
                          [val[5], val[7],0.0, 0.0],
                          [val[9], 0.0, 0.0, 0.0]]


            rangePoly = Poly2D.Poly2D()
            rangePoly.setMeanAzimuth(azmin)
            rangePoly.initPoly(rangeOrder=rangeOrder, azimuthOrder=azimuthOrder, coeffs=coeffs)

            ####Azimuth Offsets
            obj.setOffset(dy)
            obj.offsetpoly()
            val = obj.offsetPoly

            #####Unpack into 2D array
            if numCoeff == 3:
                coeffs = [[val[0], val[1]],
                          [val[2], 0.0]]

            elif numCoeff == 4:
                coeffs = [[val[0], val[1]],
                          [val[2], val[3]]]

            elif numCoeff == 6:
                coeffs = [[val[0], val[1], val[4]],
                          [val[2], val[3], 0.0],
                          [val[5], 0.0, 0.0]]

            elif numCoeff == 10:
                ####Unpacking needs to be checked.
                coeffs = [[val[0], val[1], val[4], val[8]],
                          [val[2], val[3], val[6], 0.0],
                          [val[5], val[7],0.0, 0.0],
                          [val[9], 0.0, 0.0, 0.0]]

            azimuthPoly = Poly2D.Poly2D()
            azimuthPoly.setMeanAzimuth(azmin)
            azimuthPoly.initPoly(rangeOrder=rangeOrder, azimuthOrder=azimuthOrder, coeffs=coeffs)

        ####Fallback to numpy based polynomial fitting
        else:

            x = inArr[:,0]
            y = inArr[:,2]
            dx = inArr[:,1]
            dy = inArr[:,3]
            sig = inArr[:,4]


            azimuthPoly = Poly2D.Poly2D()
            azimuthPoly.initPoly(rangeOrder=rangeOrder, azimuthOrder=azimuthOrder)
            azimuthPoly.polyfit(x,y,dy, sig=sig)
            azimuthPoly._meanAzimuth += azmin

            rangePoly = Poly2D.Poly2D()
            rangePoly.initPoly(rangeOrder=rangeOrder, azimuthOrder=azimuthOrder)
            rangePoly.polyfit(x,y,dx,sig=sig)
            rangePoly._meanAzimuth += azmin

        return (azimuthPoly, rangePoly)