ISCE_INSAR/components/isceobj/TopsProc/runIon.py

#
# Author: Cunren Liang
# Copyright 2018
# California Institute of Technology
#

import os
import shutil
import datetime
import numpy as np
import numpy.matlib

import isceobj
import logging
from isceobj.Constants import SPEED_OF_LIGHT
from isceobj.TopsProc.runBurstIfg import loadVirtualArray
from isceobj.Alos2Proc.runIonFilt import reformatMaskedAreas


logger = logging.getLogger('isce.topsinsar.ion')

#should get rid of the coherence thresholds in the future
##WARNING: when using the original full-bandwidth swath xml file, should also consider burst.image.filename
class dummy(object):
    pass


def setup(self):
    '''
    setup parameters for processing
    '''

    #initialize parameters for ionospheric correction
    ionParam = dummy()
    #The step names in the list below are exactly the function names in 'def runIon(self):'
    #when adding a new step, only put its function name (in right order) in the list,
    #and put the function (in right order) in 'def runIon(self):'
    ionParam.allSteps = ['subband', 'rawion', 'grd2ion', 'filt_gaussian', 'ionosphere_shift', 'ion2grd', 'esd']


    ###################################################################
    #users are supposed to change parameters of this section ONLY
    #SECTION 1. PROCESSING CONTROL PARAMETERS
    #1. suggested default values of the parameters
    ionParam.doIon = False
    ionParam.considerBurstProperties = False
    ionParam.startStep = ionParam.allSteps[0]
    ionParam.endStep = ionParam.allSteps[-1]

    #ionospheric layer height (km)
    ionParam.ionHeight = 200.0
    #before filtering ionosphere, if applying polynomial fitting
    #False: no fitting
    #True: with fitting
    ionParam.ionFit = True
    #window size for filtering ionosphere
    ionParam.ionFilteringWinsizeMax = 200
    ionParam.ionFilteringWinsizeMin = 100
    #window size for filtering azimuth shift caused by ionosphere
    ionParam.ionshiftFilteringWinsizeMax = 150
    ionParam.ionshiftFilteringWinsizeMin = 75
    #correct phase error caused by non-zero center frequency and azimuth shift caused by ionosphere
    #0: no correction
    #1: use mean value of a burst
    #2: use full burst
    ionParam.azshiftFlag = 1

    #better NOT try changing the following two parameters, since they are related
    #to the filtering parameters above
    #number of azimuth looks in the processing of ionosphere estimation
    ionParam.numberAzimuthLooks = 50
    #number of range looks in the processing of ionosphere estimation
    ionParam.numberRangeLooks = 200
    #number of azimuth looks of the interferogram to be unwrapped
    ionParam.numberAzimuthLooks0 = 5*2
    #number of range looks of the interferogram to be unwrapped
    ionParam.numberRangeLooks0 = 20*2


    #2. accept the above parameters from topsApp.py
    ionParam.doIon = self.ION_doIon
    ionParam.considerBurstProperties = self.ION_considerBurstProperties
    ionParam.startStep = self.ION_startStep
    ionParam.endStep = self.ION_endStep

    ionParam.ionHeight = self.ION_ionHeight
    ionParam.ionFit = self.ION_ionFit
    ionParam.ionFilteringWinsizeMax = self.ION_ionFilteringWinsizeMax
    ionParam.ionFilteringWinsizeMin = self.ION_ionFilteringWinsizeMin
    ionParam.ionshiftFilteringWinsizeMax = self.ION_ionshiftFilteringWinsizeMax
    ionParam.ionshiftFilteringWinsizeMin = self.ION_ionshiftFilteringWinsizeMin
    ionParam.azshiftFlag = self.ION_azshiftFlag
    ionParam.maskedAreas = self.ION_maskedAreas

    ionParam.numberAzimuthLooks = self.ION_numberAzimuthLooks
    ionParam.numberRangeLooks = self.ION_numberRangeLooks
    ionParam.numberAzimuthLooks0 = self.ION_numberAzimuthLooks0
    ionParam.numberRangeLooks0 = self.ION_numberRangeLooks0


    #3. check parameters
    #convert to m
    ionParam.ionHeight *= 1000.0

    #check number of looks
    if not ((ionParam.numberAzimuthLooks % ionParam.numberAzimuthLooks0 == 0) and \
       (1 <= ionParam.numberAzimuthLooks0 <= ionParam.numberAzimuthLooks)):
        raise Exception('numberAzimuthLooks must be integer multiples of numberAzimuthLooks0')
    if not ((ionParam.numberRangeLooks % ionParam.numberRangeLooks0 == 0) and \
       (1 <= ionParam.numberRangeLooks0 <= ionParam.numberRangeLooks)):
        raise Exception('numberRangeLooks must be integer multiples of numberRangeLooks0')

    #check steps for ionospheric correction
    if ionParam.startStep not in ionParam.allSteps:
        print('all steps for ionospheric correction in order: {}'.format(ionParam.allSteps))
        raise Exception('please specify the correct start step for ionospheric correction from above list')
    if ionParam.endStep not in ionParam.allSteps:
        print('all steps for ionospheric correction in order: {}'.format(ionParam.allSteps))
        raise Exception('please specify the correct start step for ionospheric correction from above list')
    if ionParam.allSteps.index(ionParam.startStep) > ionParam.allSteps.index(ionParam.endStep):
        print('correct relationship: start step <= end step')
        raise Exception('error: start step is after end step.')
    ###################################################################

    ###################################################################
    #routines that require setting parameters
    #def ionosphere(self, ionParam):
    #def ionSwathBySwath(self, ionParam):
    #def filt_gaussian(self, ionParam):
    #def ionosphere_shift(self, ionParam):
    #def ion2grd(self, ionParam):
    #def esd(self, ionParam):
    ###################################################################

    #SECTION 2. DIRECTORIES AND FILENAMES
    #directories
    ionParam.ionDirname = 'ion'
    ionParam.lowerDirname = 'lower'
    ionParam.upperDirname = 'upper'
    ionParam.ioncalDirname = 'ion_cal'
    ionParam.ionBurstDirname = 'ion_burst'
    #these are same directory names as topsApp.py/TopsProc.py
    #ionParam.referenceSlcProduct = 'reference'
    #ionParam.secondarySlcProduct = 'secondary'
    #ionParam.fineCoregDirname = 'fine_coreg'
    ionParam.fineIfgDirname = 'fine_interferogram'
    ionParam.mergedDirname = 'merged'
    #filenames
    ionParam.ionRawNoProj = 'raw_no_projection.ion'
    ionParam.ionCorNoProj = 'raw_no_projection.cor'
    ionParam.ionRaw = 'raw.ion'
    ionParam.ionCor = 'raw.cor'
    ionParam.ionFilt = 'filt.ion'
    ionParam.ionShift = 'azshift.ion'
    ionParam.warning = 'warning.txt'

    #SECTION 3. DATA PARAMETERS
    #earth's radius (m)
    ionParam.earthRadius = 6371 * 1000.0
    #reference range (m) for moving range center frequency to zero, center of center swath
    ionParam.rgRef = 875714.0
    #range bandwidth (Hz) for splitting, range processingBandwidth: [5.650000000000000e+07, 4.830000000000000e+07, 4.278991840322842e+07]
    ionParam.rgBandwidthForSplit = 40.0 * 10**6
    ionParam.rgBandwidthSub = ionParam.rgBandwidthForSplit / 3.0

    #SECTION 4. DEFINE WAVELENGTHS AND DETERMINE IF CALCULATE IONOSPHERE WITH MERGED INTERFEROGRAM
    getParamFromData = False
    referenceStartingRange = np.zeros(3)
    secondaryStartingRange = np.zeros(3)
    swathList = self._insar.getValidSwathList(self.swaths)
    for swath in swathList:
        ####Load secondary metadata
        reference = self._insar.loadProduct( os.path.join(self._insar.referenceSlcProduct, 'IW{0}.xml'.format(swath)))
        secondary = self._insar.loadProduct( os.path.join(self._insar.secondarySlcProduct, 'IW{0}.xml'.format(swath)))

        ####Indices w.r.t reference
        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
        secondaryBurstStart, secondaryBurstEnd = self._insar.commonSecondaryBurstLimits(swath-1)

        if minBurst == maxBurst:
            #print('Skipping processing of swath {0}'.format(swath))
            continue
        else:
            ii = minBurst
            jj = secondaryBurstStart + ii - minBurst

            masBurst = reference.bursts[ii]
            slvBurst = secondary.bursts[jj]

            #use the 1/3, 1/3, 1/3 scheme for splitting
            ionParam.radarWavelength = masBurst.radarWavelength
            ionParam.radarWavelengthLower = SPEED_OF_LIGHT / (SPEED_OF_LIGHT / ionParam.radarWavelength - ionParam.rgBandwidthForSplit / 3.0)
            ionParam.radarWavelengthUpper = SPEED_OF_LIGHT / (SPEED_OF_LIGHT / ionParam.radarWavelength + ionParam.rgBandwidthForSplit / 3.0)
            #use this to determine which polynomial to use to calculate a ramp when calculating ionosphere for cross A/B interferogram
            ionParam.passDirection = masBurst.passDirection.lower()

            referenceStartingRange[swath-1] = masBurst.startingRange
            secondaryStartingRange[swath-1] = slvBurst.startingRange
            getParamFromData = True

    #determine if calculate ionosphere using merged interferogram
    if np.sum(referenceStartingRange==secondaryStartingRange) != 3:
        ionParam.calIonWithMerged = False
    else:
        ionParam.calIonWithMerged = True
    #for cross Sentinel-1A/B interferogram, always not using merged interferogram
    if reference.mission != secondary.mission:
        ionParam.calIonWithMerged = False
    #there is no need to process swath by swath when there is only one swath
    #ionSwathBySwath only works when number of swaths >=2
    if len(swathList) == 1:
        ionParam.calIonWithMerged = True

    #determine if remove an empirical ramp
    if reference.mission == secondary.mission:
        ionParam.rampRemovel = 0
    else:
        #estimating ionospheric phase for cross Sentinel-1A/B interferogram
        #an empirical ramp will be removed from the estimated ionospheric phase
        if reference.mission == 'S1A' and secondary.mission == 'S1B':
            ionParam.rampRemovel = 1
        else:
            ionParam.rampRemovel = -1

    if getParamFromData == False:
        raise Exception('cannot get parameters from data')

    return ionParam


def next_pow2(a):
    x=2
    while x < a:
        x *= 2
    return x


