ISCE_INSAR/components/zerodop/GPUtopozero/test_script/gpu-topo.py

416 lines
14 KiB
Python
Raw Normal View History

2019-01-16 19:40:08 +00:00
#!/usr/bin/env python3
import argparse
import isce
import isceobj
import numpy as np
import shelve
import os
import datetime
import isceobj.Image as IF
from isceobj.Constants import SPEED_OF_LIGHT
from isceobj.Util.Poly2D import Poly2D
from iscesys import DateTimeUtil as DTU
from iscesys.Component.ProductManager import ProductManager
import gpu_topozero
def cmdLineParse():
'''
Command line parser.
'''
parser = argparse.ArgumentParser( description='Create DEM simulation for merged images')
parser.add_argument('-a','--alks', dest='alks', type=int, default=1,
help = 'Number of azimuth looks')
parser.add_argument('-r','--rlks', dest='rlks', type=int, default=1,
help = 'Number of range looks')
parser.add_argument('-d', '--dem', dest='dem', type=str, required=True,
help = 'Input DEM to use')
parser.add_argument('-m', '--reference', dest='reference', type=str, required=True,
help = 'Dir with reference frame')
2019-01-16 19:40:08 +00:00
parser.add_argument('-o', '--output', dest='outdir', type=str, required=True,
help = 'Output directory')
parser.add_argument('-n','--native', dest='nativedop', action='store_true',
default=False, help='Products in native doppler geometry instead of zero doppler')
parser.add_argument('-l','--legendre', dest='legendre', action='store_true',
default=False, help='Use legendre interpolation instead of hermite')
parser.add_argument('-f', '--full', dest='full', action='store_true',
default=False, help='Generate all topo products - masks etc')
parser.add_argument('-s', '--sentinel', dest='sntl1a', action='store_true',
default=False, help='Designate input as Sentinel data')
return parser.parse_args()
class Dummy(object):
pass
def runGPUTopo(info, demImage, dop=None, nativedop=False, legendre=False, full=False):
from zerodop.topozero import createTopozero
from isceobj.Planet.Planet import Planet
from gpu_topozero import PyTopozero
from isceobj import Constants as CN
if not os.path.isdir(info.outdir):
os.mkdir(info.outdir)
# Random vars
r0 = info.rangeFirstSample + ((info.numberRangeLooks - 1)/2) * info.slantRangePixelSpacing
tbef = info.sensingStart + datetime.timedelta(seconds = ((info.numberAzimuthLooks - 1) /2) / info.prf)
pegHdg = np.radians(info.orbit.getENUHeading(tbef))
width = info.width // info.numberRangeLooks
length = info.length // info.numberAzimuthLooks
dr = info.slantRangePixelSpacing*info.numberRangeLooks
if legendre:
omethod = 2 # LEGENDRE INTERPOLATION
else:
omethod = 0 # HERMITE INTERPOLATION
# Images
demImage.setCaster('read','FLOAT')
demImage.createImage()
if nativedop and (dop is not None):
try:
coeffs = dop._coeffs
except:
coeffs = dop
polyDoppler = Poly2D()
polyDoppler.setWidth(width)
polyDoppler.setLength(length)
polyDoppler.initPoly(rangeOrder = len(coeffs)-1, azimuthOrder=0, coeffs=[coeffs])
else:
print('Zero doppler')
#doppler = None
polyDoppler = Poly2D(name=frame.name+'_dopplerPoly')
polyDoppler.setWidth(width)
polyDoppler.setLength(length)
polyDoppler.setNormRange(1.0)
polyDoppler.setNormAzimuth(1.0)
polyDoppler.setMeanRange(0.0)
polyDoppler.setMeanAzimuth(0.0)
polyDoppler.initPoly(rangeOrder=0, azimuthOrder=0, coeffs=[[0.0]])
polyDoppler.createPoly2D()
slantRangeImage = Poly2D()
slantRangeImage.setWidth(width)
slantRangeImage.setLength(length)
slantRangeImage.setNormRange(1.0)
slantRangeImage.setNormAzimuth(1.0)
slantRangeImage.setMeanRange(0.0)
slantRangeImage.setMeanAzimuth(0.0)
slantRangeImage.initPoly(rangeOrder=1,azimuthOrder=0, coeffs=[[r0,dr]])
slantRangeImage.createPoly2D()
latImage = IF.createImage()
accessMode = 'write'
dataType = 'DOUBLE'
latImage.initImage(os.path.join(info.outdir, 'lat.rdr'),accessMode,width,dataType)
latImage.createImage()
lonImage = IF.createImage()
lonImage.initImage(os.path.join(info.outdir, 'lon.rdr'),accessMode,width,dataType)
lonImage.createImage()
losImage = IF.createImage()
dataType = 'FLOAT'
bands = 2
scheme = 'BIL'
losImage.initImage(os.path.join(info.outdir, 'los.rdr'),accessMode,width,dataType,bands=bands,scheme=scheme)
losImage.setCaster('write','DOUBLE')
losImage.createImage()
heightImage = IF.createImage()
dataType = 'DOUBLE'
heightImage.initImage(os.path.join(info.outdir, 'z.rdr'),accessMode,width,dataType)
heightImage.createImage()
if full:
incImage = IF.createImage()
dataType = 'FLOAT'
incImage.initImage(os.path.join(info.outdir, 'inc.rdr'),accessMode,width,dataType,bands=bands,scheme=scheme)
incImage.createImage()
incImagePtr = incImage.getImagePointer()
maskImage = IF.createImage()
dataType = 'BYTE'
bands = 1
maskImage.initImage(os.path.join(info.outdir, 'mask.rdr'),accessMode,width,dataType,bands=bands,scheme=scheme)
maskImage.createImage()
maskImagePtr = maskImage.getImagePointer()
else:
incImagePtr = 0
maskImagePtr = 0
elp = Planet(pname='Earth').ellipsoid
topo = PyTopozero()
topo.set_firstlat(demImage.getFirstLatitude())
topo.set_firstlon(demImage.getFirstLongitude())
topo.set_deltalat(demImage.getDeltaLatitude())
topo.set_deltalon(demImage.getDeltaLongitude())
topo.set_major(elp.a)
topo.set_eccentricitySquared(elp.e2)
topo.set_rSpace(info.slantRangePixelSpacing)
topo.set_r0(r0)
topo.set_pegHdg(pegHdg)
topo.set_prf(info.prf)
topo.set_t0(DTU.seconds_since_midnight(tbef))
topo.set_wvl(info.radarWavelength)
topo.set_thresh(.05)
topo.set_demAccessor(demImage.getImagePointer())
topo.set_dopAccessor(polyDoppler.getPointer())
topo.set_slrngAccessor(slantRangeImage.getPointer())
topo.set_latAccessor(latImage.getImagePointer())
topo.set_lonAccessor(lonImage.getImagePointer())
topo.set_losAccessor(losImage.getImagePointer())
topo.set_heightAccessor(heightImage.getImagePointer())
topo.set_incAccessor(incImagePtr)
topo.set_maskAccessor(maskImagePtr)
topo.set_numIter(25)
topo.set_idemWidth(demImage.getWidth())
topo.set_idemLength(demImage.getLength())
topo.set_ilrl(info.lookSide)
topo.set_extraIter(10)
topo.set_length(length)
topo.set_width(width)
topo.set_nRngLooks(info.numberRangeLooks)
topo.set_nAzLooks(info.numberAzimuthLooks)
topo.set_demMethod(5) # BIQUINTIC METHOD
topo.set_orbitMethod(omethod)
# Need to simplify orbit stuff later
nvecs = len(info.orbit.stateVectors.list)
topo.set_orbitNvecs(nvecs)
topo.set_orbitBasis(1) # Is this ever different?
topo.createOrbit() # Initializes the empty orbit to the right allocated size
count = 0
for sv in info.orbit.stateVectors.list:
td = DTU.seconds_since_midnight(sv.getTime())
pos = sv.getPosition()
vel = sv.getVelocity()
topo.set_orbitVector(count,td,pos[0],pos[1],pos[2],vel[0],vel[1],vel[2])
count += 1
topo.runTopo()
latImage.addDescription('Pixel-by-pixel latitude in degrees.')
latImage.finalizeImage()
latImage.renderHdr()
lonImage.addDescription('Pixel-by-pixel longitude in degrees.')
lonImage.finalizeImage()
lonImage.renderHdr()
heightImage.addDescription('Pixel-by-pixel height in meters.')
heightImage.finalizeImage()
heightImage.renderHdr()
descr = '''Two channel Line-Of-Sight geometry image (all angles in degrees). Represents vector drawn from target to platform.
Channel 1: Incidence angle measured from vertical at target (always +ve).
Channel 2: Azimuth angle measured from North in Anti-clockwise direction.'''
losImage.setImageType('bil')
losImage.addDescription(descr)
losImage.finalizeImage()
losImage.renderHdr()
demImage.finalizeImage()
if full:
descr = '''Two channel angle file.
Channel 1: Angle between ray to target and the vertical at the sensor
Channel 2: Local incidence angle accounting for DEM slope at target'''
incImage.addDescription(descr)
incImage.finalizeImage()
incImage.renderHdr()
descr = 'Radar shadow-layover mask. 1 - Radar Shadow. 2 - Radar Layover. 3 - Both.'
maskImage.addDescription(descr)
maskImage.finalizeImage()
maskImage.renderHdr()
if slantRangeImage:
try:
slantRangeImage.finalizeImage()
except:
pass
def runTopo(info, demImage, dop=None,
nativedop=False, legendre=False, full=False):
from zerodop.topozero import createTopozero
from isceobj.Planet.Planet import Planet
if not os.path.isdir(info.outdir):
os.mkdir(info.outdir)
#####Run Topo
planet = Planet(pname='Earth')
topo = createTopozero()
topo.slantRangePixelSpacing = info.slantRangePixelSpacing
topo.prf = info.prf
topo.radarWavelength = info.radarWavelength
topo.orbit = info.orbit
topo.width = info.width // info.numberRangeLooks
topo.length = info.length //info.numberAzimuthLooks
topo.wireInputPort(name='dem', object=demImage)
topo.wireInputPort(name='planet', object=planet)
topo.numberRangeLooks = info.numberRangeLooks
topo.numberAzimuthLooks = info.numberAzimuthLooks
topo.lookSide = info.lookSide
topo.sensingStart = info.sensingStart + datetime.timedelta(seconds = ((info.numberAzimuthLooks - 1) /2) / info.prf)
topo.rangeFirstSample = info.rangeFirstSample + ((info.numberRangeLooks - 1)/2) * info.slantRangePixelSpacing
topo.demInterpolationMethod='BIQUINTIC'
if legendre:
topo.orbitInterpolationMethod = 'LEGENDRE'
topo.latFilename = os.path.join(info.outdir, 'lat.rdr')
topo.lonFilename = os.path.join(info.outdir, 'lon.rdr')
topo.losFilename = os.path.join(info.outdir, 'los.rdr')
topo.heightFilename = os.path.join(info.outdir, 'z.rdr')
if full:
topo.incFilename = os.path.join(info.outdir, 'inc.rdr')
topo.maskFilename = os.path.join(info.outdir, 'mask.rdr')
if nativedop and (dop is not None):
try:
coeffs = dop._coeffs
except:
coeffs = dop
doppler = Poly2D()
doppler.setWidth(info.width // info.numberRangeLooks)
doppler.setLength(info.length // info.numberAzimuthLooks)
doppler.initPoly(rangeOrder = len(coeffs)-1, azimuthOrder=0, coeffs=[coeffs])
else:
print('Zero doppler')
doppler = None
topo.polyDoppler = doppler
topo.topo()
return
def runSimamp(outdir, hname='z.rdr'):
from iscesys.StdOEL.StdOELPy import create_writer
#####Run simamp
stdWriter = create_writer("log","",True,filename='sim.log')
objShade = isceobj.createSimamplitude()
objShade.setStdWriter(stdWriter)
hgtImage = isceobj.createImage()
hgtImage.load(os.path.join(outdir, hname) + '.xml')
hgtImage.setAccessMode('read')
hgtImage.createImage()
simImage = isceobj.createImage()
simImage.setFilename(os.path.join(outdir, 'simamp.rdr'))
simImage.dataType = 'FLOAT'
simImage.setAccessMode('write')
simImage.setWidth(hgtImage.getWidth())
simImage.createImage()
objShade.simamplitude(hgtImage, simImage, shade=3.0)
simImage.renderHdr()
hgtImage.finalizeImage()
simImage.finalizeImage()
def extractInfo(frame, inps):
'''
Extract relevant information only.
'''
info = Dummy()
ins = frame.getInstrument()
info.sensingStart = frame.getSensingStart()
info.lookSide = frame.instrument.platform.pointingDirection
info.rangeFirstSample = frame.startingRange
info.numberRangeLooks = inps.rlks
info.numberAzimuthLooks = inps.alks
fsamp = frame.rangeSamplingRate
info.slantRangePixelSpacing = 0.5 * SPEED_OF_LIGHT / fsamp
info.prf = frame.PRF
info.radarWavelength = frame.radarWavelegth
info.orbit = frame.getOrbit()
info.width = frame.getNumberOfSamples()
info.length = frame.getNumberOfLines()
info.sensingStop = frame.getSensingStop()
info.outdir = inps.outdir
return info
def extractInfoFromS1A(frame, inps):
'''
Extract relevant information only.
'''
info = Dummy()
info.sensingStart = frame.bursts[0].sensingStart
info.lookSide = -1
info.rangeFirstSample = frame.bursts[0].startingRange
info.numberRangeLooks = inps.rlks
info.numberAzimuthLooks = inps.alks
info.slantRangePixelSpacing = frame.bursts[0].rangePixelSize
info.prf = 1. / frame.bursts[0].azimuthTimeInterval
info.radarWavelength = frame.bursts[0].radarWavelength
info.orbit = frame.bursts[0].orbit
info.width = frame.bursts[0].numberOfSamples * 3
info.length = frame.bursts[0].numberOfLines * len(frame.bursts)
info.sensingStop = frame.bursts[-1].sensingStop
info.outdir = inps.outdir
return info
if __name__ == '__main__':
inps = cmdLineParse()
if (inps.sntl1a):
pm = ProductManager()
pm.configure()
frame = pm.loadProduct(inps.reference)
2019-01-16 19:40:08 +00:00
doppler = [0.]
else:
db = shelve.open(os.path.join(inps.reference, 'data'))
2019-01-16 19:40:08 +00:00
frame = db['frame']
try:
doppler = db['doppler']
except:
doppler = frame._dopplerVsPixel
db.close()
####Setup dem
demImage = isceobj.createDemImage()
demImage.load(inps.dem + '.xml')
demImage.setAccessMode('read')
if (inps.sntl1a):
info = extractInfoFromS1A(frame,inps)
else:
info = extractInfo(frame, inps)
# To revert: delete 'GPU' from method call name
runGPUTopo(info,demImage,dop=doppler,
nativedop=inps.nativedop, legendre=inps.legendre,
full=inps.full)
#runSimamp(inps.outdir)