Merge pull request #236 from lijun99/pycuampcor

Updates to PyCuAmpcor

components/isceobj/Alos2Proc/runDenseOffset.py

@@ -240,15 +240,11 @@ def runDenseOffsetGPU(self):
     print('dense offset skip width: %d' % (self.offsetSkipWidth))
     print('dense offset skip hight: %d' % (self.offsetSkipHeight))
     print('dense offset covariance surface oversample factor: %d' % (self.offsetCovarianceOversamplingFactor))
-    print('dense offset covariance surface oversample window size: %d\n' % (self.offsetCovarianceOversamplingWindowsize))
 
 
     objOffset = PyCuAmpcor.PyCuAmpcor()
     objOffset.algorithm = 0
-    objOffset.deviceID = -1
+    objOffset.derampMethod = 1 # 1=linear phase ramp, 0=take mag, 2=skip
-    objOffset.nStreams = 2
-    #original ampcor program in roi_pac uses phase gradient to deramp
-    objOffset.derampMethod = 2
     objOffset.referenceImageName = self._insar.referenceSlc
     objOffset.referenceImageHeight = m.length
     objOffset.referenceImageWidth = m.width
@@ -256,56 +252,79 @@ def runDenseOffsetGPU(self):
     objOffset.secondaryImageHeight = s.length
     objOffset.secondaryImageWidth = s.width
     objOffset.offsetImageName = self._insar.denseOffset
+    objOffset.grossOffsetImageName = self._insar.denseOffset + ".gross"
     objOffset.snrImageName = self._insar.denseOffsetSnr
+    objOffset.covImageName = self._insar.denseOffsetCov
 
     objOffset.windowSizeWidth = self.offsetWindowWidth
     objOffset.windowSizeHeight = self.offsetWindowHeight
-    #objOffset.halfSearchRangeAcross = int(self.offsetSearchWindowWidth / 2 + 0.5)
+
-    #objOffset.halfSearchRangeDown = int(self.offsetSearchWindowHeight / 2 + 0.5)
     objOffset.halfSearchRangeAcross = self.offsetSearchWindowWidth
     objOffset.halfSearchRangeDown = self.offsetSearchWindowHeight
+
     objOffset.skipSampleDown = self.offsetSkipHeight
     objOffset.skipSampleAcross = self.offsetSkipWidth
+
     #Oversampling method for correlation surface(0=fft,1=sinc)
-    objOffset.corrSufaceOverSamplingMethod = 0
+    objOffset.corrSurfaceOverSamplingMethod = 0
     objOffset.corrSurfaceOverSamplingFactor = self.offsetCovarianceOversamplingFactor
-    objOffset.corrSurfaceZoomInWindow = self.offsetCovarianceOversamplingWindowsize
+
+    # set gross offset
     objOffset.grossOffsetAcrossStatic = 0
     objOffset.grossOffsetDownStatic = 0
+    # set the margin
+    margin = 0
 
-    objOffset.referenceStartPixelDownStatic = self.offsetWindowHeight//2
+    # adjust the margin
-    objOffset.referenceStartPixelAcrossStatic = self.offsetWindowWidth//2
+    margin = max(margin, abs(objOffset.grossOffsetAcrossStatic), abs(objOffset.grossOffsetDownStatic))
 
-    objOffset.numberWindowDown = (m.length - 2*self.offsetSearchWindowHeight - self.offsetWindowHeight) // self.offsetSkipHeight
+    # set the starting pixel of the first reference window
-    objOffset.numberWindowAcross = (m.width - 2*self.offsetSearchWindowWidth - self.offsetWindowWidth) // self.offsetSkipWidth
+    objOffset.referenceStartPixelDownStatic = margin + self.offsetSearchWindowHeight
+    objOffset.referenceStartPixelAcrossStatic = margin + self.offsetSearchWindowWidth
 
-    # generic control
+    # find out the total number of windows
-    objOffset.numberWindowDownInChunk = 8
+    objOffset.numberWindowDown = (m.length - 2*margin - 2*self.offsetSearchWindowHeight - self.offsetWindowHeight) // self.offsetSkipHeight
-    objOffset.numberWindowAcrossInChunk = 8
+    objOffset.numberWindowAcross = (m.width - 2*margin - 2*self.offsetSearchWindowWidth - self.offsetWindowWidth) // self.offsetSkipWidth
+
+    # gpu job control
+    objOffset.deviceID = 0
+    objOffset.nStreams = 2
+    objOffset.numberWindowDownInChunk = 1
+    objOffset.numberWindowAcrossInChunk = 64
     objOffset.mmapSize = 16
 
+    # pass/adjust the parameters
     objOffset.setupParams()
-    objOffset.setConstantGrossOffset(0, 0)
+    # set up the starting pixels for each window, based on the gross offset
+    objOffset.setConstantGrossOffset(objOffset.grossOffsetAcrossStatic, objOffset.grossOffsetDownStatic)
+    # check whether all pixels are in image range (optional)
     objOffset.checkPixelInImageRange()
+    print('\n======================================')
+    print('Running PyCuAmpcor...')
+    print('======================================\n')
     objOffset.runAmpcor()
 
     ### Store params for later
-    self._insar.offsetImageTopoffset = objOffset.halfSearchRangeDown
+    # location of the center of the first reference window
-    self._insar.offsetImageLeftoffset = objOffset.halfSearchRangeAcross
+    self._insar.offsetImageTopoffset = objOffset.referenceStartPixelDownStatic + (objOffset.windowSizeHeight-1)//2
-
+    self._insar.offsetImageLeftoffset = objOffset.referenceStartPixelAcrossStatic +(objOffset.windowSizeWidth-1)//2
 
+    # offset image dimension,  the number of windows
     width = objOffset.numberWindowAcross
     length = objOffset.numberWindowDown
-    offsetBIP = np.fromfile(objOffset.offsetImageName.decode('utf-8'), dtype=np.float32).reshape(length, width*2)
+
+    # convert the offset image from BIP to BIL
+    offsetBIP = np.fromfile(objOffset.offsetImageName, dtype=np.float32).reshape(length, width*2)
     offsetBIL = np.zeros((length*2, width), dtype=np.float32)
     offsetBIL[0:length*2:2, :] = offsetBIP[:, 1:width*2:2]
     offsetBIL[1:length*2:2, :] = offsetBIP[:, 0:width*2:2]
-    os.remove(objOffset.offsetImageName.decode('utf-8'))
+    os.remove(objOffset.offsetImageName)
-    offsetBIL.astype(np.float32).tofile(objOffset.offsetImageName.decode('utf-8'))
+    offsetBIL.astype(np.float32).tofile(objOffset.offsetImageName)
 
+    # generate offset image description files
     outImg = isceobj.createImage()
     outImg.setDataType('FLOAT')
-    outImg.setFilename(objOffset.offsetImageName.decode('utf-8'))
+    outImg.setFilename(objOffset.offsetImageName)
     outImg.setBands(2)
     outImg.scheme = 'BIL'
     outImg.setWidth(objOffset.numberWindowAcross)
@@ -314,8 +333,11 @@ def runDenseOffsetGPU(self):
     outImg.setAccessMode('read')
     outImg.renderHdr()
 
+    # gross offset image is not needed, since all zeros
+
+    # generate snr image description files
     snrImg = isceobj.createImage()
-    snrImg.setFilename( objOffset.snrImageName.decode('utf8'))
+    snrImg.setFilename( objOffset.snrImageName)
     snrImg.setDataType('FLOAT')
     snrImg.setBands(1)
     snrImg.setWidth(objOffset.numberWindowAcross)
@@ -323,4 +345,20 @@ def runDenseOffsetGPU(self):
     snrImg.setAccessMode('read')
     snrImg.renderHdr()
 
+    # generate cov image description files
+    # covariance of azimuth/range offsets.
+    # 1st band: cov(az, az), 2nd band: cov(rg, rg), 3rd band: cov(az, rg)
+    covImg = isceobj.createImage()
+    covImg.setFilename(objOffset.covImageName)
+    covImg.setDataType('FLOAT')
+    covImg.setBands(3)
+    covImg.scheme = 'BIP'
+    covImg.setWidth(objOffset.numberWindowAcross)
+    covImg.setLength(objOffset.numberWindowDown)
+    outImg.addDescription('covariance of azimuth/range offsets')
+    covImg.setAccessMode('read')
+    covImg.renderHdr()
+
     return (objOffset.numberWindowAcross, objOffset.numberWindowDown)
+
+# end of file

components/isceobj/TopsProc/runDenseOffsets.py

@@ -156,11 +156,11 @@ def runDenseOffsetsGPU(self):
     secondary = os.path.join(self._insar.mergedDirname, sf)
 
     ####For this module currently, we need to create an actual file on disk
-
     for infile in [reference,secondary]:
         if os.path.isfile(infile):
             continue
 
+        print('Creating actual file {} ...\n'.format(infile))
         cmd = 'gdal_translate -of ENVI {0}.vrt {0}'.format(infile)
         status = os.system(cmd)
         if status:
@@ -170,21 +170,27 @@ def runDenseOffsetsGPU(self):
     m = isceobj.createSlcImage()
     m.load(reference + '.xml')
     m.setAccessMode('READ')
-#    m.createImage()
+    # re-create vrt in terms of merged full slc
+    m.renderHdr()
 
     ### Load the secondary object
     s = isceobj.createSlcImage()
     s.load(secondary + '.xml')
     s.setAccessMode('READ')
-#    s.createImage()
+    # re-create vrt in terms of merged full slc
+    s.renderHdr()
 
+    # get the dimension
     width = m.getWidth()
     length = m.getLength()
 
+    ### create the GPU processor
     objOffset = PyCuAmpcor.PyCuAmpcor()
+
+    ### Set parameters
+    # cross-correlation method, 0=Frequency domain, 1= Time domain
     objOffset.algorithm = 0
-    objOffset.deviceID = -1
+    # deramping method: 0 to take magnitude (fixed for Tops)
-    objOffset.nStreams = 2
     objOffset.derampMethod = 0
     objOffset.referenceImageName = reference + '.vrt'
     objOffset.referenceImageHeight = length
@@ -193,36 +199,59 @@ def runDenseOffsetsGPU(self):
     objOffset.secondaryImageHeight = length
     objOffset.secondaryImageWidth = width
 
-    objOffset.numberWindowDown = (length-100-self.winhgt)//self.skiphgt
+    # adjust the margin
-    objOffset.numberWindowAcross = (width-100-self.winwidth)//self.skipwidth
+    margin = max(self.margin, abs(self.azshift), abs(self.rgshift))
 