def removeHammingWindow(inputfile, outputfile, bandwidth, samplingRate, alpha, virtual=True):
    '''
    This function removes the range Hamming window imposed on the signal
    bandwidth:     range bandwidth
    samplingRate:  range sampling rate
    alpha:         alpha of the Hamming window
    '''
    #(length, width) = slc.shape


    inImg = isceobj.createSlcImage()
    inImg.load( inputfile + '.xml')

    width = inImg.getWidth()
    length = inImg.getLength()

    if not virtual:
        slc = np.memmap(inputfile, dtype=np.complex64, mode='r', shape=(length,width))
    else:
        slc = loadVirtualArray(inputfile + '.vrt')

    #fft length
    nfft = next_pow2(width)
    #Hamming window length
    nwin = np.int(np.around(bandwidth / samplingRate*nfft))
    #make it a even number, since we are going to use even fft length
    nwin = ((nwin+1)//2)*2
    #the starting and ending index of window in the spectrum
    start = np.int(np.around((nfft - nwin) / 2))
    end = np.int(np.around(start + nwin - 1))
    hammingWindow = alpha - (1.0-alpha) * np.cos(np.linspace(-np.pi, np.pi, num=nwin, endpoint=True))
    hammingWindow = 1.0/np.fft.fftshift(hammingWindow)
    spec = np.fft.fft(slc, n=nfft, axis=1)
    spec = np.fft.fftshift(spec, axes=1)
    spec[:, start:end+1] *= hammingWindow[None,:]
    spec = np.fft.fftshift(spec, axes=1)
    spec = np.fft.ifft(spec, n=nfft, axis=1)
    slcd  = spec[:, 0:width] * ((slc.real!=0) | (slc.imag!=0))
    #after these fft and ifft, the values are not scaled by constant.

    slcd.astype(np.complex64).tofile(outputfile)
    inImg.setFilename(outputfile)
    inImg.extraFilename = outputfile + '.vrt'
    inImg.setAccessMode('READ')
    inImg.renderHdr()

    return slcd


def runCmd(cmd, silent=0):

    if silent == 0:
        print("{}".format(cmd))
    status = os.system(cmd)
    if status != 0:
        raise Exception('error when running:\n{}\n'.format(cmd))


def adjustValidLineSample(reference,secondary):

    reference_lastValidLine = reference.firstValidLine + reference.numValidLines - 1
    reference_lastValidSample = reference.firstValidSample + reference.numValidSamples - 1
    secondary_lastValidLine = secondary.firstValidLine + secondary.numValidLines - 1
    secondary_lastValidSample = secondary.firstValidSample + secondary.numValidSamples - 1

    igram_lastValidLine = min(reference_lastValidLine, secondary_lastValidLine)
    igram_lastValidSample = min(reference_lastValidSample, secondary_lastValidSample)

    reference.firstValidLine = max(reference.firstValidLine, secondary.firstValidLine)
    reference.firstValidSample = max(reference.firstValidSample, secondary.firstValidSample)

    reference.numValidLines = igram_lastValidLine - reference.firstValidLine + 1
    reference.numValidSamples = igram_lastValidSample - reference.firstValidSample + 1


def multiply2(referencename, secondaryname, fact, rngname=None, ionname=None, infname=None, overlapBox=None, valid=True, virtual=True):
    '''
    This routine forms interferogram and possibly removes topographic and ionospheric phases.
    all the following indexes start from 1
    overlapBox[0]: first line
    overlapBox[1]: last line
    overlapBox[2]: first sample
    overlapBox[3]: last sample
    '''

    #use reference image
    img = isceobj.createSlcImage()
    img.load(referencename + '.xml')
    width = img.getWidth()
    length = img.getLength()

    #reference
    if not virtual:
        reference = np.memmap(referencename, dtype=np.complex64, mode='r', shape=(length,width))
    else:
        reference = loadVirtualArray(referencename + '.vrt')

    #secondary
    secondary = np.memmap(secondaryname, dtype=np.complex64, mode='r', shape=(length, width))

    #interferogram
    cJ = np.complex64(-1j)
    inf = reference[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1] \
        * np.conj(secondary[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1])

    #topography
    if rngname != None:
        rng2 = np.memmap(rngname, dtype=np.float32, mode='r', shape=(length,width))
        inf *= np.exp(cJ*fact*rng2[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1])

    #ionosphere
    if ionname != None:
        ion = np.memmap(ionname, dtype=np.float32, mode='r', shape=(length, width))
        inf *= np.exp(cJ*ion[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1])

    if valid == True:
        inf2 = inf
    else:
        inf2 = np.zeros((length,width), dtype=np.complex64)
        inf2[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1] = inf

    #inf = reference[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1] \
    #    * np.conj(secondary[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1]) \
    #    * np.exp(cJ*ion[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1]) \
    #    * np.exp(cJ*fact*rng2[overlapBox[0]-1:overlapBox[1]-1+1, overlapBox[2]-1:overlapBox[3]-1+1])

    if infname != None:
        inf2.astype(np.complex64).tofile(infname)
        img = isceobj.createIntImage()
        img.setFilename(infname)
        img.extraFilename = infname + '.vrt'
        if valid == True:
            img.setWidth(overlapBox[3]-overlapBox[2]+1)
            img.setLength(overlapBox[1]-overlapBox[0]+1)
        else:
            img.setWidth(width)
            img.setLength(length)
        img.setAccessMode('READ')
        img.renderHdr()

    return inf2


def subband(self, ionParam):
    '''
    generate subband images
    '''
    from isceobj.Sensor.TOPS import createTOPSSwathSLCProduct
    from isceobj.Util.Poly2D import Poly2D
    from contrib.alos2proc.alos2proc import rg_filter

    from isceobj.TopsProc.runFineResamp import resampSecondary
    from isceobj.TopsProc.runFineResamp import getRelativeShifts
    from isceobj.TopsProc.runFineResamp import adjustValidSampleLine
    from isceobj.TopsProc.runFineResamp import getValidLines

    #from isceobj.TopsProc.runBurstIfg import adjustValidLineSample

    print('processing subband burst interferograms')
    virtual = self.useVirtualFiles
    swathList = self._insar.getValidSwathList(self.swaths)
    for swath in swathList:
        ####Load secondary metadata
        reference = self._insar.loadProduct( os.path.join(self._insar.referenceSlcProduct, 'IW{0}.xml'.format(swath)))
        secondary = self._insar.loadProduct( os.path.join(self._insar.secondarySlcProduct, 'IW{0}.xml'.format(swath)))

        dt = secondary.bursts[0].azimuthTimeInterval
        dr = secondary.bursts[0].rangePixelSize

        ###Directory with offsets
        offdir = os.path.join(self._insar.fineOffsetsDirname, 'IW{0}'.format(swath))

        ####Indices w.r.t reference
        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
        secondaryBurstStart, secondaryBurstEnd = self._insar.commonSecondaryBurstLimits(swath-1)

        if minBurst == maxBurst:
            print('Skipping processing of swath {0}'.format(swath))
            continue

        #create dirs
        lowerDir = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname, 'IW{0}'.format(swath))
        upperDir = os.path.join(ionParam.ionDirname, ionParam.upperDirname, ionParam.fineIfgDirname, 'IW{0}'.format(swath))
        os.makedirs(lowerDir, exist_ok=True)
        os.makedirs(upperDir, exist_ok=True)

        ##############################################################
        #for resampling
        relShifts = getRelativeShifts(reference, secondary, minBurst, maxBurst, secondaryBurstStart)
        print('Shifts IW-{0}: '.format(swath), relShifts) 
   
        ####Can corporate known misregistration here
        apoly = Poly2D()
        apoly.initPoly(rangeOrder=0,azimuthOrder=0,coeffs=[[0.]])

        rpoly = Poly2D()
        rpoly.initPoly(rangeOrder=0,azimuthOrder=0,coeffs=[[0.]])

        misreg_az = self._insar.secondaryTimingCorrection / dt
        misreg_rg = self._insar.secondaryRangeCorrection / dr
        ##############################################################

        fineIfgLower = createTOPSSwathSLCProduct()
        fineIfgLower.configure()

        fineIfgUpper = createTOPSSwathSLCProduct()
        fineIfgUpper.configure()

        #only process common bursts
        for ii in range(minBurst, maxBurst):
            jj = secondaryBurstStart + ii - minBurst 
    
            masBurst = reference.bursts[ii]
            slvBurst = secondary.bursts[jj]

            print('processing reference burst: %02d, secondary burst: %02d, swath: %d'%(ii+1, jj+1, swath))
            ################################################################
            #1. removing window and subband
            for ms in ['reference', 'secondary']:
                #setup something
                if ms == 'reference':
                    burst = masBurst
                    #put the temporary file in the lower directory
                    tmpFilename = os.path.join(lowerDir, 'reference_dw_'+os.path.basename(burst.image.filename))
                    tmpFilename2 = 'reference_'+os.path.basename(burst.image.filename)
                else:
                    burst = slvBurst
                    #put the temporary file in the lower directory
                    tmpFilename = os.path.join(lowerDir, 'secondary_dw_'+os.path.basename(burst.image.filename))
                    tmpFilename2 = 'secondary_'+os.path.basename(burst.image.filename)

                #removing window
                rangeSamplingRate = SPEED_OF_LIGHT / (2.0 * burst.rangePixelSize)
                if burst.rangeWindowType == 'Hamming':
                    removeHammingWindow(burst.image.filename, tmpFilename, burst.rangeProcessingBandwidth, rangeSamplingRate, burst.rangeWindowCoefficient, virtual=virtual)  
                else:
                    raise Exception('Range weight window type: {} is not supported yet!'.format(burst.rangeWindowType))
                
                #subband
                rg_filter(tmpFilename,
                          #burst.numberOfSamples,
                          2,
                          [os.path.join(lowerDir, tmpFilename2), os.path.join(upperDir, tmpFilename2)],
                          [ionParam.rgBandwidthSub / rangeSamplingRate, ionParam.rgBandwidthSub / rangeSamplingRate],
                          [-ionParam.rgBandwidthForSplit / 3.0 / rangeSamplingRate, ionParam.rgBandwidthForSplit / 3.0 / rangeSamplingRate],
                          129,
                          512,
                          0.1,
                          0,
                          (burst.startingRange - ionParam.rgRef) / burst.rangePixelSize
                    )

                #remove temporary file
                os.remove(tmpFilename)
                os.remove(tmpFilename+'.xml')
                os.remove(tmpFilename+'.vrt')

            #2. resampling and form interferogram
            #resampling
            try:
                offset = relShifts[jj]
            except:
                raise Exception('Trying to access shift for secondary burst index {0}, which may not overlap with reference for swath {1}'.format(jj, swath))

            ####Setup initial polynomials
            ### If no misregs are given, these are zero
            ### If provided, can be used for resampling without running to geo2rdr again for fast results
            rdict = {'azpoly' : apoly,
                     'rgpoly' : rpoly,
                     'rangeOff' : os.path.join(offdir, 'range_%02d.off'%(ii+1)),
                    'azimuthOff': os.path.join(offdir, 'azimuth_%02d.off'%(ii+1))}


            ###For future - should account for azimuth and range misreg here .. ignoring for now.
            azCarrPoly, dpoly = secondary.estimateAzimuthCarrierPolynomials(slvBurst, offset = -1.0 * offset)

            rdict['carrPoly'] = azCarrPoly
            rdict['doppPoly'] = dpoly


            for lu in ['lower', 'upper']:
                masBurst2 = masBurst.clone()
                slvBurst2 = slvBurst.clone()
                slvBurstResamp2 = masBurst.clone()
                if lu == 'lower':
                    masBurst2.radarWavelength = ionParam.radarWavelengthLower
                    masBurst2.rangeProcessingBandwidth = ionParam.rgBandwidthSub
                    masBurst2.image.filename = os.path.join(lowerDir, 'reference_'+os.path.basename(masBurst.image.filename))
                    slvBurst2.radarWavelength = ionParam.radarWavelengthLower
                    slvBurst2.rangeProcessingBandwidth = ionParam.rgBandwidthSub
                    slvBurst2.image.filename = os.path.join(lowerDir, 'secondary_'+os.path.basename(slvBurst.image.filename))
                    slvBurstResamp2.radarWavelength = ionParam.radarWavelengthLower
                    slvBurstResamp2.rangeProcessingBandwidth = ionParam.rgBandwidthSub
                    slvBurstResamp2.image.filename = os.path.join(lowerDir, 'reference_'+os.path.basename(masBurst.image.filename))
                    outname = os.path.join(lowerDir, 'secondary_resamp_'+os.path.basename(slvBurst.image.filename))
                    ifgdir = lowerDir
                else:
                    masBurst2.radarWavelength = ionParam.radarWavelengthUpper
                    masBurst2.rangeProcessingBandwidth = ionParam.rgBandwidthSub
                    masBurst2.image.filename = os.path.join(upperDir, 'reference_'+os.path.basename(masBurst.image.filename))
                    slvBurst2.radarWavelength = ionParam.radarWavelengthUpper
                    slvBurst2.rangeProcessingBandwidth = ionParam.rgBandwidthSub
                    slvBurst2.image.filename = os.path.join(upperDir, 'secondary_'+os.path.basename(slvBurst.image.filename))
                    slvBurstResamp2.radarWavelength = ionParam.radarWavelengthUpper
                    slvBurstResamp2.rangeProcessingBandwidth = ionParam.rgBandwidthSub
                    slvBurstResamp2.image.filename = os.path.join(upperDir, 'reference_'+os.path.basename(masBurst.image.filename))
                    outname = os.path.join(upperDir, 'secondary_resamp_'+os.path.basename(slvBurst.image.filename))
                    ifgdir = upperDir
                outimg = resampSecondary(masBurst2, slvBurst2, rdict, outname)
                minAz, maxAz, minRg, maxRg = getValidLines(slvBurst2, rdict, outname,
                        misreg_az = misreg_az - offset, misreg_rng = misreg_rg)
                adjustValidSampleLine(slvBurstResamp2, slvBurst2,
                                             minAz=minAz, maxAz=maxAz,
                                             minRng=minRg, maxRng=maxRg)
                slvBurstResamp2.image.filename = outimg.filename
                
                #forming interferogram
                referencename = masBurst2.image.filename
                secondaryname = slvBurstResamp2.image.filename
                rngname = os.path.join(offdir, 'range_%02d.off'%(ii+1))
                infname = os.path.join(ifgdir, 'burst_%02d.int'%(ii+1))

                fact = 4.0 * np.pi * slvBurstResamp2.rangePixelSize / slvBurstResamp2.radarWavelength
                adjustValidLineSample(masBurst2,slvBurstResamp2)


                #in original runBurstIfg.py, valid samples in the interferogram are the following (indexes in the numpy matrix):
                #referenceFrame.firstValidLine:referenceFrame.firstValidLine + referenceFrame.numValidLines, referenceFrame.firstValidSample:referenceFrame.firstValidSample + referenceFrame.numValidSamples
                #after the following processing, valid samples in the interferogram are the following (indexes in the numpy matrix):
                #[masBurst.firstValidLine:masBurst.firstValidLine + masBurst.numValidLines, masBurst.firstValidSample:masBurst.firstValidSample + masBurst.numValidSamples]
                #SO THEY ARE EXACTLY THE SAME
                firstline   = masBurst2.firstValidLine + 1
                lastline    = firstline + masBurst2.numValidLines - 1
                firstcolumn = masBurst2.firstValidSample + 1
                lastcolumn  = firstcolumn + masBurst2.numValidSamples - 1
                overlapBox = [firstline, lastline, firstcolumn, lastcolumn]
                multiply2(referencename, secondaryname, fact, rngname=rngname, ionname=None, infname=infname, overlapBox=overlapBox, valid=False, virtual=virtual)

                #directly from multiply() of runBurstIfg.py
                img = isceobj.createIntImage()
                img.setFilename(infname)
                img.setWidth(masBurst2.numberOfSamples)
                img.setLength(masBurst2.numberOfLines)
                img.setAccessMode('READ')
                #img.renderHdr()

                #save it for deleting later
                masBurst2_filename = masBurst2.image.filename
                #change it for interferogram
                masBurst2.image = img

                if lu == 'lower':
                    fineIfgLower.bursts.append(masBurst2)
                else:
                    fineIfgUpper.bursts.append(masBurst2)

                #remove reference and secondary subband slcs
                os.remove(masBurst2_filename)
                os.remove(masBurst2_filename+'.xml')
                os.remove(masBurst2_filename+'.vrt')
                os.remove(slvBurst2.image.filename)
                os.remove(slvBurst2.image.filename+'.xml')
                os.remove(slvBurst2.image.filename+'.vrt')
                os.remove(slvBurstResamp2.image.filename)
                os.remove(slvBurstResamp2.image.filename+'.xml')
                os.remove(slvBurstResamp2.image.filename+'.vrt')

        fineIfgLower.numberOfBursts = len(fineIfgLower.bursts)
        fineIfgUpper.numberOfBursts = len(fineIfgUpper.bursts)
        self._insar.saveProduct(fineIfgLower, os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname, 'IW{0}.xml'.format(swath)))
        self._insar.saveProduct(fineIfgUpper, os.path.join(ionParam.ionDirname, ionParam.upperDirname, ionParam.fineIfgDirname, 'IW{0}.xml'.format(swath)))


def cal_coherence(inf, win=5, edge=0):
    '''
    compute coherence uisng only interferogram (phase).
    This routine still follows the regular equation for computing coherence,
    but assumes the amplitudes of reference and secondary are one, so that coherence
    can be computed using phase only.

    inf: interferogram
    win: window size
    edge: 0: remove all non-full convolution samples

          1: remove samples computed from less than half convolution
             (win=5 used to illustration below)
             * * *
             * * *
             * * *
             * * *
             * * *

          2: remove samples computed from less than quater convolution
             (win=5 used to illustration below)
             * * *
             * * *
             * * *

          3: remove non-full convolution samples on image edges

          4: keep all samples
    '''
    import scipy.signal as ss

    if win % 2 != 1:
        raise Exception('window size must be odd!')
    hwin = np.int(np.around((win - 1) / 2))

    filt = np.ones((win, win))
    amp  = np.absolute(inf)

    cnt = ss.convolve2d((amp!=0), filt, mode='same')
    cor = ss.convolve2d(inf/(amp + (amp==0)), filt, mode='same')
    cor = (amp!=0) * np.absolute(cor) / (cnt + (cnt==0))

    #trim edges
    if edge == 0:
        num = win * win
        cor[np.nonzero(cnt < num)] = 0.0
    elif edge == 1:
        num = win * (hwin+1)
        cor[np.nonzero(cnt < num)] = 0.0
    elif edge == 2:
        num = (hwin+1) * (hwin+1)
        cor[np.nonzero(cnt < num)] = 0.0
    elif edge == 3:
        cor[0:hwin, :] = 0.0
        cor[-hwin:, :] = 0.0
        cor[:, 0:hwin] = 0.0
        cor[:, -hwin:] = 0.0
    else:
        pass

    #print("coherence, max: {} min: {}".format(np.max(cor[np.nonzero(cor!=0)]), np.min(cor[np.nonzero(cor!=0)])))
    return cor


