ISCE_INSAR/applications/dpmApp.py

645 lines
24 KiB
Python

#!/usr/bin/env python3
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Copyright 2010 California Institute of Technology. ALL RIGHTS RESERVED.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# United States Government Sponsorship acknowledged. This software is subject to
# U.S. export control laws and regulations and has been classified as 'EAR99 NLR'
# (No [Export] License Required except when exporting to an embargoed country,
# end user, or in support of a prohibited end use). By downloading this software,
# the user agrees to comply with all applicable U.S. export laws and regulations.
# The user has the responsibility to obtain export licenses, or other export
# authority as may be required before exporting this software to any 'EAR99'
# embargoed foreign country or citizen of those countries.
#
# Author: Giangi Sacco
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
from __future__ import print_function
import time
import os
import isce
import isceobj
import iscesys
from iscesys.Compatibility import Compatibility
from isceobj.Pause import pause
from iscesys.ImageUtil.ImageUtil import ImageUtil as IU
from iscesys.Component.Application import Application
from isce.applications.insarApp import _InsarBase, logger
class Dpm(_InsarBase):
"""Dpm Application class:
Implements Dpm processing flow for a pair of scenes from
sensor raw data to geocoded correlation.
"""
def __init__(self):
super(Dpm, self).__init__()
## This indicates something has gone wrong, I must delete geocode.
del self.runGeocode
## extends _InsarBase_steps, but not in the same was as main
def _steps(self):
super(Dpm, self)._steps()
# Geocode
self.step('geocodecorifg', func=self.geocodeCorIfg)
self.step('geocodecor4rlks', func=self.geocodeCor4rlks)
self.step('renderProcDoc', func=self.renderProcDoc)
self.step('timerend', local='timeEnd',func=time.time)
self.step('logtime', func=logger.info,
delayed_args = (" 'Total Time: %i seconds'%(timeEnd-timeStart)",)
)
return None
def renderProcDoc(self):
self.insar.procDoc.renderXml()
def coherence(self):
self.runCoherence(self.correlation_method)
self.renderProcDoc()
def geocodeCorIfg(self):
corFilename = self.insar.coherenceFilename
corintFilename = corFilename.replace('.cor','.corint')
widthInt = self.insar.resampIntImage.width
rmg_to_cmplx(corFilename,corintFilename,widthInt)
corintGeocodeFilename = corintFilename+'.geo'
demGeocodeFilename = corintFilename+'.demcrop'
geo = self.runGeocode(corintFilename,
widthInt,
corintGeocodeFilename,
demGeocodeFilename)
geoWidth = geo.computeGeoImageWidth()
print("geocodecor: widthGeo = ", geoWidth)
ifgGeocodeFilename = self.insar.geocodeFilename
demCropFilename = self.insar.demCropFilename
topoflatIntFilename = self.insar.topophaseFlatFilename
widthInt = self.insar.resampIntImage.width
geo = self.runGeocode(topoflatIntFilename,
widthInt,
ifgGeocodeFilename,
demCropFilename)
geoWidth = geo.computeGeoImageWidth()
print("geocodeifg: widthGeo = ", geoWidth)
corGeocodeFilename = corFilename + '.geo'
cmplx_to_rmg(corintGeocodeFilename, corGeocodeFilename, geoWidth)
self.geo_to_rsc(ifgGeocodeFilename, corGeocodeFilename)
return None
def geocodeCor4rlks(self):
corFilename = self.insar.coherenceFilename
corintFilename = corFilename.replace('.cor','.corint')
widthInt = self.insar.resampIntImage.width
#Not the right place for this block. Create the 4rlks correlation file and
# geocode it.
#~/Util/isce/components/mroipac/looks/Nbymhgt.py
# topophase.cor topophase_4rlks.cor 9480 4 4
cor4rlksFilename = corFilename.replace('.cor','_4rlks.cor')
from mroipac.looks.Nbymhgt import Nbymhgt
nbymh = Nbymhgt()
nbymh.inputImage = corFilename
nbymh.outputImage = cor4rlksFilename
nbymh.width = widthInt
nbymh.rangeLook = 4
nbymh.azimuthLook = 4
nbymh.nbymhgt()
width4rlksInt = widthInt/4
corint4rlksFilename = cor4rlksFilename.replace('.cor','.corint')
rmg_to_cmplx(cor4rlksFilename, corint4rlksFilename, width4rlksInt)
corint4rlksGeocodeFilename = corint4rlksFilename+'.geo'
dem4rlksGeocodeFilename = corint4rlksFilename+'.demcrop'
geo4rlks = self.runGeocode4rlks(corint4rlksFilename,
width4rlksInt,
corint4rlksGeocodeFilename,
dem4rlksGeocodeFilename)
geo4rlksWidth = geo4rlks.computeGeoImageWidth()
print("geocodecor: widthGeo = ", geo4rlksWidth)
cor4rlksGeocodeFilename = cor4rlksFilename + '.geo'
cmplx_to_rmg(corint4rlksGeocodeFilename, cor4rlksGeocodeFilename, geo4rlksWidth)
# self.geo_to_rsc(ifgGeocodeFilename,corGeocodeFilename)
return None
def geo_to_rsc(self, ifgGeoFile, corGeoFile):
from isceobj.XmlUtil.XmlUtil import XmlUtil
xmlDat = {'Coordinate1':{'size':None,'startingValue':None,'delta':None},
'Coordinate2':{'size':None,'startingValue':None,'delta':None}}
rscXML = {'WIDTH':('Coordinate1','size'),
'X_FIRST':('Coordinate1','startingValue'),
'X_STEP':('Coordinate1','delta'),
'FILE_LENGTH':('Coordinate2','size'),
'Y_FIRST':('Coordinate2','startingValue'),
'Y_STEP':('Coordinate2','delta')}
rscOrder = ('WIDTH','FILE_LENGTH','X_FIRST','X_STEP','Y_FIRST','Y_STEP')
ifgGeoXmlFile = ifgGeoFile + '.xml'
xu = XmlUtil()
xuf = xu.readFile(ifgGeoXmlFile)
c = xuf.findall('component')
for cx in c:
cxn = cx.attrib['name']
p = cx.findall('property')
for e in p:
xmlDat[cxn][e.attrib['name']] = e.findall('value')[0].text
corGeoRscFile = corGeoFile + '.rsc'
with open(corGeoRscFile,'w') as RSC:
spc = " "*25
for a in rscOrder:
RSC.write(
"%s%s%s\n" % (a,spc[0:25-len(a)],xmlDat[rscXML[a][0]][rscXML[a][1]])
)
return None
def runGeocode(self, inFilename, widthIn, geoFilename, demcropFilename):
import stdproc
from isceobj import createDemImage
print("runGeocode: inFilename, widthIn = ", inFilename, widthIn)
print("runGeocode: geoFilename, demcropFilename = ", geoFilename, demcropFilename)
logger.info("Geocoding Image")
# Initialize the Dem
demImage = createDemImage()
IU.copyAttributes(self.insar.demImage, demImage)
demImage.setAccessMode('read')
demImage.createImage()
# Initialize the flattened interferogram
from isceobj import createIntImage, createImage
intImage = createIntImage()
intImage.filename = inFilename
intImage.width = widthIn
intImage.setAccessMode('read')
intImage.createImage()
minLat, maxLat, minLon, maxLon = self.insar.topo.snwe
planet = self.insar.referenceFrame.instrument.getPlatform().getPlanet()
objGeo = stdproc.createGeocode()
objGeo.listInputPorts()
objGeo.wireInputPort(name='peg',object=self.insar.peg)
objGeo.wireInputPort(name='frame',object=self.insar.referenceFrame)
objGeo.wireInputPort(name='planet',object=planet)
objGeo.wireInputPort(name='dem',object=demImage)
objGeo.wireInputPort(name='interferogram',object=intImage)
objGeo.wireInputPort(name='geoPosting', object=self.geoPosting)
print("self.geoPosting = ", self.geoPosting)
objGeo.snwe = minLat, maxLat, minLon, maxLon
objGeo.geoFilename = geoFilename
objGeo.demCropFilename = demcropFilename
#set the tag used in the outfile. each message is precided by this tag
#is the writer is not of "file" type the call has no effect
objGeo.stdWriter = self.stdWriter.set_file_tags("geocode", "log", "err", "out")
# see mocompbaseline
objFormSlc1 = self.insar.formSLC1
mocompPosition1 = objFormSlc1.getMocompPosition()
posIndx = 1
objGeo.referenceOrbit = mocompPosition1[posIndx]
prf1 = self.insar.referenceFrame.instrument.getPulseRepetitionFrequency()
dp = self.insar.dopplerCentroid.getDopplerCoefficients(inHz=False)[0]
v = self.insar.procVelocity
h = self.insar.averageHeight
objGeo.setDopplerCentroidConstantTerm(dp)
objGeo.setBodyFixedVelocity(v)
objGeo.setSpacecraftHeight(h)
objGeo.setNumberRangeLooks(self.insar.numberRangeLooks)
objGeo.setNumberAzimuthLooks(self.insar.numberAzimuthLooks)
# I have no idea what ismocomp means
goodLines = self.insar.numberValidPulses
patchSize = self.insar.patchSize
# this variable was hardcoded in geocode.f90 and was equal to (8192 - 2048)/2
is_mocomp = self.insar.is_mocomp
# is_mocomp = int((patchSize - goodLines)/2)
objGeo.setISMocomp(is_mocomp)
objGeo.geocode()
intImage.finalizeImage()
demImage.finalizeImage()
return objGeo
def runGeocode4rlks(self, inFilename, widthIn, geoFilename, demcropFilename):
import stdproc
from isceobj import createIntImage, createImage
print("runGeocode4rlks: inFilename, widthIn = ", inFilename, widthIn)
print("runGeocode4rlks: geoFilename, demcropFilename = ",
geoFilename,
demcropFilename)
pause(message="Paused in runGeocode4rlks")
logger.info("Geocoding Image")
# Initialize the Dem
from isceobj import createDemImage
demImage = createDemImage()
IU.copyAttributes(self.insar.demImage,demImage)
demImage.setAccessMode('read')
demImage.createImage()
print("demImage.firstLatitude = ", demImage.firstLatitude)
print("demImage.firstLongitude = ", demImage.firstLongitude)
print("demImage.deltaLatitude = ", demImage.deltaLatitude)
print("demImage.deltaLongitude = ", demImage.deltaLongitude)
print("demImage.width = ", demImage.width)
print("demImage.length = ", demImage.length)
demImage_lastLatitude = (
demImage.firstLatitude + (demImage.length-1)*demImage.deltaLatitude
)
demImage_lastLongitude = (
demImage.firstLongitude + (demImage.width-1)*demImage.deltaLongitude
)
print("demImage_lastLatitude = ", demImage_lastLatitude)
print("demImage_lastLongitude = ", demImage_lastLongitude)
# Initialize the input image
intImage = createIntImage()
intImage.setFilename(inFilename)
intImage.setWidth(widthIn)
intImage.setAccessMode('read')
intImage.createImage()
minLat, maxLat, minLon, maxLon = self.insar.topo.snwe
print("objTopo.minLat = ", minLat)
print("objTopo.minLon = ", minLon)
print("objTopo.maxLat = ", maxLat)
print("objTopo.maxLon = ", maxLon)
pause(message="Paused in runGeocode4rlks")
planet = self.insar.referenceFrame.instrument.getPlatform().getPlanet()
objGeo = stdproc.createGeocode()
objGeo.listInputPorts()
objGeo.wireInputPort(name='peg',object=self.insar.peg)
# objGeo.wireInputPort(name='frame',object=self.insar.referenceFrame)
objGeo.rangeFirstSample = self.insar.referenceFrame.getStartingRange()
objGeo.slantRangePixelSpacing = self.insar.referenceFrame.instrument.getRangePixelSize()*4
objGeo.prf = self.insar.referenceFrame.instrument.getPulseRepetitionFrequency()
objGeo.radarWavelength = self.insar.referenceFrame.instrument.getRadarWavelength()
objGeo.wireInputPort(name='planet',object=planet)
objGeo.wireInputPort(name='dem',object=demImage)
objGeo.wireInputPort(name='interferogram',object=intImage)
print("self.geoPosting = ",self.geoPosting)
objGeo.wireInputPort(name='geoPosting',object=self.geoPosting)
objGeo.snwe = minLat, maxLat, minLon, maxLon
objGeo.setGeocodeFilename(geoFilename)
objGeo.setDemCropFilename(demcropFilename)
#set the tag used in the outfile. each message is precided by this tag
#is the writer is not of "file" type the call has no effect
objGeo.stdWriter = self.stdWriter.set_file_tags("geocode", "log", "err", "out")
# see mocompbaseline
objFormSlc1 = self.insar.formSLC1
mocompPosition1 = objFormSlc1.getMocompPosition()
posIndx = 1
objGeo.setReferenceOrbit(mocompPosition1[posIndx])
prf1 = self.insar.referenceFrame.instrument.getPulseRepetitionFrequency()
dp = self.insar.dopplerCentroid.getDopplerCoefficients(inHz=False)[0]
v = self.insar.procVelocity
h = self.insar.averageHeight
objGeo.setDopplerCentroidConstantTerm(dp)
objGeo.setBodyFixedVelocity(v)
objGeo.setSpacecraftHeight(h)
objGeo.setNumberRangeLooks(1.0) #self.insar.numberRangeLooks)
objGeo.setNumberAzimuthLooks(1.0) #self.insar.numberAzimuthLooks)
# I have no idea what ismocomp means
goodLines = self.insar.numberValidPulses
patchSize = self.insar.patchSize
# this variable was hardcoded in geocode.f90 and was equal to (8192 - 2048)/2
is_mocomp = self.insar.is_mocomp
# is_mocomp = int((patchSize - goodLines)/2)
objGeo.setISMocomp(is_mocomp)
objGeo.geocode()
print("Input state paraemters to gecode.f90:")
print("Minimum Latitude = ", objGeo.minimumLatitude)
print("Maximum Latitude = ", objGeo.maximumLatitude)
print("Minimum Longitude = ", objGeo.minimumLongitude)
print("Maximum Longitude = ", objGeo.maximumLongitude)
print("Ellipsoid Major Semi Axis = ", objGeo.ellipsoidMajorSemiAxis)
print("Ellipsoid Eccentricity Squared = ", objGeo.ellipsoidEccentricitySquared)
print("Peg Latitude = ", objGeo.pegLatitude)
print("Peg Longitude = ", objGeo.pegLongitude)
print("Peg Heading = ", objGeo.pegHeading)
print("Range Pixel Spacing = ", objGeo.slantRangePixelSpacing)
print("Range First Sample = ", objGeo.rangeFirstSample)
print("Spacecraft Height = ", objGeo.spacecraftHeight)
print("Planet Local Radius = ", objGeo.planetLocalRadius)
print("Body Fixed Velocity = ", objGeo.bodyFixedVelocity)
print("Doppler Centroid Constant Term = ", objGeo.dopplerCentroidConstantTerm)
print("PRF = ", objGeo.prf)
print("Radar Wavelength = ", objGeo.radarWavelength)
print("S Coordinate First Line = ", objGeo.sCoordinateFirstLine)
print("Azimuth Spacing = ", objGeo.azimuthSpacing)
print("First Latitude = ", objGeo.firstLatitude)
print("First Longitude = ", objGeo.firstLongitude)
print("Delta Latitude = ", objGeo.deltaLatitude)
print("Delta Longitude = ", objGeo.deltaLongitude)
print("Length = ", objGeo.length)
print("Width = ", objGeo.width)
print("Number Range Looks = ", objGeo.numberRangeLooks)
print("Number Azimuth Looks = ", objGeo.numberAzimuthLooks)
print("Number Points Per DEM Post = ", objGeo.numberPointsPerDemPost)
print("Is Mocomp = ", objGeo.isMocomp)
print("DEM Width = ", objGeo.demWidth)
print("DEM Length = ", objGeo.demLength)
# print("Reference Orbit = ", objGeo.referenceOrbit)
print("Dim1 Reference Orbit = ", objGeo.dim1_referenceOrbit)
intImage.finalizeImage()
demImage.finalizeImage()
return objGeo
def runGeocodeCor(self):
import stdproc
logger.info("Geocoding Correlation")
objFormSlc1 = self.insar.formSLC1
# Initialize the Dem
from isceobj import createDemImage, createIntImage, createImage
demImage = createDemImage()
IU.copyAttributes(self.insar.demImage,demImage)
demImage.setAccessMode('read')
demImage.createImage()
topoflatIntFilename = self.insar.topophaseFlatFilename
widthInt = self.insar.resampIntImage.width
intImage = createIntImage()
widthInt = self.insar.resampIntImage.width
intImage.setFilename(corintFilename)
intImage.setWidth(widthInt)
intImage.setAccessMode('read')
intImage.createImage()
posIndx = 1
mocompPosition1 = objFormSlc1.getMocompPosition()
minLat, maxLat, minLon, maxLon = self.insar.topo.snwe
planet = self.insar.referenceFrame.instrument.getPlatform().getPlanet()
objGeo = stdproc.createGeocode()
objGeo.wireInputPort(name='peg',object=self.insar.peg)
objGeo.wireInputPort(name='frame',object=self.insar.referenceFrame)
objGeo.wireInputPort(name='planet',object=planet)
objGeo.wireInputPort(name='dem',object=demImage)
objGeo.wireInputPort(name='interferogram',object=intImage)
objGeo.snwe = minLat, maxLat, minLon, maxLon
corGeocodeFilename = corintFilename+'.geo'
demGeocodeFilename = corintFilename+'.demcrop'
objGeo.setGeocodeFilename(corGeocodeFilename)
objGeo.setDemCropFilename(demGeocodeFilename)
#set the tag used in the outfile. each message is precided by this tag
#is the writer is not of "file" type the call has no effect
objGeo.stdWriter = self.stdWriter.set_file_tags("geocode", "log", "err", "out")
# see mocompbaseline
objGeo.setReferenceOrbit(mocompPosition1[posIndx])
prf1 = self.insar.referenceFrame.instrument.getPulseRepetitionFrequency()
dp = self.insar.dopplerCentroid.getDopplerCoefficients(inHz=False)[0]
v = self.insar.procVelocity
h = self.insar.averageHeight
objGeo.setDopplerCentroidConstantTerm(dp)
objGeo.setBodyFixedVelocity(v)
objGeo.setSpacecraftHeight(h)
objGeo.setNumberRangeLooks(self.insar.numberRangeLooks)
objGeo.setNumberAzimuthLooks(self.insar.numberAzimuthLooks)
# I have no idea what ismocomp means
goodLines = self.insar.numberValidPulses
patchSize = self.insar.patchSize
# this variable was hardcoded in geocode.f90 and was equal to (8192 - 2048)/2
is_mocomp = int((patchSize - goodLines)/2)
objGeo.setISMocomp(is_mocomp)
objGeo.geocode()
intImage.finalizeImage()
demImage.finalizeImage()
return objGeo
def restart(self):
print("Restarting with Filtering")
return
## main() extends _InsarBase.main()
def main(self):
import time
timeStart = time.time()
super(Dpm, self).main()
# self.runCorrect()
self.runShadecpx2rg()
self.runRgoffset()
# Cull offoutliers
self.iterate_runOffoutliers()
self.runResamp_only()
self.insar.topoIntImage=self.insar.resampOnlyImage
self.runTopo()
self.runCorrect()
# Coherence ?
self.runCoherence(method=self.correlation_method)
#ouput the procDoc and pause in order to process coherence off line
#this processing should really be done using _steps.
self.insar.procDoc.renderXml()
pause(message="Paused in main")
# Filter ?
self.runFilter()
# Unwrap ?
self.verifyUnwrap()
# Geocode
self.geocodeCorIfg()
timeEnd = time.time()
logger.info("Total Time: %i seconds" %(timeEnd - timeStart))
self.insar.procDoc.renderXml()
return None
def rmgcor_ifgphs_to_cmplx(rmg,ifg,cpx,width):
import struct
import math
raise DeprecationWarning("Don't ude this function")
length = int(os.stat(ifg).st_size/8./width)
rmgFile = open(rmg,'rb')
ifgFile = open(ifg,'rb')
cpxFile = open(cpx,'wb')
w = int(width)
width2 = 2*w
fmt = "%df" % (width2,)
aCpxLine = [0.0]*width2
for iii in range(length):
anIfgLine = struct.unpack(fmt,ifgFile.read(width2*4))
aRmgLine = struct.unpack(fmt,rmgFile.read(width2*4))
for jjj in range(w):
ifgPhase = math.atan2(anIfgLine[2*jjj+1],anIfgLine[2*jjj])
# ampVal = aRmgLine[jjj]
corVal = aRmgLine[w+jjj]
aCpxLine[2*jjj] = corVal*math.cos(ifgPhase)
aCpxLine[2*jjj+1] = corVal*math.sin(ifgPhase)
cpxFile.write(struct.pack(fmt,*aCpxLine))
rmgFile.close()
ifgFile.close()
cpxFile.close()
return
def ifg1amp_ifg2amp_to_rmg(ifg1,ifg2,rmg,width):
import struct
import math
raise DeprecationWarning("Don't ude this function")
length = int(os.stat(ifg1).st_size/8./width)
ifg1File = open(ifg1,'rb')
ifg2File = open(ifg2,'rb')
rmgFile = open(rmg,'wb')
w = int(width)
width2 = 2*w
fmt = "%df" % (width2,)
aRmgLine = [0.0]*width2
for iii in range(length):
anIfg1Line = struct.unpack(fmt,ifg1File.read(width2*4))
anIfg2Line = struct.unpack(fmt,ifg2File.read(width2*4))
for jjj in range(w):
amp1 = math.sqrt(anIfg1Line[2*jjj]**2 + anIfg1Line[2*jjj+1]**2)
amp2 = math.sqrt(anIfg2Line[2*jjj]**2 + anIfg2Line[2*jjj+1]**2)
aRmgLine[jjj] = amp1
aRmgLine[w + jjj] = amp2
rmgFile.write(struct.pack(fmt,*aRmgLine))
ifg1File.close()
ifg2File.close()
rmgFile.close()
return
def rmg_to_cmplx(rmg,cpx,width):
import struct
import math
length = int(os.stat(rmg).st_size/8./width)
rmgFile = open(rmg,'rb')
cpxFile = open(cpx,'wb')
w = int(width)
width2 = 2*w
fmt = "%df" % (width2,)
aCpxLine = [0.0]*width2
for iii in range(length):
aRmgLine = struct.unpack(fmt,rmgFile.read(width2*4))
for jjj in range(w):
ampVal = aRmgLine[jjj]
corVal = aRmgLine[w+jjj]
aCpxLine[2*jjj] = ampVal
aCpxLine[2*jjj+1] = corVal
cpxFile.write(struct.pack(fmt,*aCpxLine))
rmgFile.close()
cpxFile.close()
return
def cmplx_to_rmg(ifg1,rmg,width):
import struct
import math
length = int(os.stat(ifg1).st_size/8./width)
ifg1File = open(ifg1,'rb')
rmgFile = open(rmg,'wb')
w = int(width)
width2 = 2*w
fmt = "%df" % (width2,)
aRmgLine = [0.0]*width2
for iii in range(length):
anIfg1Line = struct.unpack(fmt,ifg1File.read(width2*4))
for jjj in range(w):
amp1 = anIfg1Line[2*jjj]
amp2 = anIfg1Line[2*jjj+1]
aRmgLine[jjj] = amp1
aRmgLine[w + jjj] = amp2
rmgFile.write(struct.pack(fmt,*aRmgLine))
ifg1File.close()
rmgFile.close()
return
if __name__ == "__main__":
dpm = Dpm()
dpm.run()