diff --git a/components/isceobj/Alos2Proc/runDenseOffset.py b/components/isceobj/Alos2Proc/runDenseOffset.py index 1c87b9c..e3d3853 100644 --- a/components/isceobj/Alos2Proc/runDenseOffset.py +++ b/components/isceobj/Alos2Proc/runDenseOffset.py @@ -70,7 +70,7 @@ def runDenseOffset(self): #get water body mask wbd0=np.memmap('wbd.rdr', dtype=np.int8, mode='r', shape=(length0, width0)) - wbd0=wbd0[0+self._insar.offsetImageTopoffset:length0:self.offsetSkipHeight, + wbd0=wbd0[0+self._insar.offsetImageTopoffset:length0:self.offsetSkipHeight, 0+self._insar.offsetImageLeftoffset:width0:self.offsetSkipWidth] wbd = np.zeros((length+100, width+100), dtype=np.int8) wbd[0:wbd0.shape[0], 0:wbd0.shape[1]]=wbd0 @@ -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.nStreams = 2 - #original ampcor program in roi_pac uses phase gradient to deramp - objOffset.derampMethod = 2 + objOffset.derampMethod = 1 # 1=linear phase ramp, 0=take mag, 2=skip 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 - objOffset.referenceStartPixelAcrossStatic = self.offsetWindowWidth//2 + # adjust the margin + margin = max(margin, abs(objOffset.grossOffsetAcrossStatic), abs(objOffset.grossOffsetDownStatic)) - objOffset.numberWindowDown = (m.length - 2*self.offsetSearchWindowHeight - self.offsetWindowHeight) // self.offsetSkipHeight - objOffset.numberWindowAcross = (m.width - 2*self.offsetSearchWindowWidth - self.offsetWindowWidth) // self.offsetSkipWidth + # set the starting pixel of the first reference window + objOffset.referenceStartPixelDownStatic = margin + self.offsetSearchWindowHeight + objOffset.referenceStartPixelAcrossStatic = margin + self.offsetSearchWindowWidth - # generic control - objOffset.numberWindowDownInChunk = 8 - objOffset.numberWindowAcrossInChunk = 8 + # find out the total number of windows + objOffset.numberWindowDown = (m.length - 2*margin - 2*self.offsetSearchWindowHeight - self.offsetWindowHeight) // self.offsetSkipHeight + 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 - self._insar.offsetImageLeftoffset = objOffset.halfSearchRangeAcross + # location of the center of the first reference window + 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')) - offsetBIL.astype(np.float32).tofile(objOffset.offsetImageName.decode('utf-8')) + os.remove(objOffset.offsetImageName) + 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 diff --git a/components/isceobj/TopsProc/runDenseOffsets.py b/components/isceobj/TopsProc/runDenseOffsets.py index 75023cf..5b80ea2 100644 --- a/components/isceobj/TopsProc/runDenseOffsets.py +++ b/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 - objOffset.nStreams = 2 + # deramping method: 0 to take magnitude (fixed for Tops) 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 - objOffset.numberWindowAcross = (width-100-self.winwidth)//self.skipwidth + # adjust the margin + 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 - objOffset.referenceStartPixelAcrossStatic = 50 - + # set the skip distance between windows 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 - objOffset.numberWindowDownInChunk = 10 - objOffset.numberWindowAcrossInChunk = 10 + ### gpu control + objOffset.deviceID = 0 + 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 - ### Configure dense Ampcor object + ### print the settings 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') - - objOffset.checkPixelInImageRange() + # run ampcor 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 - self._insar.offset_left = 50 + # the center of the first reference window + 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__' : diff --git a/contrib/PyCuAmpcor/README.md b/contrib/PyCuAmpcor/README.md index 6632ea1..62ce22e 100644 --- a/contrib/PyCuAmpcor/README.md +++ b/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** diff --git a/contrib/PyCuAmpcor/examples/cuDenseOffsets.py b/contrib/PyCuAmpcor/examples/cuDenseOffsets.py index 1166f2b..96a9413 100755 --- a/contrib/PyCuAmpcor/examples/cuDenseOffsets.py +++ b/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 diff --git a/contrib/PyCuAmpcor/examples/grossOffsets.py b/contrib/PyCuAmpcor/examples/grossOffsets.py new file mode 100755 index 0000000..c9d2b80 --- /dev/null +++ b/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() diff --git a/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx b/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx index 9db1c8f..8af89d2 100644 --- a/contrib/PyCuAmpcor/src/PyCuAmpcor.pyx +++ b/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 = a.encode() - @property - def referenceImageName(self): - return self.c_cuAmpcor.param.referenceImageName - @referenceImageName.setter - def referenceImageName(self, str a): - self.c_cuAmpcor.param.referenceImageName = 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 = 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") @@ -449,4 +453,4 @@ cdef class PyCuAmpcor(object): self.c_cuAmpcor.runAmpcor() -# end of file \ No newline at end of file +# end of file diff --git a/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu b/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu index a3d4e3e..c7c894d 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorChunk.cu +++ b/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 - // Extraction of correlation surface around the peak + // estimate variance + 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_corrBatchRawZoomIn->outputToFile("r_corrBatchRawStatZoomIn", stream); - i_corrBatchZoomInValid->outputToFile("i_corrBatchStatZoomInValid", stream); + r_snrValue->outputToFile("r_snrValue", stream); + r_covValue->outputToFile("r_covValue", stream); #endif // Using the approximate estimation to adjust secondary image (half search window size becomes only 4 pixels) diff --git a/contrib/PyCuAmpcor/src/cuAmpcorController.cu b/contrib/PyCuAmpcor/src/cuAmpcorController.cu index fc9a086..625a6a4 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorController.cu +++ b/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 *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; inStreams) { @@ -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 *grossOffsetImage = new cuArrays(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 *grossOffsets = new cuArrays(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 diff --git a/contrib/PyCuAmpcor/src/cuAmpcorController.h b/contrib/PyCuAmpcor/src/cuAmpcorController.h index 9c0c77e..77973b6 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorController.h +++ b/contrib/PyCuAmpcor/src/cuAmpcorController.h @@ -27,8 +27,6 @@ public: ~cuAmpcorController(); // run interface void runAmpcor(); - // output gross offsets - void outputGrossOffsets(); }; #endif diff --git a/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu b/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu index cc9e348..d1451b4 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu +++ b/contrib/PyCuAmpcor/src/cuAmpcorParameter.cu @@ -59,6 +59,8 @@ cuAmpcorParameter::cuAmpcorParameter() useMmap = 1; // use mmap mmapSizeInGB = 1; + mergeGrossOffset = 0; // default to separate gross offset + } /** diff --git a/contrib/PyCuAmpcor/src/cuAmpcorParameter.h b/contrib/PyCuAmpcor/src/cuAmpcorParameter.h index 862284a..9116ac1 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorParameter.h +++ b/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) diff --git a/contrib/PyCuAmpcor/src/cuAmpcorUtil.h b/contrib/PyCuAmpcor/src/cuAmpcorUtil.h index 752bc88..55961f8 100644 --- a/contrib/PyCuAmpcor/src/cuAmpcorUtil.h +++ b/contrib/PyCuAmpcor/src/cuAmpcorUtil.h @@ -91,7 +91,7 @@ void cuArraysSumCorr(cuArrays *images, cuArrays *imagesValid, cuArra void cuEstimateSnr(cuArrays *corrSum, cuArrays *corrValidCount, cuArrays *maxval, cuArrays *snrValue, cudaStream_t stream); // implemented in cuEstimateStats.cu -void cuEstimateVariance(cuArrays *corrBatchRaw, cuArrays *maxloc, cuArrays *maxval, cuArrays *covValue, cudaStream_t stream); +void cuEstimateVariance(cuArrays *corrBatchRaw, cuArrays *maxloc, cuArrays *maxval, int templateSize, cuArrays *covValue, cudaStream_t stream); #endif diff --git a/contrib/PyCuAmpcor/src/cuEstimateStats.cu b/contrib/PyCuAmpcor/src/cuEstimateStats.cu index 7973f7b..18412bf 100644 --- a/contrib/PyCuAmpcor/src/cuEstimateStats.cu +++ b/contrib/PyCuAmpcor/src/cuEstimateStats.cu @@ -46,8 +46,8 @@ void cuEstimateSnr(cuArrays *corrSum, cuArrays *corrValidCount, cuAr } // cuda kernel for cuEstimateVariance -__global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX, const int NY, - const int2* maxloc, const float* maxval, float3* covValue, const int size) +__global__ void cudaKernel_estimateVar(const float* corrBatchRaw, const int NX, const int NY, const int2* maxloc, + 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; - float dyy = - ( corrBatchRaw[idx12] + corrBatchRaw[idx10] - 2*corrBatchRaw[idx11] ) * 0.5; - float dxy = - ( corrBatchRaw[idx22] + corrBatchRaw[idx00] - corrBatchRaw[idx20] - corrBatchRaw[idx02] ) *0.25; + // second-order derivatives + float dxx = - ( corrBatchRaw[idx21] + corrBatchRaw[idx01] - 2.0*corrBatchRaw[idx11] ); + 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 * winSize; - dyy = dyy * winSize; - dxy = dxy * winSize; + dxx = dxx * templateSize; + dyy = dyy * templateSize; + dxy = dxy * templateSize; 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, 99.0); - + covValue[idxImage] = make_float3(99.0, 99.0, 0.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 *corrBatchRaw, cuArrays *maxloc, cuArrays *maxval, cuArrays *covValue, cudaStream_t stream) +void cuEstimateVariance(cuArrays *corrBatchRaw, cuArrays *maxloc, cuArrays *maxval, int templateSize, cuArrays *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 \ No newline at end of file +//end of file diff --git a/contrib/PyCuAmpcor/src/cuOffset.cu b/contrib/PyCuAmpcor/src/cuOffset.cu index 0121204..ac77f13 100644 --- a/contrib/PyCuAmpcor/src/cuOffset.cu +++ b/contrib/PyCuAmpcor/src/cuOffset.cu @@ -4,7 +4,9 @@ * */ +// my module dependencies #include "cuAmpcorUtil.h" +// for FLT_MAX #include // find the max between two elements @@ -254,6 +256,7 @@ void cuDetermineSecondaryExtractOffset(cuArrays *maxLoc, cuArrays *m int blockspergrid=IDIVUP(maxLoc->size, threadsperblock); cudaKernel_determineSecondaryExtractOffset<<>> (maxLoc->devData, maxLocShift->devData, maxLoc->size, xOldRange, yOldRange, xNewRange, yNewRange); + getLastCudaError("cuDetermineSecondaryExtractOffset"); } // end of file diff --git a/contrib/PyCuAmpcor/src/cudaUtil.h b/contrib/PyCuAmpcor/src/cudaUtil.h index 5a2629c..cc174e0 100644 --- a/contrib/PyCuAmpcor/src/cudaUtil.h +++ b/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; } @@ -120,4 +120,4 @@ inline void gpuDeviceList() } #endif //__CUDAUTIL_H -//end of file \ No newline at end of file +//end of file diff --git a/contrib/stack/topsStack/cuDenseOffsets.py b/contrib/stack/topsStack/cuDenseOffsets.py deleted file mode 100755 index 1f99f83..0000000 --- a/contrib/stack/topsStack/cuDenseOffsets.py +++ /dev/null @@ -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() diff --git a/contrib/stack/topsStack/grossOffsets.py b/contrib/stack/topsStack/grossOffsets.py deleted file mode 100755 index 4d21da7..0000000 --- a/contrib/stack/topsStack/grossOffsets.py +++ /dev/null @@ -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() -