def getMergeBox(self, xmlDirname, numberRangeLooks=1, numberAzimuthLooks=1):
    '''
    xmlDirname:         directory containing xml file
    numberRangeLooks:   number of range looks to take after merging
    numberAzimuthLooks: number of azimuth looks to take after merging
    '''

    from isceobj.TopsProc.runMergeBursts import mergeBox
    from isceobj.TopsProc.runMergeBursts import adjustValidWithLooks

    swathList = self._insar.getValidSwathList(self.swaths)

    #get bursts
    frames=[]
    for swath in swathList:
        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
        if minBurst==maxBurst:
            #print('Skipping processing of swath {0}'.format(swath))
            continue
        #since burst directory does not necessarily has IW*.xml, we use the following dir
        #ifg = self._insar.loadProduct( os.path.join(self._insar.fineIfgDirname, 'IW{0}.xml'.format(swath)))
        #use lower
        #dirname = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname)
        ifg = self._insar.loadProduct( os.path.join(xmlDirname, 'IW{0}.xml'.format(swath)))
        frames.append(ifg)

    #determine merged size
    box = mergeBox(frames)
    #adjust valid with looks, 'frames' ARE CHANGED AFTER RUNNING THIS
    (burstValidBox, burstValidBox2, message) = adjustValidWithLooks(frames, box, numberAzimuthLooks, numberRangeLooks, edge=0, avalid='strict', rvalid='strict')

    return (box, burstValidBox, burstValidBox2, frames)


def merge(self, ionParam):
    '''
    merge burst interferograms and compute coherence
    '''
    from isceobj.TopsProc.runMergeBursts import mergeBox
    from isceobj.TopsProc.runMergeBursts import adjustValidWithLooks
    from isceobj.TopsProc.runMergeBursts import mergeBurstsVirtual
    from isceobj.TopsProc.runMergeBursts import multilook as multilook2

    #merge burst interferograms
    mergeFilename = self._insar.mergedIfgname
    xmlDirname = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname)
    dirs = [ionParam.lowerDirname, ionParam.upperDirname]
    for dirx in dirs:
        mergeDirname = os.path.join(ionParam.ionDirname, dirx, ionParam.mergedDirname)
        burstDirname = os.path.join(ionParam.ionDirname, dirx, ionParam.fineIfgDirname)

        frames=[]
        burstList = []
        swathList = self._insar.getValidSwathList(self.swaths)
        for swath in swathList:
            minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
            if minBurst==maxBurst:
                continue
            ifg = self._insar.loadProduct( os.path.join(xmlDirname, 'IW{0}.xml'.format(swath)))
            frames.append(ifg)
            burstList.append([os.path.join(burstDirname, 'IW{0}'.format(swath),  'burst_%02d.int'%(x+1)) for x in range(minBurst, maxBurst)])

        os.makedirs(mergeDirname, exist_ok=True)

        suffix = '.full'
        if (ionParam.numberRangeLooks0 == 1) and (ionParam.numberAzimuthLooks0 == 1):
            suffix=''

        box = mergeBox(frames)
        #adjust valid with looks, 'frames' ARE CHANGED AFTER RUNNING THIS
        #here numberRangeLooks, instead of numberRangeLooks0, is used, since we need to do next step multilooking after unwrapping. same for numberAzimuthLooks.
        (burstValidBox, burstValidBox2, message) = adjustValidWithLooks(frames, box, ionParam.numberAzimuthLooks, ionParam.numberRangeLooks, edge=0, avalid='strict', rvalid='strict')
        mergeBurstsVirtual(frames, burstList, box, os.path.join(mergeDirname, mergeFilename+suffix))
        if suffix not in ['',None]:
            multilook2(os.path.join(mergeDirname, mergeFilename+suffix),
              outname = os.path.join(mergeDirname, mergeFilename),
              alks = ionParam.numberAzimuthLooks0, rlks=ionParam.numberRangeLooks0)
        #this is never used for ionosphere correction
        else:
            print('Skipping multi-looking ....')

    #The orginal coherence calculated by topsApp.py is not good at all, use the following coherence instead
    lowerintfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname, self._insar.mergedIfgname)
    upperintfile = os.path.join(ionParam.ionDirname, ionParam.upperDirname, ionParam.mergedDirname, self._insar.mergedIfgname)
    corfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname, self._insar.correlationFilename)

    img = isceobj.createImage()
    img.load(lowerintfile + '.xml')
    width = img.width
    length = img.length
    lowerint = np.fromfile(lowerintfile, dtype=np.complex64).reshape(length, width)
    upperint = np.fromfile(upperintfile, dtype=np.complex64).reshape(length, width)

    #compute coherence only using interferogram
    #here I use differential interferogram of lower and upper band interferograms
    #so that coherence is not affected by fringes
    cord = cal_coherence(lowerint*np.conjugate(upperint), win=3, edge=4)
    cor = np.zeros((length*2, width), dtype=np.float32)
    cor[0:length*2:2, :] = np.sqrt( (np.absolute(lowerint)+np.absolute(upperint))/2.0 )
    cor[1:length*2:2, :] = cord
    cor.astype(np.float32).tofile(corfile)

    #create xml and vrt
    #img.scheme = 'BIL'
    #img.bands = 2
    #img.filename = corfile
    #img.renderHdr()

    #img = isceobj.Image.createUnwImage()
    img = isceobj.createOffsetImage()
    img.setFilename(corfile)
    img.extraFilename = corfile + '.vrt'
    img.setWidth(width)
    img.setLength(length)
    img.renderHdr()


def renameFile(oldname, newname):
    img = isceobj.createImage()
    img.load(oldname + '.xml')
    img.setFilename(newname)
    img.extraFilename = newname+'.vrt'
    img.renderHdr()

    os.rename(oldname, newname)
    os.remove(oldname + '.xml')
    os.remove(oldname + '.vrt')


def maskUnwrap(unwfile, maskfile):
    tmpfile = 'tmp.unw'
    renameFile(unwfile, tmpfile)
    cmd = "imageMath.py -e='a_0*(abs(b)!=0);a_1*(abs(b)!=0)' --a={0} --b={1} -s BIL -o={2}".format(tmpfile, maskfile, unwfile)
    runCmd(cmd)
    os.remove(tmpfile)
    os.remove(tmpfile+'.xml')
    os.remove(tmpfile+'.vrt')


def snaphuUnwrap(self, xmlDirname, wrapName, corrfile, unwrapName, nrlks, nalks, costMode = 'DEFO',initMethod = 'MST', defomax = 4.0, initOnly = False):
        #runUnwrap(self,                                           costMode = 'SMOOTH',initMethod = 'MCF', defomax = 2,   initOnly = True)
    '''
    xmlDirname:  xml dir name
    wrapName:    input interferogram
    corrfile:    input coherence file
    unwrapName:  output unwrapped interferogram
    nrlks:       number of range looks of the interferogram
    nalks:       number of azimuth looks of the interferogram
    '''

    from contrib.Snaphu.Snaphu import Snaphu
    from isceobj.Planet.Planet import Planet

    img = isceobj.createImage()
    img.load(wrapName + '.xml')
    width = img.getWidth()

    #get altitude
    swathList = self._insar.getValidSwathList(self.swaths)
    for swath in swathList[0:1]:
        ifg = self._insar.loadProduct( os.path.join(xmlDirname, 'IW{0}.xml'.format(swath)))
        wavelength = ifg.bursts[0].radarWavelength

        ####tmid 
        tstart = ifg.bursts[0].sensingStart
        tend   = ifg.bursts[-1].sensingStop
        tmid = tstart + 0.5*(tend - tstart)

        #14-APR-2018
        burst_index = np.int(np.around(len(ifg.bursts)/2))
        orbit = ifg.bursts[burst_index].orbit
        peg = orbit.interpolateOrbit(tmid, method='hermite')

        refElp = Planet(pname='Earth').ellipsoid
        llh = refElp.xyz_to_llh(peg.getPosition())
        hdg = orbit.getENUHeading(tmid)
        refElp.setSCH(llh[0], llh[1], hdg)
        earthRadius = refElp.pegRadCur
        altitude   = llh[2]

    rangeLooks = nrlks
    azimuthLooks = nalks
    azfact = 0.8
    rngfact = 0.8
    corrLooks = rangeLooks * azimuthLooks/(azfact*rngfact) 
    maxComponents = 20

    snp = Snaphu()
    snp.setInitOnly(initOnly)
    snp.setInput(wrapName)
    snp.setOutput(unwrapName)
    snp.setWidth(width)
    snp.setCostMode(costMode)
    snp.setEarthRadius(earthRadius)
    snp.setWavelength(wavelength)
    snp.setAltitude(altitude)
    snp.setCorrfile(corrfile)
    snp.setInitMethod(initMethod)
    snp.setCorrLooks(corrLooks)
    snp.setMaxComponents(maxComponents)
    snp.setDefoMaxCycles(defomax)
    snp.setRangeLooks(rangeLooks)
    snp.setAzimuthLooks(azimuthLooks)
    #snp.setCorFileFormat('FLOAT_DATA')
    snp.prepare()
    snp.unwrap()

    ######Render XML
    outImage = isceobj.Image.createUnwImage()
    outImage.setFilename(unwrapName)
    outImage.setWidth(width)
    outImage.setAccessMode('read')
    outImage.renderVRT()
    outImage.createImage()
    outImage.finalizeImage()
    outImage.renderHdr()

    #####Check if connected components was created
    if snp.dumpConnectedComponents:
        connImage = isceobj.Image.createImage()
        connImage.setFilename(unwrapName+'.conncomp')
        connImage.setWidth(width)
        connImage.setAccessMode('read')
        connImage.setDataType('BYTE')
        connImage.renderVRT()
        connImage.createImage()
        connImage.finalizeImage()
        connImage.renderHdr()

    return


def unwrap(self, ionParam):
    '''
    unwrap lower and upper band interferograms
    '''

    print('unwrapping lower and upper band interferograms')
    dirs = [ionParam.lowerDirname, ionParam.upperDirname]
    #there is only one coherence file in lower directory
    corfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname, self._insar.correlationFilename)
    for dirx in dirs:
        procdir = os.path.join(ionParam.ionDirname, dirx, ionParam.mergedDirname)
        wrapName = os.path.join(procdir, self._insar.mergedIfgname)
        unwrapName = os.path.join(procdir, self._insar.unwrappedIntFilename)
        xmlDirname = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname)
        #unwrap
        snaphuUnwrap(self, xmlDirname, wrapName, corfile, unwrapName, ionParam.numberRangeLooks0, ionParam.numberAzimuthLooks0, costMode = 'SMOOTH',initMethod = 'MCF', defomax = 2, initOnly = True)
        #remove wired things in no-data area
        maskUnwrap(unwrapName, wrapName)

    if [ionParam.numberRangeLooks0, ionParam.numberAzimuthLooks0] != [ionParam.numberRangeLooks, ionParam.numberAzimuthLooks]:
        multilook_unw(self, ionParam, ionParam.mergedDirname)