+    # set the start pixel in the reference image
+    objOffset.referenceStartPixelDownStatic = margin + self.srchgt
+    objOffset.referenceStartPixelAcrossStatic = margin + self.srcwidth
+
+    # compute the number of windows
+    objOffset.numberWindowDown = (length-2*margin-2*self.srchgt-self.winhgt)//self.skiphgt
+    objOffset.numberWindowAcross = (width-2*margin-2*self.srcwidth-self.winwidth)//self.skipwidth
+
+    # set the template window size
     objOffset.windowSizeHeight = self.winhgt
     objOffset.windowSizeWidth = self.winwidth
 
+    # set the (half) search range
     objOffset.halfSearchRangeDown = self.srchgt
     objOffset.halfSearchRangeAcross = self.srcwidth
 
-    objOffset.referenceStartPixelDownStatic = 50
+    # set the skip distance between windows
-    objOffset.referenceStartPixelAcrossStatic = 50
-
     objOffset.skipSampleDown = self.skiphgt
     objOffset.skipSampleAcross = self.skipwidth
 
-    objOffset.corrSufaceOverSamplingMethod = 0
+    # correlation surface oversampling method, # 0=FFT, 1=Sinc
+    objOffset.corrSurfaceOverSamplingMethod = 0
+    # oversampling factor
     objOffset.corrSurfaceOverSamplingFactor = self.oversample
 
-    # generic control
+    ### gpu control
-    objOffset.numberWindowDownInChunk = 10
+    objOffset.deviceID = 0
-    objOffset.numberWindowAcrossInChunk = 10
+    objOffset.nStreams = 2
+    # number of windows in a chunk/batch
+    objOffset.numberWindowDownInChunk = 1
+    objOffset.numberWindowAcrossInChunk = 64
+    # memory map cache size in GB
     objOffset.mmapSize = 16
 
+    # Modify BIL in filename to BIP if needed and store for future use
+    prefix, ext = os.path.splitext(self._insar.offsetfile)
+    if ext == '.bil':
+        ext = '.bip'
+        self._insar.offsetfile = prefix + ext
 
-    objOffset.setupParams()
+    # set the output file name
+    objOffset.offsetImageName = os.path.join(self._insar.mergedDirname, self._insar.offsetfile)
+    objOffset.grossOffsetImageName = os.path.join(self._insar.mergedDirname, self._insar.offsetfile + ".gross")
+    objOffset.snrImageName = os.path.join(self._insar.mergedDirname, self._insar.snrfile)
+    objOffset.covImageName = os.path.join(self._insar.mergedDirname, self._insar.covfile)
 
-    objOffset.setConstantGrossOffset(self.azshift,self.rgshift)
+    # merge gross offset to final offset
+    objOffset.mergeGrossOffset = 1
 
-#    objOffset.numberThreads = 1
+    ### print the settings
-    ### Configure dense Ampcor object
     print('\nReference frame: %s' % (mf))
     print('Secondary frame: %s' % (sf))
     print('Main window size width: %d' % (self.winwidth))
@@ -231,41 +260,41 @@ def runDenseOffsetsGPU(self):
     print('Search window size height: %d' % (self.srchgt))
     print('Skip sample across: %d' % (self.skipwidth))
     print('Skip sample down: %d' % (self.skiphgt))
-    print('Field margin: %d' % (self.margin))
+    print('Field margin: %d' % (margin))
     print('Oversampling factor: %d' % (self.oversample))
     print('Gross offset across: %d' % (self.rgshift))
     print('Gross offset down: %d\n' % (self.azshift))
-
-    #Modify BIL in filename to BIP if needed and store for future use
-    prefix, ext = os.path.splitext(self._insar.offsetfile)
-    if ext == '.bil':
-        ext = '.bip'
-        self._insar.offsetfile = prefix + ext
-
-    objOffset.offsetImageName = os.path.join(self._insar.mergedDirname, self._insar.offsetfile)
-    objOffset.snrImageName = os.path.join(self._insar.mergedDirname, self._insar.snrfile)
-
     print('Output dense offsets file name: %s' % (objOffset.offsetImageName))
+    print('Output gross offsets file name: %s' % (objOffset.grossOffsetImageName))
     print('Output SNR file name: %s' % (objOffset.snrImageName))
+    print('Output COV file name: %s' % (objOffset.covImageName))
+
+    # pass the parameters to C++ programs
+    objOffset.setupParams()
+    # set the (static) gross offset
+    objOffset.setConstantGrossOffset(self.azshift,self.rgshift)
+    # make sure all pixels are in range
+    objOffset.checkPixelInImageRange()
+
     print('\n======================================')
-    print('Running dense ampcor...')
+    print('Running PyCuAmpcor...')
     print('======================================\n')
 
-
+    # run ampcor
-    objOffset.checkPixelInImageRange()
     objOffset.runAmpcor()
 
-    #objOffset.denseampcor(m, s) ### Where the magic happens...
-
     ### Store params for later
+    # offset width x length, also number of windows
     self._insar.offset_width = objOffset.numberWindowAcross
     self._insar.offset_length = objOffset.numberWindowDown
-    self._insar.offset_top = 50
+    # the center of the first reference window
-    self._insar.offset_left = 50
+    self._insar.offset_top = objOffset.referenceStartPixelDownStatic + (objOffset.windowSizeHeight-1)//2
+    self._insar.offset_left = objOffset.referenceStartPixelAcrossStatic + (objOffset.windowSizeWidth-1)//2
 
+    # generate description files for output images
     outImg = isceobj.createImage()
     outImg.setDataType('FLOAT')
-    outImg.setFilename(objOffset.offsetImageName.decode('utf-8'))
+    outImg.setFilename(objOffset.offsetImageName)
     outImg.setBands(2)
     outImg.scheme = 'BIP'
     outImg.setWidth(objOffset.numberWindowAcross)
@@ -273,8 +302,19 @@ def runDenseOffsetsGPU(self):
     outImg.setAccessMode('read')
     outImg.renderHdr()
 
+    # gross offset
+    goutImg = isceobj.createImage()
+    goutImg.setDataType('FLOAT')
+    goutImg.setFilename(objOffset.grossOffsetImageName)
+    goutImg.setBands(2)
+    goutImg.scheme = 'BIP'
+    goutImg.setWidth(objOffset.numberWindowAcross)
+    goutImg.setLength(objOffset.numberWindowDown)
+    goutImg.setAccessMode('read')
+    goutImg.renderHdr()
+
     snrImg = isceobj.createImage()
-    snrImg.setFilename( objOffset.snrImageName.decode('utf-8'))
+    snrImg.setFilename(objOffset.snrImageName)
     snrImg.setDataType('FLOAT')
     snrImg.setBands(1)
     snrImg.setWidth(objOffset.numberWindowAcross)
@@ -282,6 +322,15 @@ def runDenseOffsetsGPU(self):
     snrImg.setAccessMode('read')
     snrImg.renderHdr()
 
+    covImg = isceobj.createImage()
+    covImg.setFilename(objOffset.covImageName)
+    covImg.setDataType('FLOAT')
+    covImg.setBands(3)
+    covImg.scheme = 'BIP'
+    covImg.setWidth(objOffset.numberWindowAcross)
+    covImg.setLength(objOffset.numberWindowDown)
+    covImg.setAccessMode('read')
+    covImg.renderHdr()
 
 
 if __name__ == '__main__' :

contrib/PyCuAmpcor/README.md

@@ -91,7 +91,7 @@ mm=0   # margin to be neglected
 gross=0  # whether to use a varying gross offset
 azshift=0 # constant gross offset along height/azimuth 
 rgshift=0 # constant gross offset along width/range
-deramp=0 # 0 for mag (TOPS), 1 for complex  
+deramp=0 # 0 for mag (TOPS), 1 for complex linear ramp, 2 for complex no deramping  
 oo=32  # correlation surface oversampling factor 
 outprefix=./merged/20151120_20151214/offset  # output prefix
 outsuffix=_ww64_wh64   # output suffix
@@ -368,8 +368,7 @@ RIOPAC extracts the secondary window centering at the correlation surface peak.
 | PyCuAmpcor          | CUDA variable       | ampcor.F equivalent   | Notes                     |
 | :---                | :---                | :----                 | :---                      |
 | rawDataOversamplingFactor | rawDataOversamplingFactor | i_ovs=2   | the oversampling factor for reference and secondary windows, use 2 for InSAR SLCs. |
-| derampMethod        | derampMethod        | 1 or no effect on TOPS | 0=mag for TOPS, 1=deramping (default), else=skip deramping.
+| derampMethod        | derampMethod        | 1 or no effect on TOPS | Only for complex: 0=take mag (TOPS), 1=linear deramp (default), else=skip deramp.
-
 
 **Difference to ROIPAC**
 

contrib/PyCuAmpcor/examples/cuDenseOffsets.py

@@ -89,6 +89,8 @@ def createParser():
                        help='Gross range offset (default: %(default)s).')
     gross.add_argument('--gf', '--gross-file', type=str, dest='gross_offset_file',
                        help='Varying gross offset input file')
