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/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() -