def multilook_unw(self, ionParam, mergedDirname):
    '''
    30-APR-2018
    This routine moves the original unwrapped files to a directory and takes looks
    '''
    from isceobj.TopsProc.runMergeBursts import multilook as multilook2

    oridir0 = '{}rlks_{}alks'.format(ionParam.numberRangeLooks0, ionParam.numberAzimuthLooks0)
    dirs = [ionParam.lowerDirname, ionParam.upperDirname]
    corName = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname, oridir0, self._insar.correlationFilename)
    for dirx in dirs:
        procdir = os.path.join(ionParam.ionDirname, dirx, mergedDirname)
        #create a directory for original files
        oridir = os.path.join(procdir, oridir0)
        os.makedirs(oridir, exist_ok=True)
        #move files, renameFile uses os.rename, which overwrites if file already exists in oridir. This can support re-run
        filename0 = os.path.join(procdir, self._insar.mergedIfgname)
        filename  = os.path.join(oridir, self._insar.mergedIfgname)
        if os.path.isfile(filename0):
            renameFile(filename0, filename)
        filename0 = os.path.join(procdir, self._insar.unwrappedIntFilename)
        filename  = os.path.join(oridir, self._insar.unwrappedIntFilename)
        if os.path.isfile(filename0):
            renameFile(filename0, filename)
        filename0 = os.path.join(procdir, self._insar.unwrappedIntFilename+'.conncomp')
        filename  = os.path.join(oridir, self._insar.unwrappedIntFilename+'.conncomp')
        if os.path.isfile(filename0):
            renameFile(filename0, filename)
        filename0 = os.path.join(procdir, self._insar.correlationFilename)
        filename  = os.path.join(oridir, self._insar.correlationFilename)
        if os.path.isfile(filename0):
            renameFile(filename0, filename)
        #for topophase.flat.full, move directly
        filename0 = os.path.join(procdir, self._insar.mergedIfgname+'.full.vrt')
        filename  = os.path.join(oridir, self._insar.mergedIfgname+'.full.vrt')
        if os.path.isfile(filename0):
            os.rename(filename0, filename)
        filename0 = os.path.join(procdir, self._insar.mergedIfgname+'.full.xml')
        filename  = os.path.join(oridir, self._insar.mergedIfgname+'.full.xml')
        if os.path.isfile(filename0):
            os.rename(filename0, filename)

        #multi-looking
        nrlks = np.int(np.around(ionParam.numberRangeLooks / ionParam.numberRangeLooks0))
        nalks = np.int(np.around(ionParam.numberAzimuthLooks / ionParam.numberAzimuthLooks0))
        #coherence
        if dirx == ionParam.lowerDirname:
            corName0 = os.path.join(oridir, self._insar.correlationFilename)
            corimg = isceobj.createImage()
            corimg.load(corName0 + '.xml')
            width = corimg.width
            length = corimg.length
            widthNew = np.int(width / nrlks)
            lengthNew = np.int(length / nalks)
            cor0 = (np.fromfile(corName0, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
            amp0 = (np.fromfile(corName0, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
            wgt = cor0**2
            a = multilook(wgt, nalks, nrlks)
            b = multilook(cor0, nalks, nrlks)
            c = multilook(amp0**2, nalks, nrlks)
            d = multilook((cor0!=0).astype(np.int), nalks, nrlks)
            #coherence after multiple looking
            cor = np.zeros((lengthNew*2, widthNew), dtype=np.float32)
            cor[0:lengthNew*2:2, :] = np.sqrt(c / (d + (d==0)))
            cor[1:lengthNew*2:2, :] = b / (d + (d==0))
            #output file
            corName = os.path.join(procdir, self._insar.correlationFilename)
            cor.astype(np.float32).tofile(corName)
            corimg.setFilename(corName)
            corimg.extraFilename = corName + '.vrt'
            corimg.setWidth(widthNew)
            corimg.setLength(lengthNew)
            corimg.renderHdr()
        #unwrapped file
        unwrapName0 = os.path.join(oridir, self._insar.unwrappedIntFilename)
        unwimg = isceobj.createImage()
        unwimg.load(unwrapName0 + '.xml')
        unw0 = (np.fromfile(unwrapName0, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
        amp0 = (np.fromfile(unwrapName0, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
        e = multilook(unw0*wgt, nalks, nrlks)
        f = multilook(amp0**2, nalks, nrlks)
        unw = np.zeros((lengthNew*2, widthNew), dtype=np.float32)
        unw[0:lengthNew*2:2, :] = np.sqrt(f / (d + (d==0)))
        unw[1:lengthNew*2:2, :] = e / (a + (a==0))
        #output file
        unwrapName = os.path.join(procdir, self._insar.unwrappedIntFilename)
        unw.astype(np.float32).tofile(unwrapName)
        unwimg.setFilename(unwrapName)
        unwimg.extraFilename = unwrapName + '.vrt'
        unwimg.setWidth(widthNew)
        unwimg.setLength(lengthNew)
        unwimg.renderHdr()

    #looks like the above is not a good coherence, re-calculate here
    #here I use differential interferogram of lower and upper band interferograms
    #so that coherence is not affected by fringes
    lowerIntName0 = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, mergedDirname, oridir0, self._insar.mergedIfgname)
    upperIntName0 = os.path.join(ionParam.ionDirname, ionParam.upperDirname, mergedDirname, oridir0, self._insar.mergedIfgname)
    lowerIntName = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, mergedDirname, self._insar.mergedIfgname)
    upperIntName = os.path.join(ionParam.ionDirname, ionParam.upperDirname, mergedDirname, self._insar.mergedIfgname)
    #cmd = 'looks.py -i {} -o {} -r {} -a {}'.format(lowerIntName0, lowerIntName, nrlks, nalks)
    #runCmd(cmd)
    #cmd = 'looks.py -i {} -o {} -r {} -a {}'.format(upperIntName0, upperIntName, nrlks, nalks)
    #runCmd(cmd)
    multilook2(lowerIntName0, outname = lowerIntName, alks = nalks, rlks=nrlks)
    multilook2(upperIntName0, outname = upperIntName, alks = nalks, rlks=nrlks)
    lowerint = np.fromfile(lowerIntName, dtype=np.complex64).reshape(lengthNew, widthNew)
    upperint = np.fromfile(upperIntName, dtype=np.complex64).reshape(lengthNew, widthNew)
    cor = np.zeros((lengthNew*2, widthNew), dtype=np.float32)
    cor[0:length*2:2, :] = np.sqrt( (np.absolute(lowerint)+np.absolute(upperint))/2.0 )
    cor[1:length*2:2, :] = cal_coherence(lowerint*np.conjugate(upperint), win=3, edge=4)
    cor.astype(np.float32).tofile(corName)


def create_multi_index2(width2, l1, l2):
    #for number of looks of l1 and l2
    #calculate the correponding index number of l2 in the l1 array
    #applies to both range and azimuth direction

    return ((l2 - l1) / 2.0  + np.arange(width2) * l2) / l1


def fit_surface(x, y, z, wgt, order):
    # x: x coordinate, a column vector
    # y: y coordinate, a column vector
    # z: z coordinate, a column vector
    # wgt: weight of the data points, a column vector


    #number of data points
    m = x.shape[0]
    l = np.ones((m,1), dtype=np.float64)

#    #create polynomial
#    if order == 1:
#        #order of estimated coefficents: 1, x, y
#        a1 = np.concatenate((l, x, y), axis=1)
#    elif order == 2:
#        #order of estimated coefficents: 1, x, y, x*y, x**2, y**2
#        a1 = np.concatenate((l, x, y, x*y, x**2, y**2), axis=1)
#    elif order == 3:
#        #order of estimated coefficents: 1, x, y, x*y, x**2, y**2, x**2*y, y**2*x, x**3, y**3
#        a1 = np.concatenate((l, x, y, x*y, x**2, y**2, x**2*y, y**2*x, x**3, y**3), axis=1)
#    else:
#        raise Exception('order not supported yet\n')

    if order < 1:
        raise Exception('order must be larger than 1.\n')

    #create polynomial
    a1 = l;
    for i in range(1, order+1):
        for j in range(i+1):
            a1 = np.concatenate((a1, x**(i-j)*y**(j)), axis=1)

    #number of variable to be estimated
    n = a1.shape[1]

    #do the least squares
    a = a1 * np.matlib.repmat(np.sqrt(wgt), 1, n)
    b = z * np.sqrt(wgt)
    c = np.linalg.lstsq(a, b, rcond=-1)[0]
    
    #type: <class 'numpy.ndarray'>
    return c


def cal_surface(x, y, c, order):
    #x: x coordinate, a row vector
    #y: y coordinate, a column vector
    #c: coefficients of polynomial from fit_surface
    #order: order of polynomial

    if order < 1:
        raise Exception('order must be larger than 1.\n')

    #number of lines
    length = y.shape[0]
    #number of columns, if row vector, only one element in the shape tuple
    #width = x.shape[1]
    width = x.shape[0]

    x = np.matlib.repmat(x, length, 1)
    y = np.matlib.repmat(y, 1, width)
    z = c[0] * np.ones((length,width), dtype=np.float64)

    index = 0
    for i in range(1, order+1):
        for j in range(i+1):
            index += 1
            z += c[index] * x**(i-j)*y**(j)

    return z


def weight_fitting(ionos, cor, width, length, nrli, nali, nrlo, nalo, order, coth):
    '''
    ionos:  input ionospheric phase
    cor:    coherence of the interferogram
    width:  file width
    length: file length
    nrli:   number of range looks of the input interferograms
    nali:   number of azimuth looks of the input interferograms
    nrlo:   number of range looks of the output ionosphere phase
    nalo:   number of azimuth looks of the ioutput ionosphere phase
    order:  the order of the polynomial for fitting ionosphere phase estimates
    coth:   coherence threshhold for ionosphere phase estimation
    '''

    lengthi = int(length/nali)
    widthi = int(width/nrli)
    lengtho = int(length/nalo)
    widtho = int(width/nrlo)

    #calculate output index
    rgindex = create_multi_index2(widtho, nrli, nrlo)
    azindex = create_multi_index2(lengtho, nali, nalo)

    #convert coherence to weight
    cor = cor**2/(1.009-cor**2)

    #look for data to use
    flag = (cor>coth)*(ionos!=0)
    point_index = np.nonzero(flag)
    m = point_index[0].shape[0]

    #calculate input index matrix
    x0=np.matlib.repmat(np.arange(widthi), lengthi, 1)
    y0=np.matlib.repmat(np.arange(lengthi).reshape(lengthi, 1), 1, widthi)

    x = x0[point_index].reshape(m, 1)
    y = y0[point_index].reshape(m, 1)
    z = ionos[point_index].reshape(m, 1)
    w = cor[point_index].reshape(m, 1)

    #convert to higher precision type before use
    x=np.asfarray(x,np.float64)
    y=np.asfarray(y,np.float64)
    z=np.asfarray(z,np.float64)
    w=np.asfarray(w,np.float64)
    coeff = fit_surface(x, y, z, w, order)

    #convert to higher precision type before use
    rgindex=np.asfarray(rgindex,np.float64)
    azindex=np.asfarray(azindex,np.float64)
    phase_fit = cal_surface(rgindex, azindex.reshape(lengtho, 1), coeff, order)

    #format: widtho, lengtho, single band float32
    return phase_fit


def computeIonosphere(lowerUnw, upperUnw, cor, fl, fu, adjFlag, corThresholdAdj, dispersive):
    '''
    This routine computes ionosphere and remove the relative phase unwrapping errors

    lowerUnw:        lower band unwrapped interferogram
    upperUnw:        upper band unwrapped interferogram
    cor:             coherence
    fl:              lower band center frequency
    fu:              upper band center frequency
    adjFlag:         method for removing relative phase unwrapping errors
                       0: mean value
                       1: polynomial
    corThresholdAdj: coherence threshold of samples used in removing relative phase unwrapping errors
    dispersive:      compute dispersive or non-dispersive
                       0: dispersive
                       1: non-dispersive
    '''

    #use image size from lower unwrapped interferogram
    (length, width)=lowerUnw.shape

##########################################################################################
    # ADJUST PHASE USING MEAN VALUE
    # #ajust phase of upper band to remove relative phase unwrapping errors
    # flag = (lowerUnw!=0)*(cor>=ionParam.corThresholdAdj)
    # index = np.nonzero(flag!=0)
    # mv = np.mean((lowerUnw - upperUnw)[index], dtype=np.float64)
    # print('mean value of phase difference: {}'.format(mv))
    # flag2 = (lowerUnw!=0)
    # index2 = np.nonzero(flag2)
    # #phase for adjustment
    # unwd = ((lowerUnw - upperUnw)[index2] - mv) / (2.0*np.pi)
    # unw_adj = np.around(unwd) * (2.0*np.pi)
    # #ajust phase of upper band
    # upperUnw[index2] += unw_adj
    # unw_diff = lowerUnw - upperUnw
    # print('after adjustment:')
    # print('max phase difference: {}'.format(np.amax(unw_diff)))
    # print('min phase difference: {}'.format(np.amin(unw_diff)))
##########################################################################################
    #adjust phase using mean value
    if adjFlag == 0:
        flag = (lowerUnw!=0)*(cor>=corThresholdAdj)
        index = np.nonzero(flag!=0)
        mv = np.mean((lowerUnw - upperUnw)[index], dtype=np.float64)
        print('mean value of phase difference: {}'.format(mv))
        diff = mv
    #adjust phase using a surface
    else:
        diff = weight_fitting(lowerUnw - upperUnw, cor, width, length, 1, 1, 1, 1, 2, corThresholdAdj)

    flag2 = (lowerUnw!=0)
    index2 = np.nonzero(flag2)
    #phase for adjustment
    unwd = ((lowerUnw - upperUnw) - diff)[index2] / (2.0*np.pi)
    unw_adj = np.around(unwd) * (2.0*np.pi)
    #ajust phase of upper band
    upperUnw[index2] += unw_adj

    unw_diff = (lowerUnw - upperUnw)[index2]
    print('after adjustment:')
    print('max phase difference: {}'.format(np.amax(unw_diff)))
    print('min phase difference: {}'.format(np.amin(unw_diff)))
    print('max-min: {}'.format(np.amax(unw_diff) - np.amin(unw_diff)    ))

    #ionosphere
    #fl = SPEED_OF_LIGHT / ionParam.radarWavelengthLower
    #fu = SPEED_OF_LIGHT / ionParam.radarWavelengthUpper
    f0 = (fl + fu) / 2.0
    
    #dispersive
    if dispersive == 0:
        ionos = fl * fu * (lowerUnw * fu - upperUnw * fl) / f0 / (fu**2 - fl**2)
    #non-dispersive phase
    else:
        ionos = f0 * (upperUnw*fu - lowerUnw * fl) / (fu**2 - fl**2)

    return ionos


def ionosphere(self, ionParam):

    ###################################
    #SET PARAMETERS HERE
    #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self)
    corThresholdAdj = 0.85
    ###################################

    print('computing ionosphere')
    #get files
    lowerUnwfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname, self._insar.unwrappedIntFilename)
    upperUnwfile = os.path.join(ionParam.ionDirname, ionParam.upperDirname, ionParam.mergedDirname, self._insar.unwrappedIntFilename)
    corfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname, self._insar.correlationFilename)

    #use image size from lower unwrapped interferogram
    img = isceobj.createImage()
    img.load(lowerUnwfile + '.xml')
    width = img.width
    length = img.length

    lowerUnw = (np.fromfile(lowerUnwfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    upperUnw = (np.fromfile(upperUnwfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    lowerAmp = (np.fromfile(lowerUnwfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
    upperAmp = (np.fromfile(upperUnwfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
    cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    amp = np.sqrt(lowerAmp**2+upperAmp**2)

    #masked out user-specified areas
    if ionParam.maskedAreas != None:
        maskedAreas = reformatMaskedAreas(ionParam.maskedAreas, length, width)
        for area in maskedAreas:
            lowerUnw[area[0]:area[1], area[2]:area[3]] = 0
            upperUnw[area[0]:area[1], area[2]:area[3]] = 0
            cor[area[0]:area[1], area[2]:area[3]] = 0

    #compute ionosphere
    fl = SPEED_OF_LIGHT / ionParam.radarWavelengthLower
    fu = SPEED_OF_LIGHT / ionParam.radarWavelengthUpper
    adjFlag = 1
    ionos = computeIonosphere(lowerUnw, upperUnw, cor, fl, fu, adjFlag, corThresholdAdj, 0)

    #dump ionosphere
    outDir = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname)
    os.makedirs(outDir, exist_ok=True)
    outFilename = os.path.join(outDir, ionParam.ionRawNoProj)
    ion = np.zeros((length*2, width), dtype=np.float32)
    ion[0:length*2:2, :] = amp
    ion[1:length*2:2, :] = ionos
    ion.astype(np.float32).tofile(outFilename)
    img.filename = outFilename
    img.extraFilename = outFilename + '.vrt'
    img.renderHdr()

    #dump coherence
    outFilename = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCorNoProj)
    ion[1:length*2:2, :] = cor
    ion.astype(np.float32).tofile(outFilename)
    img.filename = outFilename
    img.extraFilename = outFilename + '.vrt'
    img.renderHdr()


def cal_cross_ab_ramp(swathList, width, numberRangeLooks, passDirection):
    '''
    calculate an empirical ramp between Sentinel-1A/B
    29-JUN-2018

    swathList:        self._insar.getValidSwathList(self.swaths)
    width:            single-look image width after merging
    numberRangeLooks: number of range looks in the processing of ionosphere estimation
    passDirection:    descending/ascending
    '''
    
    #below is from processing chile_d156_160725(S1A)-160929(S1B)
    #empirical polynomial
    deg = 3
    if passDirection.lower() == 'descending':
        p = np.array([0.95381267, 2.95567604, -4.56047084, 1.05443172])
    elif passDirection.lower() == 'ascending':
        #for ascending, the polynomial is left/right flipped
        p = np.array([-0.95381267, 5.81711404, -4.21231923, 0.40344958])
    else:
        raise Exception('unknown passDirection! should be either descending or ascending')

    #ca/a166/process/160807-170305 also has the swath offset almost equal to these
    #swath offset in single-look range pixels
    swath_offset = [0, 19810, 43519]
    #total number of single-look range pixels
    tnp = 69189

    #getting x
    nswath = len(swathList)
    if nswath == 3:
        width2 = np.int(width/numberRangeLooks)
        x = np.arange(width2) / (width2 - 1.0)
    else:
        width2 = np.int(width/numberRangeLooks)
        #WARNING: what if the some swaths does not have bursts, and are not merged?
        #         here I just simply ignore this case
        offset = swath_offset[swathList[0]-1]
        x = offset / tnp + width / tnp * np.arange(width2) / (width2 - 1.0)
        
    #calculate ramp
    y_fit = x * 0.0
    for i in range(deg+1):
        y_fit += p[i] * x**[deg-i]

    return y_fit


def ionSwathBySwath(self, ionParam):
    '''
    This routine merge, unwrap and compute ionosphere swath by swath, and then
    adjust phase difference between adjacent swaths caused by relative range timing
    error between adjacent swaths.
    
    This routine includes the following steps in the merged-swath processing:
    merge(self, ionParam)
    unwrap(self, ionParam)
    ionosphere(self, ionParam)
    '''

    from isceobj.TopsProc.runMergeBursts import mergeBox
    from isceobj.TopsProc.runMergeBursts import adjustValidWithLooks
    from isceobj.TopsProc.runMergeBursts import mergeBurstsVirtual
    from isceobj.TopsProc.runMergeBursts import multilook as multilook2

    #########################################
    #SET PARAMETERS HERE
    numberRangeLooks = ionParam.numberRangeLooks
    numberAzimuthLooks = ionParam.numberAzimuthLooks
    numberRangeLooks0 = ionParam.numberRangeLooks0
    numberAzimuthLooks0 = ionParam.numberAzimuthLooks0

    #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self)
    corThresholdSwathAdj = 0.85
    corThresholdAdj = 0.85
    #########################################

    print('computing ionosphere swath by swath')
    #if ionParam.calIonWithMerged == False:
    warningInfo = '{} calculating ionosphere swath by swath, there may be slight phase error between subswaths\n'.format(datetime.datetime.now())
    with open(os.path.join(ionParam.ionDirname, ionParam.warning), 'a') as f:
        f.write(warningInfo)

    #get bursts
    numValidSwaths = 0
    swathList = self._insar.getValidSwathList(self.swaths)
    for swath in swathList:
        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
        if minBurst==maxBurst:
            #print('Skipping processing of swath {0}'.format(swath))
            continue
        numValidSwaths += 1

    if numValidSwaths <= 1:
        raise Exception('There are less than one subswaths, no need to use swath-by-swath method to compute ionosphere!')
    else:
        xmlDirname = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname)
        (box, burstValidBox, burstValidBox2, frames) = getMergeBox(self, xmlDirname, numberRangeLooks=ionParam.numberRangeLooks, numberAzimuthLooks=ionParam.numberAzimuthLooks)

    #compute ionosphere swath by swath
    corList = []
    ampList = []
    ionosList = []
    nswath = len(swathList)
    ii = -1
    for i in range(nswath):
        swath = swathList[i]
        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
        if minBurst==maxBurst:
            print('Skipping processing of swath {0}'.format(swath))
            continue
        else:
            ii += 1

        ########################################################
        #STEP 1. MERGE THE BURSTS OF A SWATH
        ########################################################
        dirs = [ionParam.lowerDirname, ionParam.upperDirname]
        for dirx in dirs:
            outputFilename = self._insar.mergedIfgname
            outputDirname = os.path.join(ionParam.ionDirname, dirx, ionParam.mergedDirname + '_IW{0}'.format(swath))
            os.makedirs(outputDirname, exist_ok=True)
            suffix = '.full'
            if (numberRangeLooks0 == 1) and (numberAzimuthLooks0 == 1):
                suffix=''

            #merge
            burstPattern = 'burst_%02d.int'
            burstDirname = os.path.join(ionParam.ionDirname, dirx, ionParam.fineIfgDirname)
            ifg = self._insar.loadProduct( os.path.join(burstDirname, 'IW{0}.xml'.format(swath)))
            bst = [os.path.join(burstDirname, 'IW{0}'.format(swath), burstPattern%(x+1)) for x in range(minBurst, maxBurst)]
            #doing adjustment before use
            adjustValidWithLooks([ifg], box, numberAzimuthLooks, numberRangeLooks, edge=0, avalid='strict', rvalid=np.int(np.around(numberRangeLooks/8.0)))
            mergeBurstsVirtual([ifg], [bst], box, os.path.join(outputDirname, outputFilename+suffix))

            #take looks
            if suffix not in ['', None]:
                multilook2(os.path.join(outputDirname, outputFilename+suffix),
                           os.path.join(outputDirname, outputFilename),
                           numberAzimuthLooks0,
                           numberRangeLooks0)
            else:
                print('skipping multilooking')

        #The orginal coherence calculated by topsApp.py is not good at all, use the following coherence instead
        lowerintfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname + '_IW{0}'.format(swath), self._insar.mergedIfgname)
        upperintfile = os.path.join(ionParam.ionDirname, ionParam.upperDirname, ionParam.mergedDirname + '_IW{0}'.format(swath), self._insar.mergedIfgname)
        corfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname + '_IW{0}'.format(swath), self._insar.correlationFilename)

        img = isceobj.createImage()
        img.load(lowerintfile + '.xml')
        width = img.width
        length = img.length
        lowerint = np.fromfile(lowerintfile, dtype=np.complex64).reshape(length, width)
        upperint = np.fromfile(upperintfile, dtype=np.complex64).reshape(length, width)


        ##########################################################################
        #slight filtering to improve the estimation accurary of swath difference
        if 1 and shutil.which('psfilt1') != None:
            cmd1 = 'mv {} tmp'.format(lowerintfile)
            cmd2 = 'psfilt1 tmp {} {} .3 32 8'.format(lowerintfile, width)
            cmd3 = 'rm tmp'
            cmd4 = 'mv {} tmp'.format(upperintfile)
            cmd5 = 'psfilt1 tmp {} {} .3 32 8'.format(upperintfile, width)
            cmd6 = 'rm tmp'

            runCmd(cmd1)
            runCmd(cmd2)
            runCmd(cmd3)
            runCmd(cmd4)
            runCmd(cmd5)
            runCmd(cmd6)
        ##########################################################################


        #compute coherence only using interferogram
        #here I use differential interferogram of lower and upper band interferograms
        #so that coherence is not affected by fringes
        cord = cal_coherence(lowerint*np.conjugate(upperint), win=3, edge=4)
        cor = np.zeros((length*2, width), dtype=np.float32)
        cor[0:length*2:2, :] = np.sqrt( (np.absolute(lowerint)+np.absolute(upperint))/2.0 )
        cor[1:length*2:2, :] = cord
        cor.astype(np.float32).tofile(corfile)

        #create xml and vrt
        #img.scheme = 'BIL'
        #img.bands = 2
        #img.filename = corfile
        #img.renderHdr()
        
        #img = isceobj.Image.createUnwImage()
        img = isceobj.createOffsetImage()
        img.setFilename(corfile)
        img.extraFilename = corfile + '.vrt'
        img.setWidth(width)
        img.setLength(length)
        img.renderHdr()

        ########################################################
        #STEP 2. UNWRAP SWATH INTERFEROGRAM
        ########################################################
        dirs = [ionParam.lowerDirname, ionParam.upperDirname]
        #there is only one coherence file in lower directory
        corfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname + '_IW{0}'.format(swath), self._insar.correlationFilename)
        for dirx in dirs:
            procdir = os.path.join(ionParam.ionDirname, dirx, ionParam.mergedDirname + '_IW{0}'.format(swath))
            wrapName = os.path.join(procdir, self._insar.mergedIfgname)
            unwrapName = os.path.join(procdir, self._insar.unwrappedIntFilename)
            xmlDirname = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname)
            #unwrap
            snaphuUnwrap(self, xmlDirname, wrapName, corfile, unwrapName, numberRangeLooks0, numberAzimuthLooks0, costMode = 'SMOOTH',initMethod = 'MCF', defomax = 2, initOnly = True)
            #remove wired things in no-data area
            maskUnwrap(unwrapName, wrapName)

        if [ionParam.numberRangeLooks0, ionParam.numberAzimuthLooks0] != [ionParam.numberRangeLooks, ionParam.numberAzimuthLooks]:
            multilook_unw(self, ionParam, ionParam.mergedDirname + '_IW{0}'.format(swath))

        ########################################################
        #STEP 3. COMPUTE IONOSPHERE
        ########################################################
        #get files
        lowerUnwfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname + '_IW{0}'.format(swath), self._insar.unwrappedIntFilename)
        upperUnwfile = os.path.join(ionParam.ionDirname, ionParam.upperDirname, ionParam.mergedDirname + '_IW{0}'.format(swath), self._insar.unwrappedIntFilename)
        corfile = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.mergedDirname + '_IW{0}'.format(swath), self._insar.correlationFilename)

        #use image size from lower unwrapped interferogram
        img = isceobj.createImage()
        img.load(lowerUnwfile + '.xml')
        width = img.width
        length = img.length

        lowerUnw = (np.fromfile(lowerUnwfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
        upperUnw = (np.fromfile(upperUnwfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
        lowerAmp = (np.fromfile(lowerUnwfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
        upperAmp = (np.fromfile(upperUnwfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
        cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
        amp = np.sqrt(lowerAmp**2+upperAmp**2)

        #masked out user-specified areas
        if ionParam.maskedAreas != None:
            maskedAreas = reformatMaskedAreas(ionParam.maskedAreas, length, width)
            for area in maskedAreas:
                lowerUnw[area[0]:area[1], area[2]:area[3]] = 0
                upperUnw[area[0]:area[1], area[2]:area[3]] = 0
                cor[area[0]:area[1], area[2]:area[3]] = 0

        #compute ionosphere
        fl = SPEED_OF_LIGHT / ionParam.radarWavelengthLower
        fu = SPEED_OF_LIGHT / ionParam.radarWavelengthUpper
        adjFlag = 1
        ionos = computeIonosphere(lowerUnw, upperUnw, cor, fl, fu, adjFlag, corThresholdAdj, 0)

        #dump result
        outDir = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname + '_IW{0}'.format(swath))
        os.makedirs(outDir, exist_ok=True)
        outFilename = os.path.join(outDir, ionParam.ionRawNoProj)
        ion = np.zeros((length*2, width), dtype=np.float32)
        ion[0:length*2:2, :] = amp
        ion[1:length*2:2, :] = ionos
        ion.astype(np.float32).tofile(outFilename)
        img.filename = outFilename
        img.extraFilename = outFilename + '.vrt'
        img.renderHdr()

        corList.append(cor)
        ampList.append(amp)
        ionosList.append(ionos)

    #do adjustment between ajacent swaths
    if numValidSwaths == 3:
        adjustList = [ionosList[0], ionosList[2]]
    else:
        adjustList = [ionosList[0]]
    for adjdata in adjustList:
        index = np.nonzero((adjdata!=0) * (ionosList[1]!=0) * (corList[1] > corThresholdSwathAdj))
        if index[0].size < 5:
            print('WARNING: too few samples available for adjustment between swaths: {} with coherence threshold: {}'.format(index[0].size, corThresholdSwathAdj))
            print('         no adjustment made')
            print('         to do ajustment, please consider using lower coherence threshold')
        else:
            print('number of samples available for adjustment in the overlap area: {}'.format(index[0].size))
            #diff = np.mean((ionosList[1] - adjdata)[index], dtype=np.float64)
            
            #use weighted mean instead
            wgt = corList[1][index]**14
            diff = np.sum((ionosList[1] - adjdata)[index] * wgt / np.sum(wgt, dtype=np.float64), dtype=np.float64)

            index2 = np.nonzero(adjdata!=0)
            adjdata[index2] = adjdata[index2] + diff

    #get merged ionosphere
    ampMerged = np.zeros((length, width), dtype=np.float32)
    corMerged = np.zeros((length, width), dtype=np.float32)
    ionosMerged = np.zeros((length, width), dtype=np.float32)
    for i in range(numValidSwaths):
        nBurst = len(burstValidBox[i])
        for j in range(nBurst):

            #index after multi-looking in merged image, index starts from 1
            first_line = np.int(np.around((burstValidBox[i][j][0] - 1) / numberAzimuthLooks + 1))
            last_line = np.int(np.around(burstValidBox[i][j][1] / numberAzimuthLooks))
            first_sample = np.int(np.around((burstValidBox[i][j][2] - 1) / numberRangeLooks + 1))
            last_sample = np.int(np.around(burstValidBox[i][j][3] / numberRangeLooks))

            corMerged[first_line-1:last_line-1+1, first_sample-1:last_sample-1+1] = \
                corList[i][first_line-1:last_line-1+1, first_sample-1:last_sample-1+1]

            ampMerged[first_line-1:last_line-1+1, first_sample-1:last_sample-1+1] = \
                ampList[i][first_line-1:last_line-1+1, first_sample-1:last_sample-1+1]

            ionosMerged[first_line-1:last_line-1+1, first_sample-1:last_sample-1+1] = \
                ionosList[i][first_line-1:last_line-1+1, first_sample-1:last_sample-1+1]

    #remove an empirical ramp
    if ionParam.rampRemovel != 0:
        warningInfo = '{} calculating ionosphere for cross S-1A/B interferogram, an empirical ramp is removed from estimated ionosphere\n'.format(datetime.datetime.now())
        with open(os.path.join(ionParam.ionDirname, ionParam.warning), 'a') as f:
            f.write(warningInfo)

        abramp = cal_cross_ab_ramp(swathList, box[1], numberRangeLooks, ionParam.passDirection)
        if ionParam.rampRemovel == -1:
            abramp *= -1.0
        #currently do not apply this
        #ionosMerged -= abramp[None, :]

    #dump ionosphere
    outDir = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname)
    os.makedirs(outDir, exist_ok=True)
    outFilename = os.path.join(outDir, ionParam.ionRawNoProj)
    ion = np.zeros((length*2, width), dtype=np.float32)
    ion[0:length*2:2, :] = ampMerged
    ion[1:length*2:2, :] = ionosMerged
    ion.astype(np.float32).tofile(outFilename)
    img.filename = outFilename
    img.extraFilename = outFilename + '.vrt'
    img.renderHdr()

    #dump coherence
    outFilename = os.path.join(outDir, ionParam.ionCorNoProj)
    ion[1:length*2:2, :] = corMerged
    ion.astype(np.float32).tofile(outFilename)
    img.filename = outFilename
    img.extraFilename = outFilename + '.vrt'
    img.renderHdr()


def multilookIndex(first, last, nl):
    '''
    create the index after multilooking
    the orginal 1-look index can start from any number such as 0, 1 or other number
    after multilooking, the index still starts from the same number.
    first: index of first pixel in the original 1-look array
    last: index of last pixel in the original 1-look array
    nl: number of looks(nl can also be 1). nl >= 1
    '''

    #number of pixels after multilooking
    num = int((last - first + 1)/nl)
    offset = (first + (first + nl - 1)) / 2.0
    index = offset + np.arange(num) * nl

    return index


def computeDopplerOffset(burst, firstline, lastline, firstcolumn, lastcolumn, nrlks=1, nalks=1):
    '''
    compute offset corresponding to center Doppler frequency
    firstline, lastline, firstcolumn, lastcolumn: index of original 1-look burst, index starts from 1.

    output: first lines > 0, last lines < 0
    '''

    Vs = np.linalg.norm(burst.orbit.interpolateOrbit(burst.sensingMid, method='hermite').getVelocity())
    Ks =   2 * Vs * burst.azimuthSteeringRate / burst.radarWavelength 

    #firstcolumn, lastcolumn: index starts from 1
    rng = multilookIndex(firstcolumn-1, lastcolumn-1, nrlks) * burst.rangePixelSize + burst.startingRange
    #firstline, lastline: index starts from 1
    eta = (  multilookIndex(firstline-1, lastline-1, nalks) - (burst.numberOfLines-1.0)/2.0) * burst.azimuthTimeInterval

    f_etac = burst.doppler(rng)
    Ka     = burst.azimuthFMRate(rng)
    eta_ref = (burst.doppler(burst.startingRange) / burst.azimuthFMRate(burst.startingRange) ) - (f_etac / Ka)
    Kt = Ks / (1.0 - Ks/Ka)

    #carr = np.pi * Kt[None,:] * ((eta[:,None] - eta_ref[None,:])**2)
    #center doppler frequency due to rotation
    dopplerOffset1 = (eta[:,None] - eta_ref[None,:]) * Kt / Ka[None,:] / (burst.azimuthTimeInterval * nalks)
    #center doppler frequency due to squint
    dopplerOffset2 = (f_etac[None,:] / Ka[None,:]) / (burst.azimuthTimeInterval * nalks)
    dopplerOffset = dopplerOffset1 + dopplerOffset2
 
    return (dopplerOffset, Ka)


def grd2ion(self, ionParam):
    from scipy import interpolate
    from scipy.interpolate import interp1d

    print('resampling ionosphere from ground to ionospheric layer')
    #get files
    corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCorNoProj)
    ionfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionRawNoProj)

    #use image size from lower unwrapped interferogram
    img = isceobj.createImage()
    img.load(corfile + '.xml')
    width = img.width
    length = img.length

    cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
    ionos = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]

    #use the satellite height of the mid burst of first swath of reference acquistion
    swathList = self._insar.getValidSwathList(self.swaths)
    reference = self._insar.loadProduct( os.path.join(self._insar.referenceSlcProduct, 'IW{0}.xml'.format(swathList[0])))
    minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swathList[0]-1)
    #no problem with this index at all
    midBurst = np.int(np.around((minBurst+ maxBurst-1) / 2.0))
    masBurst = reference.bursts[midBurst]
    #satellite height
    satHeight = np.linalg.norm(masBurst.orbit.interpolateOrbit(masBurst.sensingMid, method='hermite').getPosition())
    #orgininal doppler offset should be multiplied by this ratio
    ratio = ionParam.ionHeight/(satHeight-ionParam.earthRadius)

    xmlDirname = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname)
    (box, burstValidBox, burstValidBox2, frames) = getMergeBox(self, xmlDirname, numberRangeLooks=ionParam.numberRangeLooks, numberAzimuthLooks=ionParam.numberAzimuthLooks)

##############################################################################################################
    swathList = self._insar.getValidSwathList(self.swaths)
    frames=[]
    #for valid swaths and bursts, consistent with runMergeBursts.py
    for swath in swathList:
        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)

        if minBurst==maxBurst:
            print('Skipping processing of swath {0}'.format(swath))
            continue

        ifg = self._insar.loadProduct( os.path.join(xmlDirname, 'IW{0}.xml'.format(swath)))
        frames.append(ifg)
##############################################################################################################

    for band in [amp, ionos, cor]:
        nswath = len(frames)
        for i in range(nswath):
            nburst = len(frames[i].bursts)
            for j in range(nburst):
                #according to runBurstIfg.py, this is originally from self._insar.referenceSlcProduct, 'IW{0}.xml'
                masBurst = frames[i].bursts[j]
                (dopplerOffset, Ka) = computeDopplerOffset(masBurst, burstValidBox2[i][j][0], burstValidBox2[i][j][1], burstValidBox2[i][j][2], burstValidBox2[i][j][3], nrlks=ionParam.numberRangeLooks, nalks=ionParam.numberAzimuthLooks)
                offset = ratio * dopplerOffset

                #   0              1               2              3
                #firstlineAdj, lastlineAdj, firstcolumnAdj, lastcolumnAdj, 
                #after multiplication, index starts from 1
                firstline = np.int(np.around((burstValidBox[i][j][0] - 1) / ionParam.numberAzimuthLooks + 1))
                lastline = np.int(np.around(burstValidBox[i][j][1] / ionParam.numberAzimuthLooks))
                firstcolumn = np.int(np.around((burstValidBox[i][j][2] - 1) / ionParam.numberRangeLooks + 1))
                lastcolumn = np.int(np.around(burstValidBox[i][j][3] / ionParam.numberRangeLooks))

                #extract image
                burstImage = band[firstline-1:lastline, firstcolumn-1:lastcolumn]
                blength = lastline - firstline + 1
                bwidth = lastcolumn - firstcolumn + 1

                #interpolation
                index0 = np.linspace(0, blength-1, num=blength, endpoint=True)
                for k in range(bwidth):
                    index = index0 + offset[:, k]
                    value = burstImage[:, k]
                    f = interp1d(index, value, kind='cubic', fill_value="extrapolate")
                    
                    index_min = np.int(np.around(np.amin(index)))
                    index_max = np.int(np.around(np.amax(index)))
                    flag = index0 * 0.0
                    flag[index_min:index_max+1] = 1.0
                    #replace the original column with new column in burstImage
                    #this should also replace teh original column with new column in band
                    burstImage[:, k] = (f(index0)) * flag

    #dump ionosphere with projection
    outDir = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname)
    outFilename = os.path.join(outDir, ionParam.ionRaw)
    ion = np.zeros((length*2, width), dtype=np.float32)
    ion[0:length*2:2, :] = amp
    ion[1:length*2:2, :] = ionos
    ion.astype(np.float32).tofile(outFilename)
    img.filename = outFilename
    img.extraFilename = outFilename + '.vrt'
    img.renderHdr()

    #dump coherence with projection
    outFilename = os.path.join(outDir, ionParam.ionCor)
    ion[1:length*2:2, :] = cor
    ion.astype(np.float32).tofile(outFilename)
    img.filename = outFilename
    img.extraFilename = outFilename + '.vrt'
    img.renderHdr()


def gaussian(size, sigma, scale = 1.0):

    if size % 2 != 1:
        raise Exception('size must be odd')
    hsize = (size - 1) / 2
    x = np.arange(-hsize, hsize + 1) * scale
    f = np.exp(-x**2/(2.0*sigma**2)) / (sigma * np.sqrt(2.0*np.pi))
    f2d=np.matlib.repmat(f, size, 1) * np.matlib.repmat(f.reshape(size, 1), 1, size)

    return f2d/np.sum(f2d)


def adaptive_gaussian(ionos, wgt, size_max, size_min):
    '''
    This program performs Gaussian filtering with adaptive window size.
    ionos: ionosphere
    wgt: weight
    size_max: maximum window size
    size_min: minimum window size
    '''
    import scipy.signal as ss

    length = (ionos.shape)[0]
    width = (ionos.shape)[1]
    flag = (ionos!=0) * (wgt!=0)
    ionos *= flag
    wgt *= flag

    size_num = 100
    size = np.linspace(size_min, size_max, num=size_num, endpoint=True)
    std = np.zeros((length, width, size_num))
    flt = np.zeros((length, width, size_num))
    out = np.zeros((length, width, 1))

    #calculate filterd image and standard deviation
    #sigma of window size: size_max
    sigma = size_max / 2.0
    for i in range(size_num):
        size2 = np.int(np.around(size[i]))
        if size2 % 2 == 0:
            size2 += 1
        if (i+1) % 10 == 0:
            print('min win: %4d, max win: %4d, current win: %4d'%(np.int(np.around(size_min)), np.int(np.around(size_max)), size2))
        g2d = gaussian(size2, sigma*size2/size_max, scale=1.0)
        scale = ss.fftconvolve(wgt, g2d, mode='same')
        flt[:, :, i] = ss.fftconvolve(ionos*wgt, g2d, mode='same') / (scale + (scale==0))
        #variance of resulting filtered sample
        scale = scale**2
        var = ss.fftconvolve(wgt, g2d**2, mode='same') / (scale + (scale==0))
        #in case there is a large area without data where scale is very small, which leads to wired values in variance
        var[np.nonzero(var<0)] = 0
        std[:, :, i] = np.sqrt(var)

    std_mv = np.mean(std[np.nonzero(std!=0)], dtype=np.float64)
    diff_max = np.amax(np.absolute(std - std_mv)) + std_mv + 1
    std[np.nonzero(std==0)] = diff_max
    
    index = np.nonzero(np.ones((length, width))) + ((np.argmin(np.absolute(std - std_mv), axis=2)).reshape(length*width), )
    out = flt[index]
    out = out.reshape((length, width))

    #remove artifacts due to varying wgt
    size_smt = size_min
    if size_smt % 2 == 0:
        size_smt += 1
    g2d = gaussian(size_smt, size_smt/2.0, scale=1.0)
    scale = ss.fftconvolve((out!=0), g2d, mode='same')
    out2 = ss.fftconvolve(out, g2d, mode='same') / (scale + (scale==0))

    return out2


def filt_gaussian(self, ionParam):
    '''
    This function filters image using gaussian filter

    we projected the ionosphere value onto the ionospheric layer, and the indexes are integers.
    this reduces the number of samples used in filtering
    a better method is to project the indexes onto the ionospheric layer. This way we have orginal
    number of samples used in filtering. but this requries more complicated operation in filtering
    currently not implemented.
    a less accurate method is to use ionsphere without any projection
    '''

    #################################################
    #SET PARAMETERS HERE
    #if applying polynomial fitting
    #False: no fitting, True: with fitting
    fit = ionParam.ionFit
    #gaussian filtering window size
    size_max = ionParam.ionFilteringWinsizeMax
    size_min = ionParam.ionFilteringWinsizeMin

    #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self)
    corThresholdIon = 0.85
    #################################################

    print('filtering ionosphere')
    #I find it's better to use ionosphere that is not projected, it's mostly slowlying changing anyway.
    #this should also be better for operational use.
    ionfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionRawNoProj)
    #since I decide to use ionosphere that is not projected, I should also use coherence that is not projected.
    corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCorNoProj)

    #use ionosphere and coherence that are projected.
    #ionfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionRaw)
    #corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCor)

    outfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionFilt)

    img = isceobj.createImage()
    img.load(ionfile + '.xml')
    width = img.width
    length = img.length
    ion = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]

    #masked out user-specified areas
    if ionParam.maskedAreas != None:
        maskedAreas = reformatMaskedAreas(ionParam.maskedAreas, length, width)
        for area in maskedAreas:
            ion[area[0]:area[1], area[2]:area[3]] = 0
            cor[area[0]:area[1], area[2]:area[3]] = 0

    ########################################################################################
    #AFTER COHERENCE IS RESAMPLED AT grd2ion, THERE ARE SOME WIRED VALUES
    cor[np.nonzero(cor<0)] = 0.0
    cor[np.nonzero(cor>1)] = 0.0
    ########################################################################################

    ion_fit = weight_fitting(ion, cor, width, length, 1, 1, 1, 1, 2, corThresholdIon)

    #no fitting
    if fit == False:
        ion_fit *= 0

    ion -= ion_fit * (ion!=0)
    
    #minimize the effect of low coherence pixels
    #cor[np.nonzero( (cor<0.85)*(cor!=0) )] = 0.00001
    #filt = adaptive_gaussian(ion, cor, size_max, size_min)
    #cor**14 should be a good weight to use. 22-APR-2018
    filt = adaptive_gaussian(ion, cor**14, size_max, size_min)

    filt += ion_fit * (filt!=0)

    ion = np.zeros((length*2, width), dtype=np.float32)
    ion[0:length*2:2, :] = amp
    ion[1:length*2:2, :] = filt
    ion.astype(np.float32).tofile(outfile)
    img.filename = outfile
    img.extraFilename = outfile + '.vrt'
    img.renderHdr()