+    gross.add_argument('--mg', '--merge-gross-offset', type=int, dest='merge_gross_offset', default=0,
+                       help='Whether to merge gross offset to the output offset image (default: %(default)s).')
 
     corr = parser.add_argument_group('Correlation surface')
     corr.add_argument('--corr-stat-size', type=int, dest='corr_stat_win_size', default=21,
@@ -104,7 +106,7 @@ def createParser():
     # gpu settings
     proc = parser.add_argument_group('Processing parameters')
     proc.add_argument('--gpuid', '--gid', '--gpu-id', dest='gpuid', type=int, default=0,
-                        help='GPU ID (default: %(default)).')
+                        help='GPU ID (default: %(default)s).')
     proc.add_argument('--nstreams', dest='nstreams', type=int, default=2,
                         help='Number of cuda streams (default: %(default)s).')
     proc.add_argument('--usemmap', dest='usemmap', type=int, default=1,
@@ -238,6 +240,9 @@ def estimateOffsetField(reference, secondary, inps=None):
     print("snr: ",objOffset.snrImageName)
     print("cov: ",objOffset.covImageName)
 
+    # whether to include the gross offset in offsetImage
+    objOffset.mergeGrossOffset = inps.merge_gross_offset
+
     offsetImageName = objOffset.offsetImageName
     grossOffsetImageName = objOffset.grossOffsetImageName
     snrImageName = objOffset.snrImageName

contrib/PyCuAmpcor/examples/grossOffsets.py

@@ -0,0 +1,407 @@
+#!/usr/bin/env python3
+# Generate pixel offsets based on Antarctica velocity model (MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2 doi:https://doi.org/10.5067/D7GK8F5J8M8R)
+# Author: Minyan Zhong
+import os
+import argparse
+import isce
+import isceobj
+import gdal
+import pyproj
+import numpy as np
+import matplotlib.pyplot as plt
+
+EXAMPLE = '''
+grossOffsets.py --model_file antarctica_ice_velocity_450m_v2.nc --lon lon.rdr --lat lat.rdr --los los.rdr --los_scheme bil --ww 64 --wh 64 --sw 10 --sh 10 --mm 50 --kw 32 --kh 32 --startpixeldw 50 --startpixelac 50 --rangePixelSize 0.930 --azimuthPixelSize 2.286 --interval 1
+'''
+
+def createParser():
+    '''
+    Command line parser.
+    '''
+    
+    parser = argparse.ArgumentParser(description='Generate pixel offsets (integer pixel) based on Antarctica ice velocity model (MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2 doi:https://doi.org/10.5067/D7GK8F5J8M8R)', formatter_class=argparse.RawTextHelpFormatter, epilog=EXAMPLE)
+
+    # path to antarctica velocity model
+    parser.add_argument('--model_file', type=str, dest='model_file', required=True)
+
+    # lat, lon, los
+    parser.add_argument('--lat', type=str, dest='lat', required=True,
+                        help='latitude file')
+    parser.add_argument('--lon', type=str, dest='lon', required=True,
+                        help='longitude fie')
+
+    parser.add_argument('--los', type=str, dest='los', required=True,
+                        help='two bands raster data in float. band1: incidence angle; bands: satellite flight direction (ISCE2 convention)')
+
+    parser.add_argument('--los_scheme', type=str, dest='los_scheme', required=True,
+                        help='interleave scheme of los (bil, bsq or bip)')
+
+    # window size settings
+    parser.add_argument('--ww', type=int, dest='winwidth', default=64,
+                        help='Window width (default: %(default)s).')
+    parser.add_argument('--wh', type=int, dest='winhgt', default=64,
+                        help='Window height (default: %(default)s).')
+    parser.add_argument('--sw', type=int, dest='srcwidth', default=20,
+                        help='Half search range along width, (default: %(default)s, recommend: 4-32).')
+    parser.add_argument('--sh', type=int, dest='srchgt', default=20,
+                        help='Half search range along height (default: %(default)s, recommend: 4-32).')
+    parser.add_argument('--kw', type=int, dest='skipwidth', default=64,
+                        help='Skip across (default: %(default)s).')
+    parser.add_argument('--kh', type=int, dest='skiphgt', default=64,
+                        help='Skip down (default: %(default)s).')
+
+    # determine the number of windows
+    # either specify the starting pixel and the number of windows,
+    # or by setting them to -1, let the script to compute these parameters
+    parser.add_argument('--mm', type=int, dest='margin', default=0,
+                        help='Margin (default: %(default)s).')
+    
+    parser.add_argument('--spa','--startpixelac', dest='startpixelac', type=int, default=-1, help='Starting Pixel across of the reference image(default: %(default)s to be determined by margin and search range).')
+    
+    parser.add_argument('--spd','--startpixeldw', dest='startpixeldw', type=int, default=-1, help='Starting Pixel down of the reference image (default: %(default)s).')
+    
+    parser.add_argument('--aps', '--azimuthPixelSize', dest='azimuthPixelSize', type=float, required=True, help='azimuth pixel size')
+    
+    parser.add_argument('--rps', '--rangePixelSize', dest='rangePixelSize', type=float, required=True, help='range pixel size')
+
+    parser.add_argument('--interval', dest='interval', type=float, required=True, help='interval between reference and secondary scene (unit: day)')
+
+    parser.add_argument('--outdir', dest='outdir', type=str, default='.', help='output directory')
+
+    parser.add_argument('--outname', dest='outname', type=str, default='grossOffsets.bin', help='output name of gross pixel offsets (integer)')
+
+    return parser
+
+def cmdLineParse(iargs = None):
+    parser = createParser()
+    inps = parser.parse_args(args=iargs)
+    return inps
+
+class grossOffsets:
+    def __init__(self, inps):
+        model_path = inps.model_file
+        self.model_file = model_path
+        self.latfile = inps.lat
+        self.lonfile = inps.lon
+        self.losfile = inps.los
+
+        ds = gdal.Open(self.losfile)
+        self.XSize = ds.RasterXSize
+        self.YSize = ds.RasterYSize
+        ds = None
+        
+        self.los_scheme = inps.los_scheme.lower()
+        assert(self.los_scheme in ['bil','bsq', 'bip']), print('interleave scheme of los')
+
+        self.margin = inps.margin
+        self.winSizeHgt = inps.winhgt
+        self.winSizeWidth = inps.winwidth
+        self.searchSizeHgt = inps.srchgt
+        self.searchSizeWidth = inps.srcwidth
+        self.skipSizeHgt = inps.skiphgt
+        self.skipSizeWidth = inps.skipwidth
+
+        self.startpixelac = inps.startpixelac if inps.startpixelac != -1 else self.margin + self.searchSizeWidth
+ 
+        self.startpixeldw = inps.startpixeldw if inps.startpixeldw != -1 else self.margin + self.searchSizeHgt
+ 
+        self.azPixelSize = inps.azimuthPixelSize
+        self.rngPixelSize = inps.rangePixelSize
+
+        self.interval = inps.interval
+
+        self.outdir = inps.outdir
+        self.outname = inps.outname
+
+        self.get_veloData()
+        self.vProj = pyproj.Proj('+init=EPSG:3031')
+
+    def get_veloData(self):
+        assert os.path.exists(self.model_file), print("Please download MEaSUREs InSAR-Based Antarctica Ice Velocity Map, Version 2 at https://nsidc.org/data/NSIDC-0484/versions")
+
+        data_read = 0
+        ds = gdal.Open("NETCDF:{0}:{1}".format(self.model_file, 'VX'))
+        self.vx = ds.ReadAsArray()
+
+        ds = gdal.Open("NETCDF:{0}:{1}".format(self.model_file, 'VY'))
+        self.vy = ds.ReadAsArray()
+
+        self.vx = np.flipud(self.vx)
+        self.vy = np.flipud(self.vy) 
+
+        self.v = np.sqrt(np.multiply(self.vx,self.vx)+np.multiply(self.vy,self.vy))
+
+        self.model_spacing = 450
+        self.x0 = np.arange(-2800000,2800000,step=450)
+        self.y0 = np.arange(-2800000,2800000,step=450)+200
+
+    def runGrossOffsets(self):
+        ## Step 0: Set up projection transformers for ease of use
+        self.llhProj = pyproj.Proj('+init=EPSG:4326')
+        self.xyzProj = pyproj.Proj('+init=EPSG:4978')
+
+        # From xy to lat lon.
+        refPt = self.vProj(0.0, 0.0, inverse=True)
+       
+        ### Step 2: Cut the data
+        print('Extract the data to this radar scene...')
+        # The following code is to be consistent with "get_offset_geometry" in dense_offset.py
+
+        numWinDown = (self.YSize - self.margin*2 - self.searchSizeHgt*2 - self.winSizeHgt) // self.skipSizeHgt
+        numWinAcross = (self.XSize - self.margin*2 - self.searchSizeWidth*2 - self.winSizeWidth) // self.skipSizeWidth
+        
+        lat = np.zeros(shape=(numWinDown,numWinAcross),dtype=np.float64) 
+        lon = np.zeros(shape=(numWinDown,numWinAcross),dtype=np.float64)
+        inc = np.zeros(shape=(numWinDown,numWinAcross),dtype=np.float32)
+        azi = np.zeros(shape=(numWinDown,numWinAcross),dtype=np.float32)
+
+        self.centerOffsetHgt = self.winSizeHgt//2-1
+        self.centerOffsetWidth = self.winSizeWidth//2-1
+
+        print("Number of winows in down direction, Number of window in across direction: ")
+        print(numWinDown, numWinAcross)
+
+        cut_vx = np.zeros(shape=(numWinDown,numWinAcross))
+        cut_vy = np.zeros(shape=(numWinDown,numWinAcross))
+        cut_v = np.zeros(shape=(numWinDown,numWinAcross))
+        pixel = np.zeros(shape=(numWinDown,numWinAcross))
+        line = np.zeros(shape=(numWinDown,numWinAcross))
+
+        for iwin in range(numWinDown):
+            # Need to calculate lat lon in the interior mode.
+            print('Processing line: ',iwin, 'out of', numWinDown)
+            down = self.margin + self.skipSizeHgt * iwin  + self.centerOffsetHgt
+            off = down*self.XSize
+
+            across_indices = self.margin + np.arange(numWinAcross)*self.skipSizeWidth + self.centerOffsetWidth
+
+            # latitude
+            latline = np.memmap(filename=self.latfile,dtype='float64',offset=8*off,shape=(self.XSize))
+            # longitude
+            lonline = np.memmap(filename=self.lonfile,dtype='float64',offset=8*off,shape=(self.XSize))
+
+            # incidence angle and satellite flight direction
+            # bil
+            if self.los_scheme == "bil":
+                off2 = down * self.XSize * 2
+                losline = np.memmap(filename=self.losfile,dtype='float32',offset=4*off2,shape=(self.XSize*2))
+
+                incline = losline[0:self.XSize]
+                aziline = losline[self.XSize:self.XSize*2]
+            # bsq
+            elif self.los_scheme == 'bsq':
+                off2 = self.YSize * self.XSize + down * self.XSize
+                incline = np.memmap(filename=self.losfile,dtype='float32',offset=4*off,shape=(self.XSize))
+                aziline = np.memmap(filename=self.losfile,dtype='float32',offset=4*off2,shape=(self.XSize))
+            # bip 
+            else:
+               off2 = down * self.XSize * 2
+               losline = np.memmap(filename=self.losfile,dtype='float32',offset=4*off2,shape=(self.XSize*2))
+               incline = losline[0:self.XSize*2:2]
+               aziline = losline[1:self.XSize*2:2]
+
+            # Subset the line
+            lat[iwin,:] = latline[across_indices]
+            lon[iwin,:] = lonline[across_indices]
+            inc[iwin,:] = incline[across_indices]
+            azi[iwin,:] = aziline[across_indices]
+
+            #print(iwin,'lat: ',lat[iwin,:])
+            #print(iwin,'lon: ',lon[iwin,:])
+            #print(iwin,'inc: ',inc[iwin,:])
+            #print(iwin,'azi: ',azi[iwin,:])
+
+            #### Look up in MEaSUREs InSAR-Based Antarctica Ice Velocity Map
+
+            # Convert lat lon to grid coordinates in polar stereographic projection.
+            xyMap = pyproj.transform(self.llhProj, self.vProj, lon[iwin,:], lat[iwin,:])
+ 
+            # Extract the values in the velocity model.
+            model_spacing = self.model_spacing
+            pixel[iwin,:] = np.clip((xyMap[0]-self.x0[0])/model_spacing, 0, self.vx.shape[1]-1)
+            line[iwin,:] = np.clip((xyMap[1]-self.y0[0])/model_spacing, 0, self.vx.shape[0]-1)
+
+            pixel_int = pixel[iwin,:].astype(int)
+            line_int = line[iwin,:].astype(int)
+
+            cut_vx[iwin,:] = self.vx[line_int,pixel_int]
+            cut_vy[iwin,:] = self.vy[line_int,pixel_int]
+
+        cut_v = np.sqrt(np.multiply(cut_vx,cut_vx),np.multiply(cut_vy,cut_vy))
+        valid = np.logical_and(inc!=0, cut_v!=0)
+
+        ### Mask out invalid values ###
+        # 1. Mask out invalid values at margin.
+        cut_vx[inc==0] = np.nan
+        cut_vy[inc==0] = np.nan
+ 
+        # Get Interpolated speed.
+        cut_v = np.sqrt(np.multiply(cut_vx,cut_vx),np.multiply(cut_vy,cut_vy))
+
+        print("The speed matrix") 
+        print(cut_v)
+        print("The shape of speed matrix")
+        print(cut_v.shape)
+
+        ### Step 3: Convert XY velocity to EN velocity (clockwise rotation)
+        print('Coverting XY to EN...')
+
+        lonr = np.radians(lon - refPt[0])
+        cut_ve = np.multiply(cut_vx, np.cos(lonr)) - np.multiply(cut_vy, np.sin(lonr))
+        cut_vn = np.multiply(cut_vy, np.cos(lonr)) + np.multiply(cut_vx, np.sin(lonr))
+
+        print('Polar stereographic velocity: ', [cut_vx, cut_vy])
+        print('Local ENU velocity: ', [cut_ve, cut_vn])
+ 
+        ####Step 4: Convert EN velocity to rng and azimuth
+        #Local los and azi vector in ENU coordinate
+        print(' Coverting EN to rdr...')
+        incr = np.radians(inc)
+        azir = np.radians(azi)
+        losr = np.radians(azi-90.0)
+
+        losenu=[ np.multiply(np.sin(incr),np.cos(losr)),
+                 np.multiply(np.sin(incr),np.sin(losr)),
+                 -np.cos(incr) ]
+        
+        azienu=[ np.cos(azir),
+                 np.sin(azir),
+                 0.0 ]
+
+        # unit: pixel per day
+        grossRangeOffset = (self.interval/365.25) * (cut_ve * losenu[0] + cut_vn * losenu[1])/ self.rngPixelSize
+        grossAzimuthOffset = (self.interval/365.25) * (cut_ve * azienu[0] + cut_vn * azienu[1]) / self.azPixelSize
+
+        # Mask out invalid values at margin.
+        grossRangeOffset[inc==0] = np.nan
+        grossAzimuthOffset[inc==0] = np.nan
+ 
+        print('Gross azimuth offset: ', grossAzimuthOffset)
+        print('Gross range offset: ', grossRangeOffset)
+        print('Shape of gross offsets: ', grossRangeOffset.shape)
+
+        ### Show FLOAT results ###
+        fig=plt.figure(21,figsize=(9,9))
+        ax = fig.add_subplot(121)
+        ax.set_title('gross azimuth offset',fontsize=15)
+        cax = ax.imshow(grossAzimuthOffset,cmap=plt.cm.coolwarm)
+        cbar = fig.colorbar(cax,shrink=0.8)
+        cbar.set_label("pixel",fontsize=15)
+
+        ax = fig.add_subplot(122)
+        ax.set_title('gross range offset',fontsize=15)
+        cax = ax.imshow(grossRangeOffset,cmap=plt.cm.coolwarm)
+        cbar = fig.colorbar(cax,shrink=0.8)
+        cbar.set_label("pixel",fontsize=15)
+
+        figname = os.path.join(self.outdir,'pixel_offsets.png')
+        fig.savefig(figname,format='png')
+        plt.close()
+
+        # Save grossRangeOffset and grossAzimuthOffset as ISCE supported images.
+        # Range
+        rangeFileName = os.path.join(self.outdir, 'grossRange.off')
+        driver = gdal.GetDriverByName('ENVI')
+        dst_ds = driver.Create(rangeFileName, xsize=grossRangeOffset.shape[1], ysize=grossRangeOffset.shape[0], bands=1, eType=gdal.GDT_Float32)
+        dst_ds.GetRasterBand(1).WriteArray(grossRangeOffset,0,0)
+        dst_ds = None
+
+        outImage = isceobj.createImage()
+        outImage.setDataType('FLOAT')
+        outImage.setFilename(rangeFileName)
+        outImage.setBands(1)
+        outImage.scheme='BIL'
+        outImage.setLength(grossRangeOffset.shape[0])
+        outImage.setWidth(grossRangeOffset.shape[1])
+        outImage.setAccessMode('read')
+        outImage.renderHdr()
+
+        # Azimuth
+        azimuthFileName = os.path.join(self.outdir, 'grossAzimuth.off')
+        driver = gdal.GetDriverByName('ENVI')
+        dst_ds = driver.Create(azimuthFileName, xsize=grossAzimuthOffset.shape[1], ysize=grossAzimuthOffset.shape[0], bands=1, eType=gdal.GDT_Float32)
+        dst_ds.GetRasterBand(1).WriteArray(grossAzimuthOffset,0,0)
+        dst_ds = None
+
+        outImage = isceobj.createImage()
+        outImage.setDataType('FLOAT')
+        outImage.setFilename(azimuthFileName)
+        outImage.setBands(1)
+        outImage.scheme='BIL'
+        outImage.setLength(grossAzimuthOffset.shape[0])
+        outImage.setWidth(grossAzimuthOffset.shape[1])
+        outImage.setAccessMode('read')
+        outImage.renderHdr()
+
+        ### Round to integer ###
+        grossAzimuthOffset_int = np.rint(grossAzimuthOffset).astype(np.int32)
+        grossRangeOffset_int = np.rint(grossRangeOffset).astype(np.int32)
+
+        ### Show Integer results ###
+        fig=plt.figure(22,figsize=(9,9))
+        ax = fig.add_subplot(121)
+        ax.set_title('gross azimuth offset (int)',fontsize=15)
+        cax = ax.imshow(grossAzimuthOffset_int,cmap=plt.cm.coolwarm)
+        cbar = fig.colorbar(cax,shrink=0.8)
+        cbar.set_label("pixel",fontsize=15)
+
+        ax = fig.add_subplot(122)
+        ax.set_title('gross range offset (int)',fontsize=15)
+        cax = ax.imshow(grossRangeOffset_int,cmap=plt.cm.coolwarm)
+        cbar = fig.colorbar(cax,shrink=0.8)
+        cbar.set_label("pixel",fontsize=15)
+
+        figname = os.path.join(self.outdir,'pixel_offsets_int.png')
+        fig.savefig(figname,format='png')
+        plt.close()
+
+        # Save grossRangeOffset and grossAzimuthOffset as ISCE supported images.
+        # Range
+        rangeFileName = os.path.join(self.outdir, 'grossRange_int.off')
+        driver = gdal.GetDriverByName('ENVI')
+        dst_ds = driver.Create(rangeFileName, xsize=grossRangeOffset.shape[1], ysize=grossRangeOffset.shape[0], bands=1, eType=gdal.GDT_Int32)
+        dst_ds.GetRasterBand(1).WriteArray(grossRangeOffset_int,0,0)
+        dst_ds = None
+
+        outImage = isceobj.createImage()
+        outImage.setDataType('INT')
+        outImage.setFilename(rangeFileName)
+        outImage.setBands(1)
+        outImage.scheme='BIL'
+        outImage.setLength(grossRangeOffset.shape[0])
+        outImage.setWidth(grossRangeOffset.shape[1])
+        outImage.setAccessMode('read')
+        outImage.renderHdr()
+
+        # Azimuth
+        azimuthFileName = os.path.join(self.outdir, 'grossAzimuth_int.off')
+        driver = gdal.GetDriverByName('ENVI')
+        dst_ds = driver.Create(azimuthFileName, xsize=grossAzimuthOffset.shape[1], ysize=grossAzimuthOffset.shape[0], bands=1, eType=gdal.GDT_Int32)
+        dst_ds.GetRasterBand(1).WriteArray(grossAzimuthOffset_int,0,0)
+        dst_ds = None
+
+        outImage = isceobj.createImage()
+        outImage.setDataType('INT')
+        outImage.setFilename(azimuthFileName)
+        outImage.setBands(1)
+        outImage.scheme='BIL'
+        outImage.setLength(grossAzimuthOffset.shape[0])
+        outImage.setWidth(grossAzimuthOffset.shape[1])
+        outImage.setAccessMode('read')
+        outImage.renderHdr()
+
+        # Round to integer and write to raw binary file
+        numTotal = numWinDown * numWinAcross
+        grossOffsets_int = np.hstack((grossAzimuthOffset_int.reshape(numTotal,1), grossRangeOffset_int.reshape(numTotal,1)))
+        print("grossOffsets: \n", grossOffsets_int, grossOffsets_int.dtype)
+        grossOffsets_int.tofile(os.path.join(self.outdir, self.outname))
+
+        return 0
+
+def main(iargs=None):
+    inps = cmdLineParse(iargs)
+    grossObj = grossOffsets(inps)
+    grossObj.runGrossOffsets()
+
+if __name__=='__main__':
+    main()

contrib/PyCuAmpcor/src/PyCuAmpcor.pyx

@@ -87,6 +87,8 @@ cdef extern from "cuAmpcorParameter.h":
         int *grossOffsetAcross              ## Gross offsets between reference and secondary windows (across)
         int grossOffsetDown0                ## constant gross offset (down)
         int grossOffsetAcross0              ## constant gross offset (across)
+        int mergeGrossOffset                ## whether to merge gross offset to final offset
+
         int referenceStartPixelDown0        ## the first pixel of reference image (down), be adjusted with margins and gross offset
         int referenceStartPixelAcross0      ## the first pixel of reference image (across)
         int *referenceChunkStartPixelDown   ## array of starting pixels for all reference chunks (down)
@@ -259,12 +261,7 @@ cdef class PyCuAmpcor(object):
     @secondaryImageName.setter
     def secondaryImageName(self, str a):
         self.c_cuAmpcor.param.secondaryImageName = <string> a.encode()
-    @property
+
-    def referenceImageName(self):
-        return self.c_cuAmpcor.param.referenceImageName
-    @referenceImageName.setter
-    def referenceImageName(self, str a):
-        self.c_cuAmpcor.param.referenceImageName = <string> a.encode()
     @property
     def referenceImageHeight(self):
         return self.c_cuAmpcor.param.referenceImageHeight
@@ -342,6 +339,13 @@ cdef class PyCuAmpcor(object):
     def offsetImageName(self, str a):
         self.c_cuAmpcor.param.offsetImageName = <string> a.encode()
 
+    @property
+    def mergeGrossOffset(self):
+        return self.c_cuAmpcor.param.mergeGrossOffset
+    @mergeGrossOffset.setter
+    def mergeGrossOffset(self, int a):
+        self.c_cuAmpcor.param.mergeGrossOffset = a
+
     @property
     def snrImageName(self):
         return self.c_cuAmpcor.param.snrImageName.decode("utf-8")

contrib/PyCuAmpcor/src/cuAmpcorChunk.cu

@@ -3,7 +3,7 @@
 
 /**
  * Run ampcor process for a batch of images (a chunk)
- * @param[in] idxDown_  index oIDIVUP(i,j) ((i+j-1)/j)f the chunk along Down/Azimuth direction
+ * @param[in] idxDown_  index of the chunk along Down/Azimuth direction
  * @param[in] idxAcross_ index of the chunk along Across/Range direction
  */
 void cuAmpcorChunk::run(int idxDown_, int idxAcross_)
@@ -65,29 +65,31 @@ void cuAmpcorChunk::run(int idxDown_, int idxAcross_)
     // 41 x 41, if halfsearchrange=20
     cuArraysMaxloc2D(r_corrBatchRaw, offsetInit, r_maxval, stream);
 
-    // Estimation of statistics
+    // estimate variance
-    // Extraction of correlation surface around the peak
+    cuEstimateVariance(r_corrBatchRaw, offsetInit, r_maxval, r_referenceBatchRaw->size, r_covValue, stream);
+
+    // estimate SNR
+    // step1: extraction of correlation surface around the peak
     cuArraysCopyExtractCorr(r_corrBatchRaw, r_corrBatchRawZoomIn, i_corrBatchZoomInValid, offsetInit, stream);
 
-    // Summation of correlation and data point values
+    // step2: summation of correlation and data point values
     cuArraysSumCorr(r_corrBatchRawZoomIn, i_corrBatchZoomInValid, r_corrBatchSum, i_corrBatchValidCount, stream);
 
 #ifdef CUAMPCOR_DEBUG
+    r_maxval->outputToFile("r_maxval", stream);
+    r_corrBatchRawZoomIn->outputToFile("r_corrBatchRawStatZoomIn", stream);
     i_corrBatchZoomInValid->outputToFile("i_corrBatchZoomInValid", stream);
     r_corrBatchSum->outputToFile("r_corrBatchSum", stream);
+    i_corrBatchValidCount->outputToFile("i_corrBatchValidCount", stream);
 #endif
 
-    // SNR
+    // step3: divide the peak value by the mean of surrounding values
     cuEstimateSnr(r_corrBatchSum, i_corrBatchValidCount, r_maxval, r_snrValue, stream);
 
-    // Variance
-    cuEstimateVariance(r_corrBatchRaw, offsetInit, r_maxval, r_covValue, stream);
-
 #ifdef CUAMPCOR_DEBUG
     offsetInit->outputToFile("i_offsetInit", stream);
-    r_maxval->outputToFile("r_maxval", stream);
+    r_snrValue->outputToFile("r_snrValue", stream);
-    r_corrBatchRawZoomIn->outputToFile("r_corrBatchRawStatZoomIn", stream);
+    r_covValue->outputToFile("r_covValue", stream);
-    i_corrBatchZoomInValid->outputToFile("i_corrBatchStatZoomInValid", stream);
 #endif
 
     // Using the approximate estimation to adjust secondary image (half search window size becomes only 4 pixels)

contrib/PyCuAmpcor/src/cuAmpcorController.cu

@@ -41,7 +41,9 @@ void cuAmpcorController::runAmpcor()
     GDALAllRegister();
     // reference and secondary images; use band=1 as default
     // TODO: selecting band
+    std::cout << "Opening reference image " << param->referenceImageName << "...\n";
     GDALImage *referenceImage = new GDALImage(param->referenceImageName, 1, param->mmapSizeInGB);
+    std::cout << "Opening secondary image " << param->secondaryImageName << "...\n";
     GDALImage *secondaryImage = new GDALImage(param->secondaryImageName, 1, param->mmapSizeInGB);
 
     cuArrays<float2> *offsetImage, *offsetImageRun;
@@ -100,9 +102,12 @@ void cuAmpcorController::runAmpcor()
         << nChunksDown << " x " << nChunksAcross  << std::endl;
 
     // iterative over chunks down
+    int message_interval = nChunksDown/10;
     for(int i = 0; i<nChunksDown; i++)
     {
-        std::cout << "Processing chunk (" << i <<", x" << ") out of " << nChunksDown << std::endl;
+        if(i%message_interval == 0)
+            std::cout << "Processing chunks (" << i+1 <<", x) - (" << std::min(nChunksDown, i+message_interval )
+                << ", x) out of " << nChunksDown << std::endl;
         // iterate over chunks across
         for(int j=0; j<nChunksAcross; j+=param->nStreams)
         {
@@ -124,12 +129,32 @@ void cuAmpcorController::runAmpcor()
     cuArraysCopyExtract(offsetImageRun, offsetImage, make_int2(0,0), streams[0]);
     cuArraysCopyExtract(snrImageRun, snrImage, make_int2(0,0), streams[0]);
     cuArraysCopyExtract(covImageRun, covImage, make_int2(0,0), streams[0]);
-    // save outputs to files
+
-    offsetImage->outputToFile(param->offsetImageName, streams[0]);
+    /* save the offsets and gross offsets */
+    // copy the offset to host
+    offsetImage->allocateHost();
+    offsetImage->copyToHost(streams[0]);
+    // construct the gross offset
+    cuArrays<float2> *grossOffsetImage = new cuArrays<float2>(param->numberWindowDown, param->numberWindowAcross);
+    grossOffsetImage->allocateHost();
+    for(int i=0; i< param->numberWindows; i++)
+        grossOffsetImage->hostData[i] = make_float2(param->grossOffsetDown[i], param->grossOffsetAcross[i]);
+
+    // check whether to merge gross offset
+    if (param->mergeGrossOffset)
+    {
+        // if merge, add the gross offsets to offset
+        for(int i=0; i< param->numberWindows; i++)
+            offsetImage->hostData[i] += grossOffsetImage->hostData[i];
+    }
+    // output both offset and gross offset
+    offsetImage->outputHostToFile(param->offsetImageName);
+    grossOffsetImage->outputHostToFile(param->grossOffsetImageName);
+    delete grossOffsetImage;
+
+    // save the snr/cov images
     snrImage->outputToFile(param->snrImageName, streams[0]);
     covImage->outputToFile(param->covImageName, streams[0]);
-    // also save the gross offsets
-    outputGrossOffsets();
 
     // Delete arrays.
     delete offsetImage;
@@ -150,19 +175,4 @@ void cuAmpcorController::runAmpcor()
     delete secondaryImage;
 
 }
-
-/**
- * Output gross offset fields
- */
-void cuAmpcorController::outputGrossOffsets()
-{
-    cuArrays<float2> *grossOffsets = new cuArrays<float2>(param->numberWindowDown, param->numberWindowAcross);
-    grossOffsets->allocateHost();
-
-    for(int i=0; i< param->numberWindows; i++)
-        grossOffsets->hostData[i] = make_float2(param->grossOffsetDown[i], param->grossOffsetAcross[i]);
-    grossOffsets->outputHostToFile(param->grossOffsetImageName);
-    delete grossOffsets;
-}
-
 // end of file

contrib/PyCuAmpcor/src/cuAmpcorController.h

@@ -27,8 +27,6 @@ public:
     ~cuAmpcorController();
     // run interface
     void runAmpcor();
-    // output gross offsets
-    void outputGrossOffsets();
 };
 #endif
 

contrib/PyCuAmpcor/src/cuAmpcorParameter.cu

@@ -59,6 +59,8 @@ cuAmpcorParameter::cuAmpcorParameter()
     useMmap = 1; // use mmap
     mmapSizeInGB = 1;
 
+    mergeGrossOffset = 0; // default to separate gross offset
+
 }
 
 /**

contrib/PyCuAmpcor/src/cuAmpcorParameter.h

@@ -112,6 +112,7 @@ public:
     int grossOffsetAcross0;     ///< gross offset static component (across)
     int *grossOffsetDown;		///< Gross offsets between reference and secondary windows (down)
     int *grossOffsetAcross;     ///< Gross offsets between reference and secondary windows (across)
+    int mergeGrossOffset;       ///< whether to merge gross offsets into the final offsets
 
     int *referenceChunkStartPixelDown;    ///< reference starting pixels for each chunk (down)
     int *referenceChunkStartPixelAcross;  ///< reference starting pixels for each chunk (across)

contrib/PyCuAmpcor/src/cuAmpcorUtil.h

@@ -91,7 +91,7 @@ void cuArraysSumCorr(cuArrays<float> *images, cuArrays<int> *imagesValid, cuArra
 void cuEstimateSnr(cuArrays<float> *corrSum, cuArrays<int> *corrValidCount, cuArrays<float> *maxval, cuArrays<float> *snrValue, cudaStream_t stream);
 
 // implemented in cuEstimateStats.cu
-void cuEstimateVariance(cuArrays<float> *corrBatchRaw, cuArrays<int2> *maxloc, cuArrays<float> *maxval, cuArrays<float3> *covValue, cudaStream_t stream);
+void cuEstimateVariance(cuArrays<float> *corrBatchRaw, cuArrays<int2> *maxloc, cuArrays<float> *maxval, int templateSize, cuArrays<float3> *covValue, cudaStream_t stream);
 
 #endif
 

contrib/PyCuAmpcor/src/cuEstimateStats.cu

@@ -46,8 +46,8 @@ void cuEstimateSnr(cuArrays<float> *corrSum, cuArrays<int> *corrValidCount, cuAr
 }
 
 // cuda kernel for cuEstimateVariance
-__global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX, const int NY,
+__global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX, const int NY, const int2* maxloc,
-    const int2* maxloc, const float* maxval, float3* covValue, const int size)
+        const float* maxval, const int templateSize, float3* covValue, const int size)
 {
 
     // Find image id.
@@ -63,7 +63,7 @@ __global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX,
     // Check if maxval is on the margin.
     if (px-1 < 0 || py-1 <0 || px + 1 >=NX || py+1 >=NY)  {
 
-        covValue[idxImage] = make_float3(99.0, 99.0, 99.0);
+        covValue[idxImage] = make_float3(99.0, 99.0, 0.0);
 
     }
     else {
@@ -78,29 +78,27 @@ __global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX,
         int idx21 = offset + (px + 1) * NY + py    ;
         int idx22 = offset + (px + 1) * NY + py + 1;
 
-        float dxx = - ( corrBatchRaw[idx21] + corrBatchRaw[idx01] - 2*corrBatchRaw[idx11] ) * 0.5;
+        // second-order derivatives
-        float dyy = - ( corrBatchRaw[idx12] + corrBatchRaw[idx10] - 2*corrBatchRaw[idx11] ) * 0.5;
+        float dxx = - ( corrBatchRaw[idx21] + corrBatchRaw[idx01] - 2.0*corrBatchRaw[idx11] );
-        float dxy = - ( corrBatchRaw[idx22] + corrBatchRaw[idx00] - corrBatchRaw[idx20] - corrBatchRaw[idx02] ) *0.25;
+        float dyy = - ( corrBatchRaw[idx12] + corrBatchRaw[idx10] - 2.0*corrBatchRaw[idx11] ) ;
+        float dxy = ( corrBatchRaw[idx22] + corrBatchRaw[idx00] - corrBatchRaw[idx20] - corrBatchRaw[idx02] ) *0.25;
 
-        float n2 = fmaxf(1 - peak, 0.0);
+        float n2 = fmaxf(1.0 - peak, 0.0);
 
-        int winSize = NX*NY;
+        dxx = dxx * templateSize;
-
+        dyy = dyy * templateSize;
-        dxx = dxx * winSize;
+        dxy = dxy * templateSize;
-        dyy = dyy * winSize;
-        dxy = dxy * winSize;
 
         float n4 = n2*n2;
         n2 = n2 * 2;
-        n4 = n4 * 0.5 * winSize;
+        n4 = n4 * 0.5 * templateSize;
 
         float u = dxy * dxy - dxx * dyy;
         float u2 = u*u;
 
+        // if the Gaussian curvature is too small
         if (fabsf(u) < 1e-2) {
-
+            covValue[idxImage] = make_float3(99.0, 99.0, 0.0);
-            covValue[idxImage] = make_float3(99.0, 99.0, 99.0);
-
         }
         else {
                 float cov_xx = (- n2 * u * dyy + n4 * ( dyy*dyy + dxy*dxy) ) / u2;
@@ -113,18 +111,19 @@ __global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX,
 
 /**
  * Estimate the variance of the correlation surface
+ * @param[in] templateSize size of reference chip
  * @param[in] corrBatchRaw correlation surface
  * @param[in] maxloc maximum location
  * @param[in] maxval maximum value
  * @param[out] covValue variance value
  * @param[in] stream cuda stream
  */
