ISCE_INSAR/components/isceobj/Alos2Proc/runIonFilt.py

#
# Author: Cunren Liang
# Copyright 2015-present, NASA-JPL/Caltech
#

import os
import logging
import numpy as np
import numpy.matlib

import isceobj

logger = logging.getLogger('isce.alos2insar.runIonFilt')

def runIonFilt(self):
    '''compute and filter ionospheric phase
    '''
    catalog = isceobj.Catalog.createCatalog(self._insar.procDoc.name)
    self.updateParamemetersFromUser()

    if not self.doIon:
        catalog.printToLog(logger, "runIonFilt")
        self._insar.procDoc.addAllFromCatalog(catalog)
        return

    masterTrack = self._insar.loadTrack(master=True)
    slaveTrack = self._insar.loadTrack(master=False)

    from isceobj.Alos2Proc.runIonSubband import defineIonDir
    ionDir = defineIonDir()
    subbandPrefix = ['lower', 'upper']

    ionCalDir = os.path.join(ionDir['ion'], ionDir['ionCal'])
    os.makedirs(ionCalDir, exist_ok=True)
    os.chdir(ionCalDir)


    ############################################################
    # STEP 1. compute ionospheric phase
    ############################################################
    from isceobj.Constants import SPEED_OF_LIGHT
    from isceobj.Alos2Proc.Alos2ProcPublic import create_xml

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

    print('\ncomputing ionosphere')
    #get files
    ml2 = '_{}rlks_{}alks'.format(self._insar.numberRangeLooks1*self._insar.numberRangeLooksIon, 
                              self._insar.numberAzimuthLooks1*self._insar.numberAzimuthLooksIon)

    lowerUnwfile = subbandPrefix[0]+ml2+'.unw'
    upperUnwfile = subbandPrefix[1]+ml2+'.unw'
    corfile = 'diff'+ml2+'.cor'

    #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, :]
    cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
    #amp = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]

    #masked out user-specified areas
    if self.maskedAreasIon != None:
        maskedAreas = reformatMaskedAreas(self.maskedAreasIon, 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

    #remove possible wired values in coherence
    cor[np.nonzero(cor<0)] = 0.0
    cor[np.nonzero(cor>1)] = 0.0

    #compute ionosphere
    fl = SPEED_OF_LIGHT / self._insar.subbandRadarWavelength[0]
    fu = SPEED_OF_LIGHT / self._insar.subbandRadarWavelength[1]
    adjFlag = 1
    cor[np.nonzero(cor<corThresholdAdj)] = 0.0
    ionos = computeIonosphere(lowerUnw, upperUnw, cor**corOrderAdj, fl, fu, adjFlag, 0)

    #dump ionosphere
    ionfile = 'ion'+ml2+'.ion'
    # 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(ionfile)
    # img.filename = ionfile
    # img.extraFilename = ionfile + '.vrt'
    # img.renderHdr()

    ionos.astype(np.float32).tofile(ionfile)
    create_xml(ionfile, width, length, 'float')


    ############################################################
    # STEP 2. filter ionospheric phase
    ############################################################

    #################################################
    #SET PARAMETERS HERE
    #if applying polynomial fitting
    #False: no fitting, True: with fitting
    fit = self.fitIon
    #gaussian filtering window size
    size_max = self.filteringWinsizeMaxIon
    size_min = self.filteringWinsizeMinIon

    if size_min >= size_max:
        print('\n\nWARNING: minimum window size for filtering ionosphere phase {} >= maximum window size {}'.format(size_min, size_max))
        print('         resetting maximum window size to {}\n\n'.format(size_min+5))
        size_max = size_min + 5

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

    #Now changed to use lower band coherence. crl, 23-apr-2020.
    useDiffCoherence = False
    if useDiffCoherence:
        #parameters for using diff coherence
        corfile = 'diff'+ml2+'.cor'
        corThresholdFit = 0.95
        # 1 is not good for low coherence case, changed to 20
        #corOrderFit = 1
        corOrderFit = 20
        corOrderFilt = 20
    else:
        #parameters for using lower/upper band coherence
        corfile = subbandPrefix[0]+ml2+'.cor'
        corThresholdFit = 0.4
        corOrderFit = 10
        corOrderFilt = 10

    #################################################

    print('\nfiltering ionosphere')
    ionfile = 'ion'+ml2+'.ion'
    #corfile = 'diff'+ml2+'.cor'
    ionfiltfile = 'filt_ion'+ml2+'.ion'

    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, :]
    ion = np.fromfile(ionfile, dtype=np.float32).reshape(length, width)
    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 self.maskedAreasIon != None:
        maskedAreas = reformatMaskedAreas(self.maskedAreasIon, 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

    #remove possible wired values in coherence
    cor[np.nonzero(cor<0)] = 0.0
    cor[np.nonzero(cor>1)] = 0.0

    # #applying water body mask here
    # waterBodyFile = 'wbd'+ml2+'.wbd'
    # if os.path.isfile(waterBodyFile):
    #     print('applying water body mask to coherence used to compute ionospheric phase')
    #     wbd = np.fromfile(waterBodyFile, dtype=np.int8).reshape(length, width)
    #     cor[np.nonzero(wbd!=0)] = 0.00001

    if fit:
        import copy
        wgt = copy.deepcopy(cor)
        wgt[np.nonzero(wgt<corThresholdFit)] = 0.0
        ion_fit = weight_fitting(ion, wgt**corOrderFit, width, length, 1, 1, 1, 1, 2)
        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**corOrderFilt, size_max, size_min)

    if fit:
        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(ionfiltfile)
    # img.filename = ionfiltfile
    # img.extraFilename = ionfiltfile + '.vrt'
    # img.renderHdr()

    filt.astype(np.float32).tofile(ionfiltfile)
    create_xml(ionfiltfile, width, length, 'float')


    ############################################################
    # STEP 3. resample ionospheric phase
    ############################################################
    from contrib.alos2proc_f.alos2proc_f import rect
    from isceobj.Alos2Proc.Alos2ProcPublic import create_xml
    from scipy.interpolate import interp1d
    import shutil

    #################################################
    #SET PARAMETERS HERE
    #interpolation method
    interpolationMethod = 1
    #################################################

    print('\ninterpolate ionosphere')

    ml3 = '_{}rlks_{}alks'.format(self._insar.numberRangeLooks1*self._insar.numberRangeLooks2, 
                              self._insar.numberAzimuthLooks1*self._insar.numberAzimuthLooks2)

    ionfiltfile = 'filt_ion'+ml2+'.ion'
    #ionrectfile = 'filt_ion'+ml3+'.ion'
    ionrectfile = self._insar.multilookIon

    img = isceobj.createImage()
    img.load(ionfiltfile + '.xml')
    width2 = img.width
    length2 = img.length

    img = isceobj.createImage()
    img.load(os.path.join('../../', ionDir['insar'], self._insar.multilookDifferentialInterferogram) + '.xml')
    width3 = img.width
    length3 = img.length

    #number of range looks output
    nrlo = self._insar.numberRangeLooks1*self._insar.numberRangeLooks2
    #number of range looks input
    nrli = self._insar.numberRangeLooks1*self._insar.numberRangeLooksIon
    #number of azimuth looks output
    nalo = self._insar.numberAzimuthLooks1*self._insar.numberAzimuthLooks2
    #number of azimuth looks input
    nali = self._insar.numberAzimuthLooks1*self._insar.numberAzimuthLooksIon

    if (self._insar.numberRangeLooks2 != self._insar.numberRangeLooksIon) or \
       (self._insar.numberAzimuthLooks2 != self._insar.numberAzimuthLooksIon):
        #this should be faster using fortran
        if interpolationMethod == 0:
            rect(ionfiltfile, ionrectfile,
                width2,length2,
                width3,length3,
                nrlo/nrli, 0.0,
                0.0, nalo/nali,
                (nrlo-nrli)/(2.0*nrli),
                (nalo-nali)/(2.0*nali),
                'REAL','Bilinear')
        #finer, but slower method
        else:
            ionfilt = np.fromfile(ionfiltfile, dtype=np.float32).reshape(length2, width2)
            index2 = np.linspace(0, width2-1, num=width2, endpoint=True)
            index3 = np.linspace(0, width3-1, num=width3, endpoint=True) * nrlo/nrli + (nrlo-nrli)/(2.0*nrli)
            ionrect = np.zeros((length3, width3), dtype=np.float32)
            for i in range(length2):
                f = interp1d(index2, ionfilt[i,:], kind='cubic', fill_value="extrapolate")
                ionrect[i, :] = f(index3)
            
            index2 = np.linspace(0, length2-1, num=length2, endpoint=True)
            index3 = np.linspace(0, length3-1, num=length3, endpoint=True) * nalo/nali + (nalo-nali)/(2.0*nali)
            for j in range(width3):
                f = interp1d(index2, ionrect[0:length2, j], kind='cubic', fill_value="extrapolate")
                ionrect[:, j] = f(index3)
            ionrect.astype(np.float32).tofile(ionrectfile)
            del ionrect
        create_xml(ionrectfile, width3, length3, 'float')

        os.rename(ionrectfile, os.path.join('../../insar', ionrectfile))
        os.rename(ionrectfile+'.vrt', os.path.join('../../insar', ionrectfile)+'.vrt')
        os.rename(ionrectfile+'.xml', os.path.join('../../insar', ionrectfile)+'.xml')
        os.chdir('../../insar')
    else:
        shutil.copyfile(ionfiltfile, os.path.join('../../insar', ionrectfile))
        os.chdir('../../insar')
        create_xml(ionrectfile, width3, length3, 'float')
    #now we are in 'insar'


    ############################################################
    # STEP 4. correct interferogram
    ############################################################
    from isceobj.Alos2Proc.Alos2ProcPublic import renameFile
    from isceobj.Alos2Proc.Alos2ProcPublic import runCmd

    if self.applyIon:
        print('\ncorrect interferogram')
        if os.path.isfile(self._insar.multilookDifferentialInterferogramOriginal):
            print('original interferogram: {} is already here, do not rename: {}'.format(self._insar.multilookDifferentialInterferogramOriginal, self._insar.multilookDifferentialInterferogram))
        else:
            print('renaming {} to {}'.format(self._insar.multilookDifferentialInterferogram, self._insar.multilookDifferentialInterferogramOriginal))
            renameFile(self._insar.multilookDifferentialInterferogram, self._insar.multilookDifferentialInterferogramOriginal)

        cmd = "imageMath.py -e='a*exp(-1.0*J*b)' --a={} --b={} -s BIP -t cfloat -o {}".format(
            self._insar.multilookDifferentialInterferogramOriginal,
            self._insar.multilookIon,
            self._insar.multilookDifferentialInterferogram)
        runCmd(cmd)
    else:
        print('\nionospheric phase estimation finished, but correction of interfeorgram not requested')

    os.chdir('../')

    catalog.printToLog(logger, "runIonFilt")
    self._insar.procDoc.addAllFromCatalog(catalog)



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

    lowerUnw:        lower band unwrapped interferogram
    upperUnw:        upper band unwrapped interferogram
    wgt:             weight
    fl:              lower band center frequency
    fu:              upper band center frequency
    adjFlag:         method for removing relative phase unwrapping errors
                       0: mean value
                       1: polynomial
    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)*(wgt!=0)
        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, wgt, width, length, 1, 1, 1, 1, 2)

    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 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, weight, width, length, nrli, nali, nrlo, nalo, order):
    '''
    ionos:  input ionospheric phase
    weight: weight
    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
    '''

    from isceobj.Alos2Proc.Alos2ProcPublic import create_multi_index2

    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)

    #look for data to use
    flag = (weight!=0)*(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 = weight[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 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 reformatMaskedAreas(maskedAreas, length, width):
    '''
    reformat masked areas coordinates that are ready to use
    'maskedAreas' is a 2-D list. Each element in the 2-D list is a four-element list: [firstLine,
    lastLine, firstColumn, lastColumn], with line/column numbers starting with 1. If one of the
    four elements is specified with -1, the program will use firstLine/lastLine/firstColumn/
    lastColumn instead.

    output is a 2-D list containing the corresponding python-list/array-format indexes.
    '''
    numberOfAreas = len(maskedAreas)
    maskedAreasReformated = [[0, length, 0, width] for i in range(numberOfAreas)]

    for i in range(numberOfAreas):
        if maskedAreas[i][0] != -1:
            maskedAreasReformated[i][0] = maskedAreas[i][0] - 1
        if maskedAreas[i][1] != -1:
            maskedAreasReformated[i][1] = maskedAreas[i][1]
        if maskedAreas[i][2] != -1:
            maskedAreasReformated[i][2] = maskedAreas[i][2] - 1
        if maskedAreas[i][3] != -1:
            maskedAreasReformated[i][3] = maskedAreas[i][3]
        if (not (0 <= maskedAreasReformated[i][0] <= length-1)) or \
           (not (1 <= maskedAreasReformated[i][1] <= length)) or \
           (not (0 <= maskedAreasReformated[i][2] <= width-1)) or \
           (not (1 <= maskedAreasReformated[i][3] <= width)) or \
           (not (maskedAreasReformated[i][1]-maskedAreasReformated[i][0]>=1)) or \
           (not (maskedAreasReformated[i][3]-maskedAreasReformated[i][2]>=1)):
            raise Exception('area {} masked out in ionospheric phase estimation not correct'.format(i+1))

    return maskedAreasReformated