def ionosphere_shift(self, ionParam):
    '''
    calculate azimuth shift caused by ionosphere using ionospheric phase
    '''

    #################################################
    #SET PARAMETERS HERE
    #gaussian filtering window size
    #size = np.int(np.around(width / 12.0))
    #size = ionParam.ionshiftFilteringWinsize
    size_max = ionParam.ionshiftFilteringWinsizeMax
    size_min = ionParam.ionshiftFilteringWinsizeMin

    #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self)
    #if applying polynomial fitting
    #0: no fitting, 1: with fitting
    fit = 0
    corThresholdIonshift = 0.85
    #################################################


####################################################################
    #STEP 1. GET DERIVATIVE OF IONOSPHERE
####################################################################

    #get files
    ionfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionFilt)
    #we are using filtered ionosphere, so we should use coherence file that is not projected.
    #corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCor)
    corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCorNoProj)
    img = isceobj.createImage()
    img.load(ionfile + '.xml')
    width = img.width
    length = img.length
    amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
    ion = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]

    ########################################################################################
    #AFTER COHERENCE IS RESAMPLED AT grd2ion, THERE ARE SOME WIRED VALUES
    cor[np.nonzero(cor<0)] = 0.0
    cor[np.nonzero(cor>1)] = 0.0
    ########################################################################################

    #get the azimuth derivative of ionosphere
    dion = np.diff(ion, axis=0)
    dion = np.concatenate((dion, np.zeros((1,width))), axis=0)

    #remove the samples affected by zeros
    flag_ion0 = (ion!=0)
    #moving down by one line
    flag_ion1 = np.roll(flag_ion0, 1, axis=0)
    flag_ion1[0,:] = 0
    #moving up by one line
    flag_ion2 = np.roll(flag_ion0, -1, axis=0)
    flag_ion2[-1,:] = 0
    #now remove the samples affected by zeros
    flag_ion = flag_ion0 * flag_ion1 * flag_ion2
    dion *= flag_ion

    flag = flag_ion * (cor>corThresholdIonshift)
    index = np.nonzero(flag)