-void cuEstimateVariance(cuArrays<float> *corrBatchRaw, cuArrays<int2> *maxloc, cuArrays<float> *maxval, cuArrays<float3> *covValue, cudaStream_t stream)
+void cuEstimateVariance(cuArrays<float> *corrBatchRaw, cuArrays<int2> *maxloc, cuArrays<float> *maxval, int templateSize, cuArrays<float3> *covValue, cudaStream_t stream)
 {
     int size = corrBatchRaw->count;
     // One dimensional launching parameters to loop over every correlation surface.
     cudaKernel_estimateVar<<< IDIVUP(size, NTHREADS), NTHREADS, 0, stream>>>
-        (corrBatchRaw->devData, corrBatchRaw->height, corrBatchRaw->width, maxloc->devData, maxval->devData, covValue->devData, size);
+        (corrBatchRaw->devData, corrBatchRaw->height, corrBatchRaw->width, maxloc->devData, maxval->devData, templateSize, covValue->devData, size);
     getLastCudaError("cudaKernel_estimateVar error\n");
 }
 //end of file

contrib/PyCuAmpcor/src/cuOffset.cu

@@ -4,7 +4,9 @@
  *
  */
 
+// my module dependencies
 #include "cuAmpcorUtil.h"
+// for FLT_MAX
 #include <cfloat>
 
 // find the max between two elements