####################################################################
    #STEP 2. FIT A POLYNOMIAL TO THE DERIVATIVE OF IONOSPHERE
####################################################################

    order = 3

    #look for data to use
    point_index = np.nonzero(flag)
    m = point_index[0].shape[0]

    #calculate input index matrix
    x0=np.matlib.repmat(np.arange(width), length, 1)
    y0=np.matlib.repmat(np.arange(length).reshape(length, 1), 1, width)

    x = x0[point_index].reshape(m, 1)
    y = y0[point_index].reshape(m, 1)
    z = dion[point_index].reshape(m, 1)
    w = cor[point_index].reshape(m, 1)

    #convert to higher precision type before use
    x=np.asfarray(x,np.float64)
    y=np.asfarray(y,np.float64)
    z=np.asfarray(z,np.float64)
    w=np.asfarray(w,np.float64)
    coeff = fit_surface(x, y, z, w, order)

    rgindex = np.arange(width)
    azindex = np.arange(length).reshape(length, 1)
    #convert to higher precision type before use
    rgindex=np.asfarray(rgindex,np.float64)
    azindex=np.asfarray(azindex,np.float64)
    dion_fit = cal_surface(rgindex, azindex, coeff, order)

    #no fitting
    if fit == 0:
        dion_fit *= 0
    dion_res = (dion - dion_fit)*(dion!=0)


####################################################################
    #STEP 3. FILTER THE RESIDUAL OF THE DERIVATIVE OF IONOSPHERE
####################################################################

    #this will be affected by low coherence areas like water, so not use this.
    #filter the derivation of ionosphere
    #if size % 2 == 0:
    #    size += 1
    #sigma = size / 2.0

    #g2d = gaussian(size, sigma, scale=1.0)
    #scale = ss.fftconvolve((dion_res!=0), g2d, mode='same')
    #dion_filt = ss.fftconvolve(dion_res, g2d, mode='same') / (scale + (scale==0))

    #minimize the effect of low coherence pixels
    cor[np.nonzero( (cor<0.85)*(cor!=0) )] = 0.00001
    dion_filt = adaptive_gaussian(dion_res, cor, size_max, size_min)

    dion = (dion_fit + dion_filt)*(dion!=0)

    #return dion


####################################################################
    #STEP 4. CONVERT TO AZIMUTH SHIFT
####################################################################

    #use the satellite height of the mid burst of first swath of reference acquistion
    swathList = self._insar.getValidSwathList(self.swaths)
    reference = self._insar.loadProduct( os.path.join(self._insar.referenceSlcProduct, 'IW{0}.xml'.format(swathList[0])))
    minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swathList[0]-1)
    #no problem with this index at all
    midBurst = np.int(np.around((minBurst+ maxBurst-1) / 2.0))
    masBurst = reference.bursts[midBurst]

    #shift casued by ionosphere [unit: masBurst.azimuthTimeInterval]
    rng = masBurst.rangePixelSize * ((np.arange(width))*ionParam.numberRangeLooks + (ionParam.numberRangeLooks - 1.0) / 2.0) + masBurst.startingRange
    Ka = masBurst.azimuthFMRate(rng)
    ionShift = dion / (masBurst.azimuthTimeInterval * ionParam.numberAzimuthLooks) / (4.0 * np.pi) / Ka[None, :] / masBurst.azimuthTimeInterval

    #output
    outfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionShift)
    tmp = np.zeros((length*2, width), dtype=np.float32)
    tmp[0:length*2:2, :] = amp
    tmp[1:length*2:2, :] = ionShift
    tmp.astype(np.float32).tofile(outfile)
    img.filename = outfile
    img.extraFilename = outfile + '.vrt'
    img.renderHdr()