@@ -254,6 +256,7 @@ void cuDetermineSecondaryExtractOffset(cuArrays<int2> *maxLoc, cuArrays<int2> *m
     int blockspergrid=IDIVUP(maxLoc->size, threadsperblock);
     cudaKernel_determineSecondaryExtractOffset<<<blockspergrid, threadsperblock, 0, stream>>>
         (maxLoc->devData, maxLocShift->devData, maxLoc->size, xOldRange, yOldRange, xNewRange, yNewRange);
+    getLastCudaError("cuDetermineSecondaryExtractOffset");
 }
 
 // end of file

contrib/PyCuAmpcor/src/cudaUtil.h

@@ -81,7 +81,7 @@ inline int gpuDeviceInit(int devID)
     }
 
     checkCudaErrors(cudaSetDevice(devID));
-    printf("gpuDeviceInit() Using CUDA Device %d ...\n", devID);
+    printf("Using CUDA Device %d ...\n", devID);
 
     return devID;
 }

contrib/stack/topsStack/cuDenseOffsets.py

@@ -1,312 +0,0 @@
-#!/usr/bin/env python3
-
-# author: Minyan Zhong
-
-import numpy as np
-import argparse
-import os
-import isce
-import isceobj
-import shelve
-import datetime
-from isceobj.Location.Offset import OffsetField
-from iscesys.StdOEL.StdOELPy import create_writer
-#from mroipac.ampcor.DenseAmpcor import DenseAmpcor
-
-from contrib.PyCuAmpcor import PyCuAmpcor
-from grossOffsets import grossOffsets
-
-#from isceobj.Utils.denseoffsets import denseoffsets
-from isceobj.Util.decorators import use_api
-
-from pprint import pprint
-
-def createParser():
-    '''
-    Command line parser.
-    '''
-
-    parser = argparse.ArgumentParser( description='Generate offset field between two Sentinel slc')
-    parser.add_argument('-m','--reference', type=str, dest='reference', required=True,
-            help='Reference image')
-    parser.add_argument('-s', '--secondary',type=str, dest='secondary', required=True,
-            help='Secondary image')
-    parser.add_argument('-l', '--lat',type=str, dest='lat', required=False,
-           help='Latitude')
-    parser.add_argument('-L', '--lon',type=str, dest='lon', required=False,
-           help='Longitude')
-    parser.add_argument('--los',type=str, dest='los', required=False,
-           help='Line of Sight')
-    parser.add_argument('--referencexml',type=str, dest='referencexml', required=False,
-           help='Reference Image Xml File')
-    parser.add_argument('--ww', type=int, dest='winwidth', default=64,
-            help='Window Width')
-    parser.add_argument('--wh', type=int, dest='winhgt', default=64,
-            help='Window height')
-    parser.add_argument('--sw', type=int, dest='srcwidth', default=20,
-            help='Search window width')
-    parser.add_argument('--sh', type=int, dest='srchgt', default=20,
-            help='Search window height')
-    parser.add_argument('--mm', type=int, dest='margin', default=50,
-            help='Margin')
-    parser.add_argument('--kw', type=int, dest='skipwidth', default=64,
-            help='Skip across')
-    parser.add_argument('--kh', type=int, dest='skiphgt', default=64,
-            help='Skip down')
-
-    parser.add_argument('--nwa', type=int, dest='numWinAcross', default=-1,
-            help='Number of Window Across')
-    parser.add_argument('--nwd', type=int, dest='numWinDown', default=-1,
-            help='Number of Window Down')
-
-    parser.add_argument('-op','--outprefix', type=str, dest='outprefix', default='dense_ampcor',
-            help='Output prefix')
-
-    parser.add_argument('-os','--outsuffix', type=str, dest='outsuffix', default='dense_ampcor',
-            help='Output suffix')
-
-    parser.add_argument('-g','--gross', type=int, dest='gross', default=0,
-            help='Use gross offset or not')
-
-    parser.add_argument('--aa', type=int, dest='azshift', default=0,
-            help='Gross azimuth offset')
-
-    parser.add_argument('--rr', type=int, dest='rgshift', default=0,
-            help='Gross range offset')
-
-    parser.add_argument('--oo', type=int, dest='oversample', default=32,
-            help = 'Oversampling factor')
-
-    parser.add_argument('-r', '--redo', dest='redo', type=int, default=0
-       , help='To redo or not')
-
-    parser.add_argument('-drmp', '--deramp', dest='deramp', type=int, default=0
-       , help='deramp method (0: mag, 1: complex)')
-
-    parser.add_argument('-gid', '--gpuid', dest='gpuid', type=int, default=-1
-       , help='GPU ID')
-
-
-    return parser
-
-def cmdLineParse(iargs = None):
-    parser = createParser()
-    inps =  parser.parse_args(args=iargs)
-
-    return inps
-
-@use_api
-def estimateOffsetField(reference, secondary, inps=None):
-
-    ###Loading the secondary image object
-    sim = isceobj.createSlcImage()
-    sim.load(secondary+'.xml')
-    sim.setAccessMode('READ')
-    sim.createImage()
-
-
-    ###Loading the reference image object
-    sar = isceobj.createSlcImage()
-    sar.load(reference + '.xml')
-    sar.setAccessMode('READ')
-    sar.createImage()
-
-
-    width = sar.getWidth()
-    length = sar.getLength()
-
-    objOffset = PyCuAmpcor.PyCuAmpcor()
-
-    objOffset.algorithm = 0
-    objOffset.deviceID = inps.gpuid  # -1:let system find the best GPU
-    objOffset.nStreams =   2 #cudaStreams
-    objOffset.derampMethod = inps.deramp
-    print(objOffset.derampMethod)
-
-    objOffset.referenceImageName = reference + '.vrt'
-    objOffset.referenceImageHeight = length
-    objOffset.referenceImageWidth = width
-    objOffset.secondaryImageName = secondary + '.vrt'
-    objOffset.secondaryImageHeight = length
-    objOffset.secondaryImageWidth = width
-
-    print("image length:",length)
-    print("image width:",width)
-
-    objOffset.numberWindowDown = (length-2*inps.margin-2*inps.srchgt-inps.winhgt)//inps.skiphgt
-    objOffset.numberWindowAcross = (width-2*inps.margin-2*inps.srcwidth-inps.winwidth)//inps.skipwidth
-
-    if (inps.numWinDown != -1):
-        objOffset.numberWindowDown = inps.numWinDown
-
-    if (inps.numWinAcross != -1):
-        objOffset.numberWindowAcross = inps.numWinAcross
-
-
-    print("nlines: ",objOffset.numberWindowDown)
-    print("ncolumns: ",objOffset.numberWindowAcross)
-
-
-    # window size
-    objOffset.windowSizeHeight = inps.winhgt
-    objOffset.windowSizeWidth = inps.winwidth
-
-    print(objOffset.windowSizeHeight)
-    print(objOffset.windowSizeWidth)
-
-    # search range
-    objOffset.halfSearchRangeDown = inps.srchgt
-    objOffset.halfSearchRangeAcross = inps.srcwidth
-    print(inps.srchgt,inps.srcwidth)
-
-    # starting pixel
-    objOffset.referenceStartPixelDownStatic = inps.margin
-    objOffset.referenceStartPixelAcrossStatic = inps.margin
-
-    # skip size
-    objOffset.skipSampleDown = inps.skiphgt
-    objOffset.skipSampleAcross = inps.skipwidth
-
-    # oversampling
-    objOffset.corrSufaceOverSamplingMethod = 0
-    objOffset.corrSurfaceOverSamplingFactor = inps.oversample
-    #objOffset.rawDataOversamplingFactor = 4
-
-    # output filenames
-    objOffset.offsetImageName = str(inps.outprefix) + str(inps.outsuffix) + '.bip'
-    objOffset.grossOffsetImageName = str(inps.outprefix) + str(inps.outsuffix) + '_gross.bip'
-    objOffset.snrImageName = str(inps.outprefix) + str(inps.outsuffix) + '_snr.bip'
-
-    print("offsetfield: ",objOffset.offsetImageName)
-    print("gross offsetfield: ",objOffset.grossOffsetImageName)
-    print("snr: ",objOffset.snrImageName)
-
-    offsetImageName = objOffset.offsetImageName.decode('utf8')
-    #print(type(offsetImageName))
-    #print(offsetImageName)
-    #print(type(objOffset.numberWindowAcross))
-    grossOffsetImageName = objOffset.grossOffsetImageName.decode('utf8')
-    snrImageName = objOffset.snrImageName.decode('utf8')
-
-    print(offsetImageName)
-    print(inps.redo)
-    if os.path.exists(offsetImageName) and inps.redo==0:
-        print('offsetfield file exists')
-    else:
-        # generic control
-        objOffset.numberWindowDownInChunk = 5
-        objOffset.numberWindowAcrossInChunk = 5
-        objOffset.mmapSize = 16
-
-        objOffset.setupParams()
-
-        ## Set Gross Offset ###
-
-        if inps.gross == 0:
-            objOffset.setConstantGrossOffset(0, 0)
-        else:
-
-            print("Setting up grossOffset...")
-
-            objGrossOff = grossOffsets()
-
-            objGrossOff.setXSize(width)
-            objGrossOff.setYize(length)
-            objGrossOff.setMargin(inps.margin)
-            objGrossOff.setWinSizeHgt(inps.winhgt)
-            objGrossOff.setWinSizeWidth(inps.winwidth)
-            objGrossOff.setSearchSizeHgt(inps.srchgt)
-            objGrossOff.setSearchSizeWidth(inps.srcwidth)
-            objGrossOff.setSkipSizeHgt(inps.skiphgt)
-            objGrossOff.setSkipSizeWidth(inps.skipwidth)
-            objGrossOff.setLatFile(inps.lat)
-            objGrossOff.setLonFile(inps.lon)
-            objGrossOff.setLosFile(inps.los)
-            objGrossOff.setReferenceFile(inps.referencexml)
-            objGrossOff.setbTemp(inps.bTemp)
-
-
-
-            grossDown, grossAcross = objGrossOff.runGrossOffsets()
-
-            # change nan to 0
-            grossDown = np.nan_to_num(grossDown)
-            grossAcross = np.nan_to_num(grossAcross)
-
-            print("Before plotting the gross offsets (min and max): ", np.nanmin(grossDown),np.nanmax(grossDown))
-            print("Before plotting the gross offsets (min and max): ", np.rint(np.nanmin(grossDown)),np.rint(np.nanmax(grossDown)))
-
-            grossDown = np.int32(np.rint(grossDown.ravel()))
-            grossAcross = np.int32(np.rint(grossAcross.ravel()))
-
-            print(np.amin(grossDown), np.amax(grossDown))
-            print(np.amin(grossAcross), np.amax(grossAcross))
-
-            print(grossDown.shape)
-            print(grossDown.shape)
-
-            objOffset.setVaryingGrossOffset(grossDown, grossAcross)
-            #objOffset.setVaryingGrossOffset(np.zeros(shape=grossDown.shape,dtype=np.int32), np.zeros(shape=grossAcross.shape,dtype=np.int32))
-
-        # check
-        objOffset.checkPixelInImageRange()
-
-        # Run the code
-        print('Running PyCuAmpcor')
-
-        objOffset.runAmpcor()
-
-        print('Finished')
-
-        sar.finalizeImage()
-        sim.finalizeImage()
-
-        # Finalize the results
-        # offsetfield
-
-        outImg = isceobj.createImage()
-        outImg.setDataType('FLOAT')
-        outImg.setFilename(offsetImageName)
-        outImg.setBands(2)
-        outImg.scheme = 'BIP'
-        outImg.setWidth(objOffset.numberWindowAcross)
-        outImg.setLength(objOffset.numberWindowDown)
-        outImg.setAccessMode('read')
-        outImg.renderHdr()
-
-
-        # gross offsetfield
-        outImg = isceobj.createImage()
-        outImg.setDataType('FLOAT')
-        outImg.setFilename(grossOffsetImageName)
-        outImg.setBands(2)
-        outImg.scheme = 'BIP'
-        outImg.setWidth(objOffset.numberWindowAcross)
-        outImg.setLength(objOffset.numberWindowDown)
-        outImg.setAccessMode('read')
-        outImg.renderHdr()
-
-        # snr
-        snrImg = isceobj.createImage()
-        snrImg.setFilename(snrImageName)
-        snrImg.setDataType('FLOAT')
-        snrImg.setBands(1)
-        snrImg.setWidth(objOffset.numberWindowAcross)
-        snrImg.setLength(objOffset.numberWindowDown)
-        snrImg.setAccessMode('read')
-        snrImg.renderHdr()
-
-    return objOffset
-
-def main(iargs=None):
-
-    inps = cmdLineParse(iargs)
-    outDir = os.path.dirname(inps.outprefix)
-    print(inps.outprefix)
-    os.makedirs(outDir, exist_ok=True)
-
-    objOffset = estimateOffsetField(inps.reference, inps.secondary, inps)
-
-if __name__ == '__main__':
-
-    main()