def ion2grd(self, ionParam):
    from scipy import interpolate
    from scipy.interpolate import interp1d

    #################################################
    #SET PARAMETERS HERE
    #correct phase error caused by non-zero center frequency
    #and azimuth shift caused by ionosphere
    #0: no correction
    #1: use mean value of a burst
    #2: use full burst
    azshiftFlag = ionParam.azshiftFlag
    #################################################

    print('resampling ionosphere from ionospheric layer to ground')
    #get files
    ionFiltFile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionFilt)
    dionfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionShift)
    corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCorNoProj)
    img = isceobj.createImage()
    img.load(ionFiltFile + '.xml')
    width = img.width
    length = img.length
    ion = (np.fromfile(ionFiltFile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    dion = (np.fromfile(dionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]

    print('resampling ionosphere in range')
    #in the following, column index of burst (one look) will never exceed merged image index (one look) on the left side.
    #so we only add one multi-looked sample on the right side in case it exceeds on this side
    #index starts from 0
    ionOneRangeLook = np.zeros((length, (width+1)*ionParam.numberRangeLooks), dtype=np.float32)
    if azshiftFlag == 2:
        dionOneRangeLook = np.zeros((length, (width+1)*ionParam.numberRangeLooks), dtype=np.float32)
    indexRange = np.linspace(1-1, (width+1)*ionParam.numberRangeLooks-1, num=(width+1)*ionParam.numberRangeLooks, endpoint=True)
    indexRange2 = multilookIndex(1-1, width*ionParam.numberRangeLooks-1, ionParam.numberRangeLooks)
    for i in range(length):
        f = interp1d(indexRange2, ion[i, :], kind='cubic', fill_value="extrapolate")
        ionOneRangeLook[i, :] = f(indexRange)
        if azshiftFlag == 2:
            f2 = interp1d(indexRange2, dion[i, :], kind='cubic', fill_value="extrapolate")
            dionOneRangeLook[i, :] = f2(indexRange)
 
    #use the satellite height of the mid burst of first swath of reference acquistion
    swathList = self._insar.getValidSwathList(self.swaths)
    reference = self._insar.loadProduct( os.path.join(self._insar.referenceSlcProduct, 'IW{0}.xml'.format(swathList[0])))
    minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swathList[0]-1)
    #no problem with this index at all
    midBurst = np.int(np.around((minBurst+ maxBurst-1) / 2.0))
    masBurst = reference.bursts[midBurst]
    #satellite height
    satHeight = np.linalg.norm(masBurst.orbit.interpolateOrbit(masBurst.sensingMid, method='hermite').getPosition())
    #orgininal doppler offset should be multiplied by this ratio
    ratio = ionParam.ionHeight/(satHeight-ionParam.earthRadius)

    xmlDirname = os.path.join(ionParam.ionDirname, ionParam.lowerDirname, ionParam.fineIfgDirname)
    (box, burstValidBox, burstValidBox2, frames) = getMergeBox(self, xmlDirname, numberRangeLooks=ionParam.numberRangeLooks, numberAzimuthLooks=ionParam.numberAzimuthLooks)

    ##############################################################################################################
    swathList = self._insar.getValidSwathList(self.swaths)
    frames=[]
    swathList2 = []
    minBurst2 =[]
    #for valid swaths and bursts, consistent with runMergeBursts.py
    for swath in swathList:
        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)

        if minBurst==maxBurst:
            print('Skipping processing of swath {0}'.format(swath))
            continue

        ifg = self._insar.loadProduct( os.path.join(xmlDirname, 'IW{0}.xml'.format(swath)))
        frames.append(ifg)
        swathList2.append(swath)
        minBurst2.append(minBurst)
    ##############################################################################################################

    print('resampling ionosphere in azimuth')
    nswath = len(frames)
    for i in range(nswath):
        nburst = len(frames[i].bursts)
        ###output directory for burst ionosphere
        outdir = os.path.join(ionParam.ionDirname, ionParam.ionBurstDirname, 'IW{0}'.format(swathList2[i]))
        os.makedirs(outdir, exist_ok=True)

        for j in range(nburst):
            #according to runBurstIfg.py, this is originally from self._insar.referenceSlcProduct, 'IW{0}.xml'
            masBurst = frames[i].bursts[j]
            (dopplerOffset, Ka) = computeDopplerOffset(masBurst, 1, masBurst.numberOfLines, 1, masBurst.numberOfSamples, nrlks=1, nalks=1)
            offset = ratio * dopplerOffset
            #output ionosphere for this burst
            burstIon = np.zeros((masBurst.numberOfLines, masBurst.numberOfSamples), dtype=np.float32)
            burstDion = np.zeros((masBurst.numberOfLines, masBurst.numberOfSamples), dtype=np.float32)

            #              index in merged           index in burst
            lineOff = burstValidBox[i][j][0] - burstValidBox2[i][j][0]
            columnOff = burstValidBox[i][j][2] - burstValidBox2[i][j][2]
            #use index starts from 0
            #1-look index of burst in the 1-look merged image
            indexBurst0 = np.linspace(0+lineOff, masBurst.numberOfLines-1+lineOff, num=masBurst.numberOfLines, endpoint=True)
            #1-look index of multi-looked merged image in the 1-look merged image
            indexMerged = multilookIndex(1-1, length*ionParam.numberAzimuthLooks-1, ionParam.numberAzimuthLooks)
            for k in range(masBurst.numberOfSamples):
                index = indexMerged
                value = ionOneRangeLook[:, k+columnOff]
                f = interp1d(index, value, kind='cubic', fill_value="extrapolate")

                indexBurst = indexBurst0 + offset[:, k]
                burstIon[:, k] = f(indexBurst)

                if azshiftFlag == 2:
                    value2 = dionOneRangeLook[:, k+columnOff]
                    f2 = interp1d(index, value2, kind='cubic', fill_value="extrapolate")
                    burstDion[:, k] = f2(indexBurst)

            #calculate phase caused by ionospheric shift and non-zero center frequency
            #index after multi-looking in merged image, index starts from 1
            first_line = np.int(np.around((burstValidBox[i][j][0] - 1) / ionParam.numberAzimuthLooks + 1))
            last_line = np.int(np.around(burstValidBox[i][j][1] / ionParam.numberAzimuthLooks))
            first_sample = np.int(np.around((burstValidBox[i][j][2] - 1) / ionParam.numberRangeLooks + 1))
            last_sample = np.int(np.around(burstValidBox[i][j][3] / ionParam.numberRangeLooks))

            burstDionMultilook = dion[first_line-1:last_line-1+1, first_sample-1:last_sample-1+1]
            #for avoid areas with strong decorrelation like water
            burstCorMultilook = cor[first_line-1:last_line-1+1, first_sample-1:last_sample-1+1]
            #index = np.nonzero(burstDionMultilook!=0)
            index = np.nonzero(burstCorMultilook>0.85)
            if len(index[0]) < 10:
                dionMean = 0.0
            else:
                dionMean = np.mean(burstDionMultilook[index], dtype=np.float64)

            if azshiftFlag == 0:
                #no correction
                burstIonShift = 0
            elif azshiftFlag == 1:
                #use mean value
                burstIonShift =  2.0 * np.pi * (dopplerOffset * Ka[None,:] * (masBurst.azimuthTimeInterval)) * (dionMean*masBurst.azimuthTimeInterval)
            elif azshiftFlag == 2:
                #use full burst
                burstIonShift =  2.0 * np.pi * (dopplerOffset * Ka[None,:] * (masBurst.azimuthTimeInterval)) * (burstDion*masBurst.azimuthTimeInterval)
            else:
                raise Exception('unknown option for correcting azimuth shift caused by ionosphere!')

            burstIon += burstIonShift
            print('resampling burst %02d of swath %d, azimuth shift caused by ionosphere: %8.5f azimuth lines'%(minBurst2[i]+j+1, swathList2[i], dionMean))

            #create xml and vrt files
            filename = os.path.join(outdir, '%s_%02d.ion'%('burst', minBurst2[i]+j+1))
            burstIon.astype(np.float32).tofile(filename)
            burstImg = isceobj.createImage()
            burstImg.setDataType('FLOAT')
            burstImg.setFilename(filename)
            burstImg.extraFilename = filename + '.vrt'
            burstImg.setWidth(masBurst.numberOfSamples)
            burstImg.setLength(masBurst.numberOfLines)
            burstImg.renderHdr()
        print('')


def multilook(data, nalks, nrlks):
    '''
    doing multiple looking
    
    ATTENTION:
    NO AVERAGING BY DIVIDING THE NUMBER OF TOTAL SAMPLES IS DONE.
    '''

    (length, width)=data.shape
    width2 = np.int(width/nrlks)
    length2 = np.int(length/nalks)

    tmp2 = np.zeros((length2, width), dtype=data.dtype)
    data2 = np.zeros((length2, width2), dtype=data.dtype)
    for i in range(nalks):
        tmp2 += data[i:length2*nalks:nalks, :]
    for i in range(nrlks):
        data2 += tmp2[:, i:width2*nrlks:nrlks]

    return data2


def get_overlap_box(swath, minBurst, maxBurst):

    #number of burst
    nBurst = maxBurst - minBurst
    if nBurst <= 1:
        print('number of burst: {}, no need to get overlap box'.format(nBurst))
        return None

    overlapBox = []
    overlapBox.append([])
    for ii in range(minBurst+1, maxBurst):
        topBurst = swath.bursts[ii-1]
        curBurst = swath.bursts[ii]

        #overlap lines, line index starts from 1
        offLine = np.int(np.round( (curBurst.sensingStart - topBurst.sensingStart).total_seconds() / curBurst.azimuthTimeInterval))
        firstLineTop = topBurst.firstValidLine + 1
        lastLineTop = topBurst.firstValidLine + topBurst.numValidLines
        firstLineCur = offLine + curBurst.firstValidLine + 1
        lastLineCur = offLine + curBurst.firstValidLine + curBurst.numValidLines

        if lastLineTop < firstLineCur:
            raise Exception('there is not enough overlap between burst {} and burst {}\n'.format(ii-1+1, ii+1))

        firstLine = firstLineCur
        lastLine = lastLineTop

        #overlap samples, sample index starts from 1
        offSample = np.int(np.round(       (curBurst.startingRange - topBurst.startingRange) / curBurst.rangePixelSize         ))
        firstSampleTop = topBurst.firstValidSample + 1
        lastSampleTop = topBurst.firstValidSample + topBurst.numValidSamples
        firstSampleCur = offSample + curBurst.firstValidSample + 1
        lastSampleCur = offSample + curBurst.firstValidSample + curBurst.numValidSamples

        firstSample = max(firstSampleTop, firstSampleCur)
        lastSample = min(lastSampleTop, lastSampleCur)

        #overlap area index. all indexes start from 1.
        #                 |                 top burst                  |                                   current burst                                 |
        #                    0           1          2           3              4                   5                    6                     7
        overlapBox.append([firstLine, lastLine, firstSample, lastSample, firstLine-offLine, lastLine-offLine, firstSample-offSample, lastSample-offSample])

    return overlapBox


def esd(self, ionParam):
    '''
    esd after ionosphere correction
    '''
    ######################################
    #SET PARAMETERS HERE
    #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self)
    nalks = 5
    nrlks = 30
    corThreshold = 0.75
    ######################################

    print('applying ESD to compensate phase error caused by residual misregistration')

    virtual = self.useVirtualFiles
    swathList = self._insar.getValidSwathList(self.swaths)
    for swath in swathList:

        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
        nBurst = maxBurst - minBurst

        if nBurst <= 1:
            continue
    
        ####Load relevant products
        reference = self._insar.loadProduct( os.path.join(self._insar.referenceSlcProduct, 'IW{0}.xml'.format(swath)))
        secondary = self._insar.loadProduct( os.path.join(self._insar.fineCoregDirname, 'IW{0}.xml'.format(swath)))

        #get overlap area
        for ii in range(minBurst, maxBurst):
            jj = ii - minBurst
            ####Process the top bursts
            masBurst = reference.bursts[ii] 
            slvBurst = secondary.bursts[jj]
            adjustValidLineSample(masBurst,slvBurst)
        overlapBox = get_overlap_box(reference, minBurst, maxBurst)
        
        #using esd to calculate mis-registration
        misreg = np.array([])
        totalSamples = 0
        for ii in range(minBurst+1, maxBurst):
            jj = ii - minBurst
            ####Process the top bursts
            masBurstTop = reference.bursts[ii-1] 
            slvBurstTop = secondary.bursts[jj-1]

            masBurstCur = reference.bursts[ii] 
            slvBurstCur = secondary.bursts[jj]

            #get info
            referencename = masBurstTop.image.filename
            secondaryname = slvBurstTop.image.filename
            ionname = os.path.join(ionParam.ionDirname, ionParam.ionBurstDirname, 'IW{0}'.format(swath), '%s_%02d.ion'%('burst',ii+1-1))
            rngname = os.path.join(self._insar.fineOffsetsDirname, 'IW{0}'.format(swath), 'range_%02d.off'%(ii+1-1))
            fact = 4.0 * np.pi * slvBurstTop.rangePixelSize / slvBurstTop.radarWavelength
            #infTop = multiply2(referencename, secondaryname, ionname, rngname, fact, overlapBox[jj][0:4], virtual=virtual)
            infTop = multiply2(referencename, secondaryname, fact, rngname=rngname, ionname=ionname, infname=None, overlapBox=overlapBox[jj][0:4], valid=True, virtual=virtual)
            (dopTop, Ka) = computeDopplerOffset(masBurstTop, overlapBox[jj][0], overlapBox[jj][1], overlapBox[jj][2], overlapBox[jj][3], nrlks=nrlks, nalks=nalks)
            #rng = multilookIndex(overlapBox[jj][2]-1, overlapBox[jj][3]-1, nrlks) * masBurstTop.rangePixelSize + masBurstTop.startingRange
            #Ka  = masBurstTop.azimuthFMRate(rng)
            frqTop = dopTop * Ka[None,:] * (masBurstTop.azimuthTimeInterval * nalks)

            referencename = masBurstCur.image.filename
            secondaryname = slvBurstCur.image.filename
            ionname = os.path.join(ionParam.ionDirname, ionParam.ionBurstDirname, 'IW{0}'.format(swath), '%s_%02d.ion'%('burst',ii+1))
            rngname = os.path.join(self._insar.fineOffsetsDirname, 'IW{0}'.format(swath), 'range_%02d.off'%(ii+1))
            fact = 4.0 * np.pi * slvBurstCur.rangePixelSize / slvBurstCur.radarWavelength
            #infCur = multiply2(referencename, secondaryname, ionname, rngname, fact, overlapBox[jj][4:8], virtual=virtual)
            infCur = multiply2(referencename, secondaryname, fact, rngname=rngname, ionname=ionname, infname=None, overlapBox=overlapBox[jj][4:8], valid=True, virtual=virtual)
            (dopCur, Ka) = computeDopplerOffset(masBurstCur, overlapBox[jj][4], overlapBox[jj][5], overlapBox[jj][6], overlapBox[jj][7], nrlks=nrlks, nalks=nalks)
            #rng = multilookIndex(overlapBox[jj][6]-1, overlapBox[jj][7]-1, nrlks) * masBurstCur.rangePixelSize + masBurstCur.startingRange
            #Ka  = masBurstCur.azimuthFMRate(rng)
            frqCur = dopCur * Ka[None,:] * (masBurstCur.azimuthTimeInterval * nalks)

            infTop = multilook(infTop, nalks, nrlks)
            infCur = multilook(infCur, nalks, nrlks)
            infDif = infTop * np.conjugate(infCur)
            cor    = cal_coherence(infDif, win=3, edge=4)
            index = np.nonzero(cor > corThreshold)
            totalSamples += infTop.size

            if index[0].size:
                #misregistration in sec. it should be OK to only use reference frequency to compute ESD
                misreg0 = np.angle(infDif[index]) / (2.0 * np.pi * (frqTop[index]-frqCur[index]))
                misreg=np.append(misreg, misreg0.flatten())
                print("misregistration at burst %02d and burst %02d of swath %d: %10.5f azimuth lines"%(ii+1-1, ii+1, swath, np.mean(misreg0, dtype=np.float64)/masBurstCur.azimuthTimeInterval))
            else:
                print("no samples available for ESD at burst %02d and burst %02d of swath %d"%(ii+1-1, ii+1, swath))

        percentage = 100.0 * len(misreg) / totalSamples
        #number of samples per overlap: 100/5*23334/150 = 3111.2
        print("samples available for ESD at swath %d: %d out of %d available, percentage: %5.1f%%"%(swath, len(misreg), totalSamples, percentage))
        if len(misreg) < 1000:
            print("too few samples available for ESD, no ESD correction will be applied\n")
            misreg = 0
            continue
        else:
            misreg = np.mean(misreg, dtype=np.float64)
            print("misregistration from ESD: {} sec, {} azimuth lines\n".format(misreg, misreg/reference.bursts[minBurst].azimuthTimeInterval))

        #use mis-registration estimated from esd to compute phase error
        for ii in range(minBurst, maxBurst):
            jj = ii - minBurst
            ####Process the top bursts
            masBurst = reference.bursts[ii] 
            slvBurst = secondary.bursts[jj]

            ionname = os.path.join(ionParam.ionDirname, ionParam.ionBurstDirname, 'IW{0}'.format(swath), '%s_%02d.ion'%('burst',ii+1))
            ion = np.fromfile(ionname, dtype=np.float32).reshape(masBurst.numberOfLines, masBurst.numberOfSamples)
            (dopplerOffset, Ka) = computeDopplerOffset(masBurst, 1, masBurst.numberOfLines, 1, masBurst.numberOfSamples, nrlks=1, nalks=1)
            centerFrequency = dopplerOffset * Ka[None,:] * (masBurst.azimuthTimeInterval)

            ion += 2.0 * np.pi * centerFrequency * misreg
            #overwrite
            ion.astype(np.float32).tofile(ionname)


def esd_noion(self, ionParam):
    '''
    esd after ionosphere correction
    '''
    ######################################
    #SET PARAMETERS HERE
    #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self)
    nalks = 5
    nrlks = 30
    corThreshold = 0.75
    ######################################

    print('applying ESD to compensate phase error caused by residual misregistration')


    esddir = 'esd'


    virtual = self.useVirtualFiles
    swathList = self._insar.getValidSwathList(self.swaths)
    for swath in swathList:

        minBurst, maxBurst = self._insar.commonReferenceBurstLimits(swath-1)
        nBurst = maxBurst - minBurst

        if nBurst <= 1:
            continue
    
        ####Load relevant products
        reference = self._insar.loadProduct( os.path.join(self._insar.referenceSlcProduct, 'IW{0}.xml'.format(swath)))
        secondary = self._insar.loadProduct( os.path.join(self._insar.fineCoregDirname, 'IW{0}.xml'.format(swath)))

        #get overlap area
        for ii in range(minBurst, maxBurst):
            jj = ii - minBurst
            ####Process the top bursts
            masBurst = reference.bursts[ii] 
            slvBurst = secondary.bursts[jj]
            adjustValidLineSample(masBurst,slvBurst)
        overlapBox = get_overlap_box(reference, minBurst, maxBurst)
        
        #using esd to calculate mis-registration
        misreg = np.array([])
        totalSamples = 0
        for ii in range(minBurst+1, maxBurst):
            jj = ii - minBurst
            ####Process the top bursts
            masBurstTop = reference.bursts[ii-1] 
            slvBurstTop = secondary.bursts[jj-1]

            masBurstCur = reference.bursts[ii] 
            slvBurstCur = secondary.bursts[jj]

            #get info
            referencename = masBurstTop.image.filename
            secondaryname = slvBurstTop.image.filename
            ionname = os.path.join(ionParam.ionDirname, ionParam.ionBurstDirname, 'IW{0}'.format(swath), '%s_%02d.ion'%('burst',ii+1-1))
            rngname = os.path.join(self._insar.fineOffsetsDirname, 'IW{0}'.format(swath), 'range_%02d.off'%(ii+1-1))
            fact = 4.0 * np.pi * slvBurstTop.rangePixelSize / slvBurstTop.radarWavelength
            #infTop = multiply2(referencename, secondaryname, ionname, rngname, fact, overlapBox[jj][0:4], virtual=virtual)
            infTop = multiply2(referencename, secondaryname, fact, rngname=rngname, ionname=None, infname=None, overlapBox=overlapBox[jj][0:4], valid=True, virtual=virtual)
            (dopTop, Ka) = computeDopplerOffset(masBurstTop, overlapBox[jj][0], overlapBox[jj][1], overlapBox[jj][2], overlapBox[jj][3], nrlks=nrlks, nalks=nalks)
            #rng = multilookIndex(overlapBox[jj][2]-1, overlapBox[jj][3]-1, nrlks) * masBurstTop.rangePixelSize + masBurstTop.startingRange
            #Ka  = masBurstTop.azimuthFMRate(rng)
            frqTop = dopTop * Ka[None,:] * (masBurstTop.azimuthTimeInterval * nalks)

            referencename = masBurstCur.image.filename
            secondaryname = slvBurstCur.image.filename
            ionname = os.path.join(ionParam.ionDirname, ionParam.ionBurstDirname, 'IW{0}'.format(swath), '%s_%02d.ion'%('burst',ii+1))
            rngname = os.path.join(self._insar.fineOffsetsDirname, 'IW{0}'.format(swath), 'range_%02d.off'%(ii+1))
            fact = 4.0 * np.pi * slvBurstCur.rangePixelSize / slvBurstCur.radarWavelength
            #infCur = multiply2(referencename, secondaryname, ionname, rngname, fact, overlapBox[jj][4:8], virtual=virtual)
            infCur = multiply2(referencename, secondaryname, fact, rngname=rngname, ionname=None, infname=None, overlapBox=overlapBox[jj][4:8], valid=True, virtual=virtual)
            (dopCur, Ka) = computeDopplerOffset(masBurstCur, overlapBox[jj][4], overlapBox[jj][5], overlapBox[jj][6], overlapBox[jj][7], nrlks=nrlks, nalks=nalks)
            #rng = multilookIndex(overlapBox[jj][6]-1, overlapBox[jj][7]-1, nrlks) * masBurstCur.rangePixelSize + masBurstCur.startingRange
            #Ka  = masBurstCur.azimuthFMRate(rng)
            frqCur = dopCur * Ka[None,:] * (masBurstCur.azimuthTimeInterval * nalks)

            infTop = multilook(infTop, nalks, nrlks)
            infCur = multilook(infCur, nalks, nrlks)
            infDif = infTop * np.conjugate(infCur)
            cor    = cal_coherence(infDif, win=3, edge=4)
            index = np.nonzero(cor > corThreshold)
            totalSamples += infTop.size

            if index[0].size:
                #misregistration in sec. it should be OK to only use reference frequency to compute ESD
                misreg0 = np.angle(infDif[index]) / (2.0 * np.pi * (frqTop[index]-frqCur[index]))
                misreg=np.append(misreg, misreg0.flatten())
                print("misregistration at burst %02d and burst %02d of swath %d: %10.5f azimuth lines"%(ii+1-1, ii+1, swath, np.mean(misreg0, dtype=np.float64)/masBurstCur.azimuthTimeInterval))
            else:
                print("no samples available for ESD at burst %02d and burst %02d of swath %d"%(ii+1-1, ii+1, swath))

        percentage = 100.0 * len(misreg) / totalSamples
        #number of samples per overlap: 100/5*23334/150 = 3111.2
        print("samples available for ESD at swath %d: %d out of %d available, percentage: %5.1f%%"%(swath, len(misreg), totalSamples, percentage))
        if len(misreg) < 1000:
            print("too few samples available for ESD, no ESD correction will be applied\n")
            misreg = 0
            continue
        else:
            misreg = np.mean(misreg, dtype=np.float64)
            print("misregistration from ESD: {} sec, {} azimuth lines\n".format(misreg, misreg/reference.bursts[minBurst].azimuthTimeInterval))


        sdir = os.path.join(ionParam.ionDirname, esddir, 'IW{0}'.format(swath))
        os.makedirs(sdir, exist_ok=True)

        #use mis-registration estimated from esd to compute phase error
        for ii in range(minBurst, maxBurst):
            jj = ii - minBurst
            ####Process the top bursts
            masBurst = reference.bursts[ii] 
            slvBurst = secondary.bursts[jj]

            #ionname = os.path.join(ionParam.ionDirname, ionParam.ionBurstDirname, 'IW{0}'.format(swath), '%s_%02d.ion'%('burst',ii+1))
            #ion = np.fromfile(ionname, dtype=np.float32).reshape(masBurst.numberOfLines, masBurst.numberOfSamples)
            
            (dopplerOffset, Ka) = computeDopplerOffset(masBurst, 1, masBurst.numberOfLines, 1, masBurst.numberOfSamples, nrlks=1, nalks=1)
            centerFrequency = dopplerOffset * Ka[None,:] * (masBurst.azimuthTimeInterval)

            ion = 2.0 * np.pi * centerFrequency * misreg
            #overwrite
            ionname = os.path.join(ionParam.ionDirname, esddir, 'IW{0}'.format(swath), '%s_%02d.esd'%('burst',ii+1))
            ion.astype(np.float32).tofile(ionname)



            #create xml and vrt files
            burstImg = isceobj.createImage()
            burstImg.setDataType('FLOAT')
            burstImg.setFilename(ionname)
            burstImg.extraFilename = ionname + '.vrt'
            burstImg.setWidth(masBurst.numberOfSamples)
            burstImg.setLength(masBurst.numberOfLines)
            burstImg.renderHdr()


def rawion(self, ionParam):
    '''
    a simple wrapper
    '''

    if ionParam.calIonWithMerged == True:
        #merge bursts
        merge(self, ionParam)

        #unwrap
        unwrap(self, ionParam)

        #compute ionosphere
        ionosphere(self, ionParam)
    else:
        #an alternative of the above steps: processing swath by swath
        ionSwathBySwath(self, ionParam)


def run_step(currentStep, ionParam):
    return ionParam.allSteps.index(ionParam.startStep) <= ionParam.allSteps.index(currentStep) <= ionParam.allSteps.index(ionParam.endStep)


def runIon(self):

    #get processing parameters
    ionParam = setup(self)

    #if do ionospheric correction
    if ionParam.doIon == False:
        return

    #form subband interferograms
    if run_step('subband', ionParam):
        subband(self, ionParam)

    #compute ionosphere (raw_no_projection.ion) and coherence (raw_no_projection.cor) without projection
    if run_step('rawion', ionParam):
        rawion(self, ionParam)
    #next we move to 'ion_cal' to do the remaining processing

    #resample ionosphere from the ground layer to ionospheric layer
    if run_step('grd2ion', ionParam):
        grd2ion(self, ionParam)

    #filter ionosphere
    if run_step('filt_gaussian', ionParam):
        filt_gaussian(self, ionParam)

    #only do the following steps when considering burst properties
    #ionosphere shift
    if run_step('ionosphere_shift', ionParam) and ionParam.considerBurstProperties:
        ionosphere_shift(self, ionParam)

    #resample from ionospheric layer to ground layer, get ionosphere for each burst
    if run_step('ion2grd', ionParam) and ionParam.considerBurstProperties:
        ion2grd(self, ionParam)

    #esd
    if run_step('esd', ionParam) and ionParam.considerBurstProperties:
        esd(self, ionParam)

    #pure esd without applying ionospheric correction
    #esd_noion(self, ionParam)

    return