contrib/stack/topsStack/grossOffsets.py

@@ -1,379 +0,0 @@
-#!/usr/bin/env python3
-# Generate grossOffsets (pixel) based on velocity field
-# Author: Minyan Zhong
-
-
-import numpy as np
-import pyproj
-import subprocess
-import isce
-import isceobj
-from iscesys.Component.ProductManager import ProductManager as PM
-import numpy as np
-from netCDF4 import Dataset
-#from mpl_toolkits.basemap import Basemap
-from osgeo import gdal
-
-from scipy.interpolate import interp2d, griddata
-
-import matplotlib.pyplot as plt
-
-
-class grossOffsets:
-
-    def __init__(self):
-        
-        self.nfig = 1
-        self.figsize = (10,10)
-
-
-        ## Antarctica Velocity File
-        self.vel_file = '/net/jokull/nobak/mzzhong/Ant_Plot/Data/antarctica_ice_velocity_900m.nc'
-        self.vProj = pyproj.Proj('+init=EPSG:3031')
-
-    def setMode(self,mode):
-
-        if mode == 'interior' or mode == 'exterior':
-            self.mode = mode
-        else:
-            raise Exception('Wrong gross offset mode')
-
-    def setLatFile(self,val):
-        self.latfile = val
-
-    def setLonFile(self,val):
-        self.lonfile = val
-
-    def setLosFile(self,val):
-        self.losfile = val
-
-    def setXSize(self,val):
-        self.XSize = val
-
-    def setYize(self,val):
-        self.YSize = val
-
-    def setMargin(self,val):
-        self.margin = val
-
-    def setWinSizeHgt(self,val):
-        self.winSizeHgt = val
-
-    def setWinSizeWidth(self,val):
-        self.winSizeWidth = val
-
-    def setSearchSizeHgt(self,val):
-        self.searchSizeHgt = val
-
-    def setSearchSizeWidth(self,val):
-        self.searchSizeWidth = val
-
-    def setSkipSizeHgt(self,val):
-        self.skipSizeHgt = val
-
-    def setSkipSizeWidth(self,val):
-        self.skipSizeWidth = val
-
-    # exterior mode
-    def setOffsetLat(self,lat):
-        self.lat = lat
-
-    def setOffsetLon(self,lon):
-        self.lon = lon
-
-    def setOffsetInc(self,inc):
-        self.inc = inc
-
-    def setOffsetAzi(self,azi):
-        self.azi = azi
-
-    def setNumWinDown(self,numWinDown):
-        self.numWinDown = numWinDown
-
-    def setNumWinAcross(self,numWinAcross):
-        self.numWinAcross = numWinAcross
-
-    def setbTemp(self,val):
-        self.bTemp = val
-
-    def setPixelSize(self,azPixelSize,rngPixelSize):
-        
-        self.azPixelSize = azPixelSize
-        self.rngPixelSize = rngPixelSize
-
-    def get_veloData(self):
-        
-        print("getting velocity data...")
-        fh=Dataset(self.vel_file,mode='r')
-        self.vx = fh.variables['vx'][:]
-        self.vy = fh.variables['vy'][:]
-        self.vx = np.flipud(self.vx)
-        self.vy = np.flipud(self.vy) 
-        self.v = np.sqrt(np.multiply(self.vx,self.vx)+np.multiply(self.vy,self.vy))
-        print(self.v.shape)
-
-        self.x0 = np.arange(-2800000,2800000,step=900)
-        self.y0 = np.arange(-2800000,2800000,step=900)+200
-
-        #x,y = np.meshgrid(self.x0,self.y0)
-
-        #from mpl_toolkits.basemap import Basemap
-        #self.AntVeloDataMap = Basemap(width=5600000,height=5600000,\
-        #                        resolution='l',projection='stere',\
-        #                        lat_ts=-71,lat_0=-90,lon_0=0)
-
-        #self.vel_lon, self.vel_lat= self.vProj(x,y,inverse="true")
-
-    def runGrossOffsets(self):
-
-        ###Pieces of information needed
-        
-        ###These pieces of information come from the output of "topo" module from ISCE
-        ### llh - size(3) - lat,lon,hgt of pixel under consideration
-        ### los - size(2) - inc, azi LOS angles
-        
-        
-        ###These pieces of information come from an external velocity product, e.g from NSIDC
-        
-        ###  vx - scalar  - Velocity in x direction at pixel under consideration
-        ###  vy - scalar  - Velocity in y direction at pixel under consideration
-        ###  vproj - string - Projection system of the velocity field
-        ###        - EPSG:3031 for Antarctica
-        ###        - EPSG:3413 for Greenland
-        
-        
-        #### The equations below describe the operations needed for a single pixel
-        #### I will use Greenland as an example. Easy to change for Antarctica by changing the coordinate system.
-        
-        ### Step 0: Set up projection transformers for ease of use
-        self.llhProj = pyproj.Proj('+init=EPSG:4326')       ##Standard lat,lon, hgt
-        self.xyzProj = pyproj.Proj('+init=EPSG:4978')    ##Standard xyz (ECEF)
-        refPt = self.vProj(0.0, 0.0, inverse=True)
-
-        print(refPt)
-        
-        ### Step 1: Set up radar image information
-        azPixelSize = self.azPixelSize
-        rngPixelSize = self.rngPixelSize
-        
-        ### Step 2: Cut the data
-       
-        print('Obtain the velocity data...')
-        self.get_veloData()
-
-        print('Extract the data to this radar scene...')
-
-        if self.mode == 'interior':
-
-            numWinDown = (self.YSize - self.margin*2 - self.searchSizeHgt*2 - self.winSizeHgt) // self.skipSizeHgt
-            numWinAcross = (self.XSize - self.margin*2 - self.searchSizeWidth*2 - self.winSizeWidth) // self.skipSizeWidth
-
-            lat = np.zeros(shape=(numWinDown,numWinAcross)) 
-            lon = np.zeros(shape=(numWinDown,numWinAcross))
-            inc = np.zeros(shape=(numWinDown,numWinAcross))
-            azi = np.zeros(shape=(numWinDown,numWinAcross))
-
-            self.centerOffsetHgt = self.searchSizeHgt+self.skipSizeHgt//2-1
-            self.centerOffsetWidth = self.searchSizeWidth+self.skipSizeWidth//2-1
-
-        elif self.mode == 'exterior':
-
-            numWinDown = self.numWinDown 
-            numWinAcross = self.numWinAcross
-
-            lat = self.lat
-            lon = self.lon
-            inc = self.inc
-            azi = self.azi
-           
-        print(numWinDown)
-        print(numWinAcross)
-
-        cut_vx = np.zeros(shape=(numWinDown,numWinAcross))
-        cut_vy = np.zeros(shape=(numWinDown,numWinAcross))
-        cut_v = np.zeros(shape=(numWinDown,numWinAcross))
-        pixel = np.zeros(shape=(numWinDown,numWinAcross))
-        line = np.zeros(shape=(numWinDown,numWinAcross))
- 
-
-        for iwin in range(numWinDown):
-
-            if self.mode == 'interior':
-
-                print('Processing line: ',iwin)
-                down = self.margin + self.skipSizeHgt * iwin  + self.centerOffsetHgt
-                off = down*self.XSize
-
-                # Warning: depend on the ENVI format. This is for BSQ
-                off2 = self.YSize * self.XSize + down*self.XSize
-    
-                start = self.margin + self.centerOffsetWidth
-                end = self.margin + self.skipSizeWidth * numWinAcross
-    
-                latline = np.memmap(filename=self.latfile,dtype='float64',offset=8*off,shape=(self.XSize))
-                lonline = np.memmap(filename=self.lonfile,dtype='float64',offset=8*off,shape=(self.XSize))
-                incline = np.memmap(filename=self.losfile,dtype='float32',offset=4*off,shape=(self.XSize))
-                aziline = np.memmap(filename=self.losfile,dtype='float32',offset=4*off2,shape=(self.XSize))
-    
-                lat[iwin,:] = latline[start:end:self.skipSizeWidth]
-                lon[iwin,:] = lonline[start:end:self.skipSizeWidth]
-                inc[iwin,:] = incline[start:end:self.skipSizeWidth]
-                azi[iwin,:] = aziline[start:end:self.skipSizeWidth]
-
-            #print(iwin,': ',lon[iwin,:])
-            #print(iwin,': ',lat[iwin,:])
-            #print(iwin,': ',inc[iwin,:])
-            #print(iwin,': ',azi[iwin,:])
-
-            #### Look up in MEaSUREs InSAR-Based Antarctica Ice Velocity Map
-
-            xyMap = pyproj.transform(self.llhProj, self.vProj, lon[iwin,:], lat[iwin,:])
-            #xyMap = self.vProj(lon[iwin,:],lat[iwin,:])
- 
-            pixel[iwin,:] = np.clip((xyMap[0]-self.x0[0])/900, 0, self.vx.shape[1]-1)
-            line[iwin,:] = np.clip((xyMap[1]-self.y0[0])/900, 0, self.vx.shape[0]-1)
-
-            pixel_int = pixel[iwin,:].astype(int)
-            line_int = line[iwin,:].astype(int)
-
-            # For Debug
-            #print(iwin,': ', 'location: ', xyMap[0],xyMap[1]) 
-            #print(iwin,': ', 'location: ', lon[iwin,:],lat[iwin,:])
-            
-            cut_vx[iwin,:] = self.vx[line_int,pixel_int]
-            cut_vy[iwin,:] = self.vy[line_int,pixel_int]
-            cut_v[iwin,:] = np.sqrt(np.multiply(cut_vx[iwin,:],cut_vx[iwin,:]),np.multiply(cut_vy[iwin,:],cut_vy[iwin,:]))
-
-        ## Interpolate offsetfield
-        # Mask out invalid value based on the value of lat (or lon) (only work for polar region)
-        # Mask out zero velocity
-        print('Interpolating velocity field...')
-        valid = np.logical_and(lat!=0,cut_v!=0)
-        x0=np.arange(numWinAcross)
-        y0=np.arange(numWinDown)
-        xx,yy=np.meshgrid(x0,y0)
-
-        grid_x,grid_y=[grid.ravel() for grid in np.meshgrid(x0,y0)]
-        
-        points = np.column_stack((xx[valid],yy[valid]))
-        print(points.shape)
-        in_dat = cut_vx[valid]
-        cut_vx_new = griddata(points, in_dat, (grid_x, grid_y), method='linear')
-        in_dat = cut_vy[valid]
-        cut_vy_new = griddata(points, in_dat, (grid_x, grid_y), method='linear')
-        
-        cut_vx_new = cut_vx_new.reshape(numWinDown,numWinAcross)
-        cut_vy_new = cut_vy_new.reshape(numWinDown,numWinAcross)
-        
-        # mask out invalid values at margin
-        cut_vx_new[lat==0] = np.nan
-        cut_vy_new[lat==0] = np.nan
-        
-        cut_v_new = np.sqrt(np.multiply(cut_vx_new,cut_vx_new),np.multiply(cut_vy_new,cut_vy_new))
-        
-        print(cut_v_new)
-        print(cut_v_new.shape)
-
-        ##########
-        #fig=plt.figure(10,figsize=(10,10))
-        #ax = fig.add_subplot(111)
-        #print(cut_v.shape)
-        #ax.imshow(np.clip(cut_v_new,0,1000),cmap=plt.cm.viridis)
-        #fig.savefig('10.png',format='png')
-
-        ### Step 3: Convert XY velocity to EN velocity (clockwise rotation)
-
-        print('Coverting XY to EN...')
-        lonr = np.radians(lon - refPt[0])
-
-        cut_ve = np.multiply(cut_vx_new, np.cos(lonr)) - np.multiply(cut_vy_new, np.sin(lonr))
-        cut_vn = np.multiply(cut_vy_new, np.cos(lonr)) + np.multiply(cut_vx_new, np.sin(lonr))
-        
-        print('Polar stereographic velocity: ', [cut_vx, cut_vy])
-        print('Local ENU velocity: ', [cut_ve, cut_vn])
- 
-        ####Step 4: Convert EN velocity to rng and azimuth
-        
-        #Local los and azi vector in ENU coordinate
-
-        print(' Coverting EN to rdr...')
-        incr = np.radians(inc)
-        azir = np.radians(azi)
-        losr = np.radians(azi-90.0)
-
-        losenu=[ np.multiply(np.sin(incr),np.cos(losr)),
-                 np.multiply(np.sin(incr),np.sin(losr)),
-                 -np.cos(incr) ]
-        
-        azienu=[ np.cos(azir),
-                 np.sin(azir),
-                 0.0 ]
-        
-        grossRangeOffset = (self.bTemp/365.25) * (cut_ve * losenu[0] + cut_vn * losenu[1])/ rngPixelSize
-        
-        grossAzimuthOffset = (self.bTemp/365.25) * (cut_ve * azienu[0] + cut_vn * azienu[1]) / azPixelSize
-        
-        print('Gross azimuth offset: ', grossAzimuthOffset)
-        print('Gross range offset: ', grossRangeOffset)
-
-        # Float
-        fig=plt.figure(21,figsize=(16,9))
-        #ax = fig.add_subplot(121)
-        ax = fig.add_axes([0.05,0.05,0.4,0.9])
-        ax.set_title('gross azimuth offset',fontsize=15)
-        print(grossRangeOffset.shape)
-        cax = ax.imshow(grossAzimuthOffset,cmap=plt.cm.coolwarm)
-        cbar = fig.colorbar(cax,fraction=0.035,pad=0.04,ticks=np.arange(np.rint(np.nanmin(grossAzimuthOffset)),np.rint(np.nanmax(grossAzimuthOffset))+0.1))
-        cbar.set_label("pixel",fontsize=15)
-
-        #ax = fig.add_subplot(122)
-        ax = fig.add_axes([0.55,0.05,0.4,0.9])
-        ax.set_title('gross range offset',fontsize=15)
-        print(grossRangeOffset.shape)
-        cax = ax.imshow(grossRangeOffset,cmap=plt.cm.coolwarm)
-
-        cbar = fig.colorbar(cax,fraction=0.035,pad=0.04,ticks=np.arange(np.rint(np.nanmin(grossRangeOffset)),np.rint(np.nanmax(grossRangeOffset))+0.1))
-        cbar.set_label("pixel",fontsize=15)
-
-        #fig.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.1)
-        #fig.tight_layout()
-        fig.savefig('21.png',format='png')
-
-
-
-
-        ## Round to integer
-        fig=plt.figure(22,figsize=(16,9))
-        #ax = fig.add_subplot(121)
-        ax = fig.add_axes([0.05,0.05,0.4,0.9])
-        ax.set_title('gross azimuth offset',fontsize=15)
-        print(grossRangeOffset.shape)
-        cax = ax.imshow(np.rint(grossAzimuthOffset),cmap=plt.cm.coolwarm)
-        cbar = fig.colorbar(cax,fraction=0.035,pad=0.04,ticks=np.arange(np.rint(np.nanmin(grossAzimuthOffset)),np.rint(np.nanmax(grossAzimuthOffset))+0.1))
-        cbar.set_label("pixel",fontsize=15)
-
-        #ax = fig.add_subplot(122)
-        ax = fig.add_axes([0.55,0.05,0.4,0.9])
-        ax.set_title('gross range offset',fontsize=15)
-        print(grossRangeOffset.shape)
-        cax = ax.imshow(np.rint(grossRangeOffset),cmap=plt.cm.coolwarm)
-
-        print("Before plotting the gross offsets (min and max): ", np.rint(np.nanmin(grossAzimuthOffset)),np.rint(np.nanmax(grossAzimuthOffset)))
-
-        cbar = fig.colorbar(cax,fraction=0.035,pad=0.04,ticks=np.arange(np.rint(np.nanmin(grossRangeOffset)),np.rint(np.nanmax(grossRangeOffset))+0.1))
-        cbar.set_label("pixel",fontsize=15)
-
-        #fig.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.1)
-        #fig.tight_layout()
-        fig.savefig('22.png',format='png')
-
-        return grossAzimuthOffset, grossRangeOffset
-
-def main():
-
-    grossObj = grossOffsets()
-    grossObj.runGrossOffsets()
-
-if __name__=='__main__':
-    main()
-