diff --git a/README.md b/README.md index cc7b5bb..c57de73 100644 --- a/README.md +++ b/README.md @@ -168,7 +168,7 @@ ISCE application such as 'stripmapApp.py as "> stripmapApp.py" rather than as ### License required for dependencies to enable some workflows in ISCE -Some of the applications, or workflows (such as stripmapApp.py and isceApp.py), +Some of the applications, or workflows (such as insarApp.py and isceApp.py), in ISCE that may be familiar to users will not work with this open source version of ISCE without obtaining licensed components. WinSAR users who have downloaded ISCE from the UNAVCO website (https://winsar.unavco.org/software/isce) have signed diff --git a/components/isceobj/Image/Image.py b/components/isceobj/Image/Image.py index 5de1abd..60b1dfd 100755 --- a/components/isceobj/Image/Image.py +++ b/components/isceobj/Image/Image.py @@ -804,7 +804,7 @@ class ImageCoordinate(Configurable): doc="Starting value of the coordinate.") self._coordEnd = self.parameter('coordEnd', public_name='endingValue', default=None, units='', type=float, mandatory=False, - doc="Starting value of the coordinate.") + doc="Ending value of the coordinate.") self._coordDelta = self.parameter('coordDelta', public_name='delta', default=1, units='', type=float, mandatory=False, doc="Coordinate quantization.") diff --git a/components/isceobj/Sensor/SConscript b/components/isceobj/Sensor/SConscript index e44efe4..c5afbc3 100644 --- a/components/isceobj/Sensor/SConscript +++ b/components/isceobj/Sensor/SConscript @@ -25,7 +25,7 @@ listFiles = ['ALOS.py','CEOS.py','COSMO_SkyMed.py','COSMO_SkyMed_SLC.py', 'Radarsat2.py','TerraSARX.py','Polarimetry.py','Sensor.py', 'ROI_PAC.py','Sentinel1.py','TanDEMX.py','KOMPSAT5.py', 'Risat1.py', 'Risat1_SLC.py', 'UAVSAR_RPI.py', 'UAVSAR_Stack.py', 'UAVSAR_Polsar.py', - 'ALOS2.py', 'ERS_SLC.py', 'ALOS_SLC.py', 'EnviSAT_SLC.py', 'ERS_EnviSAT_SLC.py', '__init__.py'] + 'ALOS2.py', 'ERS_SLC.py', 'ALOS_SLC.py', 'EnviSAT_SLC.py', 'ERS_EnviSAT_SLC.py', 'SICD_RGZERO.py','__init__.py'] helpList,installHelp = envSensor['HELP_BUILDER'](envSensor,'__init__.py',install) envSensor.Install(installHelp,helpList) diff --git a/components/isceobj/Sensor/SICD_RGZERO.py b/components/isceobj/Sensor/SICD_RGZERO.py new file mode 100644 index 0000000..88460f6 --- /dev/null +++ b/components/isceobj/Sensor/SICD_RGZERO.py @@ -0,0 +1,247 @@ +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Piyush Agram +# Copyright 2010, by the California Institute of Technology. ALL RIGHTS +# RESERVED. United States Government Sponsorship acknowledged. Any commercial +# use must be negotiated with the Office of Technology Transfer at the +# California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this +# software, the user agrees to comply with all applicable U.S. export laws and +# regulations. User has the responsibility to obtain export licenses, or other +# export authority as may be required before exporting such information to +# foreign countries or providing access to foreign persons. +# +# NASA Jet Propulsion Laboratory +# California Institute of Technology +# (C) 2010 All Rights Reserved +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +import datetime +import isceobj +from isceobj.Orbit.Orbit import StateVector +from isceobj.Planet.AstronomicalHandbook import Const +from isceobj.Planet.Planet import Planet +from isceobj.Scene.Frame import Frame +from iscesys.DateTimeUtil.DateTimeUtil import DateTimeUtil as DTUtil +from iscesys.Component.Component import Component +from isceobj.Sensor.Sensor import Sensor + + +SICD = Component.Parameter( + 'sicd', + public_name='SICD', + default=None, + type=str, + mandatory=True, + intent='input', + doc='SICD input file') + +class SICD_RGZERO(Sensor): + """ + A class to parse SICD RGZERO metadata + """ + + parameter_list = (SICD,) + Sensor.parameter_list + logging_name = "isce.sensor.SICD_RGZERO" + family_name = "sicd_rgzero" + + def __init__(self): + super(SICD_RGZERO,self).__init__() + self._sicdmeta = None + + return None + + def getFrame(self): + return self.frame + + + def parse(self): + try: + import sarpy.io.complex as cf + except ImportError: + raise Exception('You need to install sarpy from NGA - https://github.com/ngageoint/sarpy to work with SICD data') + self._sicdmeta = cf.open(self.sicd).sicdmeta + self.populateMetadata() + + + def _populatePlatform(self): + mdict = self._sicdmeta + platform = self.frame.getInstrument().getPlatform() + + platform.setMission(mdict.CollectionInfo.CollectorName) + platform.setPlanet(Planet(pname="Earth")) + side = mdict.SCPCOA.SideOfTrack + if side.startswith('R'): + side = -1 + else: + side = 1 + platform.setPointingDirection(side) + + if mdict.CollectionInfo.RadarMode.ModeType.upper() != 'STRIPMAP': + raise Exception('SICD ModeType should be STRIPMAP') + + if mdict.CollectionInfo.CollectType.upper() != 'MONOSTATIC': + raise Exception('SICD ModeType should be MONOSTATIC') + + + def _populateInstrument(self, mdict=None): + if mdict is None: + mdict = self._sicdmeta + + instrument = self.frame.getInstrument() + + ###Ensure that data is actually SICD RGZERO + if (mdict.Grid.Type != 'RGZERO'): + raise Exception('Input data must be SICD RGZERO') + + if (mdict.Grid.ImagePlane != 'SLANT'): + raise Exception('Input data must be SICD RGZERO in Slant Range plane') + + rangePixelSize = mdict.Grid.Row.SS + azimuthPixelSize = mdict.Grid.Col.SS + fs = Const.c/(2*rangePixelSize) + + fc = mdict.RMA.INCA.FreqZero + prf = mdict.Timeline.IPP.Set.IPPPoly[1] * mdict.ImageFormation.RcvChanProc.PRFScaleFactor + + instrument.setRadarWavelength(Const.c/fc) + instrument.setPulseRepetitionFrequency(prf) + instrument.setRangePixelSize(rangePixelSize) + instrument.setPulseLength(mdict.RadarCollection.Waveform.WFParameters[0].TxPulseLength) + instrument.setChirpSlope(mdict.RadarCollection.Waveform.WFParameters[0].TxRFBandwidth / mdict.RadarCollection.Waveform.WFParameters[0].TxPulseLength ) + instrument.setRangeSamplingRate(fs) + instrument.setInPhaseValue(0.) + instrument.setQuadratureValue(0.) + instrument.platform.setAntennaLength(2.2 * azimuthPixelSize) + + def _populateFrame(self, mdict=None): + if mdict is None: + mdict = self._sicdmeta + + startRange = mdict.RMA.INCA.R_CA_SCP - (mdict.ImageData.SCPPixel.Row * mdict.Grid.Row.SS) + + ####Compute the UTC times + zd_t_scp = mdict.RMA.INCA.TimeCAPoly[0] + ss_zd_s = 1 /self.frame.PRF + sensingStart = mdict.Timeline.CollectStart + datetime.timedelta(seconds = (zd_t_scp - mdict.ImageData.SCPPixel.Col * ss_zd_s)) + sensingStop = sensingStart + datetime.timedelta(seconds = (mdict.ImageData.NumCols-1) / self.frame.PRF) + sensingMid = sensingStart + 0.5 * (sensingStop - sensingStart) + + self.frame.setStartingRange(startRange) + if mdict.SCPCOA.ARPVel.Z > 0: + self.frame.setPassDirection('ASCENDING') + else: + self.frame.setPassDirection('DESCENDING') + + self.frame.setOrbitNumber(9999) + self.frame.setProcessingFacility(mdict.ImageCreation.Site) + self.frame.setProcessingSoftwareVersion(mdict.ImageCreation.Application) + + pol = mdict.ImageFormation.TxRcvPolarizationProc + self.frame.setPolarization(pol[0] + pol[2]) + self.frame.setNumberOfLines(mdict.ImageData.NumCols) + self.frame.setNumberOfSamples(mdict.ImageData.NumRows) + + + self.frame.setSensingStart(sensingStart) + self.frame.setSensingMid(sensingMid) + self.frame.setSensingStop(sensingStop) + + rangePixelSize = self.frame.getInstrument().getRangePixelSize() + farRange = startRange + self.frame.getNumberOfSamples()*rangePixelSize + self.frame.setFarRange(farRange) + + + def _populateOrbit(self, mdict=None): + import numpy.polynomial.polynomial as poly + if mdict is None: + mdict = self._sicdmeta + + raw_start_time = mdict.Timeline.CollectStart + + tmin = self.frame.sensingStart - datetime.timedelta(seconds=5) + tmax = self.frame.sensingStop + datetime.timedelta(seconds=5) + + + orbit = self.frame.getOrbit() + orbit.setReferenceFrame('ECEF') + orbit.setOrbitSource('Header') + + posX = mdict.Position.ARPPoly.X + posY = mdict.Position.ARPPoly.Y + posZ = mdict.Position.ARPPoly.Z + velX = poly.polyder(posX) + velY = poly.polyder(posY) + velZ = poly.polyder(posZ) + + tinp = tmin + while tinp <= tmax: + + deltaT = (tinp - raw_start_time).total_seconds() + vec = StateVector() + vec.setTime(tinp) + vec.setPosition([poly.polyval(deltaT, posX), + poly.polyval(deltaT, posY), + poly.polyval(deltaT, posZ)]) + vec.setVelocity([poly.polyval(deltaT, velX), + poly.polyval(deltaT, velY), + poly.polyval(deltaT, velZ)]) + + orbit.addStateVector(vec) + tinp = tinp + datetime.timedelta(seconds=1) + + + def populateImage(self): + import sarpy.io.complex as cf + + img = cf.open(self.sicd) + data = img.read_chip() + data.T.tofile(self.output) + + rawImage = isceobj.createSlcImage() + rawImage.setByteOrder('l') + rawImage.setFilename(self.output) + rawImage.setAccessMode('read') + rawImage.setWidth(self.frame.getNumberOfSamples()) + rawImage.setXmax(self.frame.getNumberOfSamples()) + rawImage.setXmin(0) + self.getFrame().setImage(rawImage) + #rawImage.renderHdr() + + def _populateExtras(self): + """ + Populate some extra fields. + """ + from sarpy.geometry.point_projection import coa_projection_set + import numpy as np + mdict = self._sicdmeta + + ###Imagesize + rows = np.linspace(0., mdict.ImageData.NumRows*1.0, num=3) + rdot = [] + + for grow in rows: + pt = coa_projection_set(mdict,[grow,0]) + rdot.append( pt[1][0]) + + self.frame._dopplerVsPixel = list(np.polyfit(rows, rdot, 2)[::-1]) + + + def extractImage(self): + """Extract the raw image data""" + self.parse() + self._populateExtras() + self.populateImage() + + def extractDoppler(self): + """ + Return the doppler centroid as defined in the HDF5 file. + """ + dopp = self.frame._dopplerVsPixel + quadratic = {} + quadratic['a'] = dopp[0] + quadratic['b'] = dopp[1] + quadratic['c'] = dopp[2] + return quadratic + diff --git a/components/isceobj/Sensor/__init__.py b/components/isceobj/Sensor/__init__.py index 39d5eda..39ac127 100755 --- a/components/isceobj/Sensor/__init__.py +++ b/components/isceobj/Sensor/__init__.py @@ -95,6 +95,7 @@ createALOS_SLC = partial(factory_template, 'ALOS_SLC') createEnviSAT_SLC = partial(factory_template, 'EnviSAT_SLC') createERS_ENVISAT = partial(factory_template, 'ERS_ENVISAT') createERS_EnviSAT_SLC = partial(factory_template, 'ERS_EnviSAT_SLC') +createSICD_RGZERO = partial(factory_template, 'SICD_RGZERO') SENSORS = {'ALOS' : createALOS, 'ALOS_SLC' : createALOS_SLC, @@ -117,7 +118,8 @@ SENSORS = {'ALOS' : createALOS, 'SENTINEL1' : createSentinel1, 'ENVISAT_SLC': createEnviSAT_SLC, 'ERS_ENVISAT' : createERS_ENVISAT, - 'ERS_ENVISAT_SLC' : createERS_EnviSAT_SLC} + 'ERS_ENVISAT_SLC' : createERS_EnviSAT_SLC, + 'SICD_RGZERO' : createSICD_RGZERO} #These are experimental and can be added in as they become ready # 'JERS': createJERS, diff --git a/components/isceobj/StripmapProc/Factories.py b/components/isceobj/StripmapProc/Factories.py index be7beb1..9f65fc1 100644 --- a/components/isceobj/StripmapProc/Factories.py +++ b/components/isceobj/StripmapProc/Factories.py @@ -52,7 +52,7 @@ def isRawSensor(sensor): ''' Check if input data is raw / slc. ''' - if str(sensor).lower() in ["terrasarx","cosmo_skymed_slc","radarsat2",'tandemx', 'kompsat5','risat1_slc','sentinel1', 'alos2','ers_slc','alos_slc','envisat_slc', 'uavsar_rpi','ers_envisat_slc']: + if str(sensor).lower() in ["terrasarx","cosmo_skymed_slc","radarsat2",'tandemx', 'kompsat5','risat1_slc','sentinel1', 'alos2','ers_slc','alos_slc','envisat_slc', 'uavsar_rpi','ers_envisat_slc','sicd_rgzero']: return False else: return True @@ -63,7 +63,7 @@ def isZeroDopplerSLC(sensor): Check if SLC is zero doppler / native doppler. ''' - if str(sensor).lower() in ["terrasarx","cosmo_skymed_slc","radarsat2",'tandemx', 'kompsat5','risat1_slc','sentinel1', 'alos2','ers_slc','envisat_slc','ers_envisat_slc']: + if str(sensor).lower() in ["terrasarx","cosmo_skymed_slc","radarsat2",'tandemx', 'kompsat5','risat1_slc','sentinel1', 'alos2','ers_slc','envisat_slc','ers_envisat_slc','sicd_rgzero']: return True elif sensor.lower() in ['alos_slc', 'uavsar_rpi']: return False @@ -76,7 +76,7 @@ def getDopplerMethod(sensor): Return appropriate doppler method based on user input. ''' - if str(sensor).lower() in ["terrasarx","cosmo_skymed_slc","radarsat2",'tandemx', 'kompsat5','risat1_slc','sentinel1', 'alos2','ers_slc','alos_slc','envisat_slc', 'uavsar_rpi','cosmo_skymed','ers_envisat_slc']: + if str(sensor).lower() in ["terrasarx","cosmo_skymed_slc","radarsat2",'tandemx', 'kompsat5','risat1_slc','sentinel1', 'alos2','ers_slc','alos_slc','envisat_slc', 'uavsar_rpi','cosmo_skymed','ers_envisat_slc','sicd_rgzero']: res = 'useDEFAULT' else: res = 'useDOPIQ' diff --git a/components/isceobj/StripmapProc/runDenseOffsets.py b/components/isceobj/StripmapProc/runDenseOffsets.py index 4887b85..eb27366 100644 --- a/components/isceobj/StripmapProc/runDenseOffsets.py +++ b/components/isceobj/StripmapProc/runDenseOffsets.py @@ -65,7 +65,7 @@ def estimateOffsetField(master, slave, denseOffsetFileName, objOffset.offsetImageName = denseOffsetFileName + '.bil' objOffset.snrImageName = denseOffsetFileName +'_snr.bil' - + objOffset.covImageName = denseOffsetFileName +'_cov.bil' objOffset.denseampcor(sar, sim) diff --git a/components/isceobj/TopsProc/TopsProc.py b/components/isceobj/TopsProc/TopsProc.py index 5f7778d..9e250e8 100644 --- a/components/isceobj/TopsProc/TopsProc.py +++ b/components/isceobj/TopsProc/TopsProc.py @@ -320,6 +320,14 @@ OFFSET_SNR_FILE = Component.Parameter( mandatory=False, doc='Filename for gross dense offsets SNR. Used in runDenseOffsets.') +OFFSET_COV_FILE = Component.Parameter( + 'covfile', + public_name='Offset covariance filename', + default='dense_offsets_cov.bil', + type=str, + mandatory=False, + doc='Filename for gross dense offsets covariance. Used in runDenseOffsets.') + FILT_OFFSET_OUTPUT_FILE = Component.Parameter( 'filt_offsetfile', public_name='Filtered offset filename', @@ -333,6 +341,7 @@ OFFSET_GEOCODE_LIST = Component.Parameter('off_geocode_list', public_name='offset geocode list', default = [OFFSET_OUTPUT_FILE, OFFSET_SNR_FILE, + OFFSET_COV_FILE, FILT_OFFSET_OUTPUT_FILE], container = list, type=str, @@ -385,6 +394,7 @@ class TopsProc(Component): OFFSET_WIDTH, OFFSET_OUTPUT_FILE, OFFSET_SNR_FILE, + OFFSET_COV_FILE, FILT_OFFSET_OUTPUT_FILE, OFFSET_GEOCODE_LIST) diff --git a/components/isceobj/TopsProc/runDenseOffsets.py b/components/isceobj/TopsProc/runDenseOffsets.py index ff288d3..f04b571 100644 --- a/components/isceobj/TopsProc/runDenseOffsets.py +++ b/components/isceobj/TopsProc/runDenseOffsets.py @@ -118,9 +118,11 @@ def runDenseOffsetsCPU(self): objOffset.offsetImageName = os.path.join(self._insar.mergedDirname, self._insar.offsetfile) objOffset.snrImageName = os.path.join(self._insar.mergedDirname, self._insar.snrfile) + objOffset.covImageName = os.path.join(self._insar.mergedDirname, self._insar.covfile) print('Output dense offsets file name: %s' % (objOffset.offsetImageName)) print('Output SNR file name: %s' % (objOffset.snrImageName)) + print('Output covariance file name: %s' % (objOffset.covImageName)) print('\n======================================') print('Running dense ampcor...') print('======================================\n') diff --git a/components/mroipac/ampcor/DenseAmpcor.py b/components/mroipac/ampcor/DenseAmpcor.py index f11c2d5..eb4202f 100755 --- a/components/mroipac/ampcor/DenseAmpcor.py +++ b/components/mroipac/ampcor/DenseAmpcor.py @@ -253,6 +253,13 @@ SNR_IMAGE_NAME = Component.Parameter('snrImageName', mandatory=False, doc = 'File name for output SNR') +COV_IMAGE_NAME = Component.Parameter('covImageName', + public_name = 'COV_IMAGE_NAME', + default = 'dense_ampcor_cov.bil', + type=str, + mandatory=False, + doc = 'File name for output covariance') + MARGIN = Component.Parameter('margin', public_name = 'MARGIN', default = 50, @@ -298,6 +305,7 @@ class DenseAmpcor(Component): BAND2, OFFSET_IMAGE_NAME, SNR_IMAGE_NAME, + COV_IMAGE_NAME, MARGIN, NUMBER_THREADS) @@ -374,10 +382,16 @@ class DenseAmpcor(Component): self.locationAcross = np.frombuffer(mp.Array('i', numlen).get_obj(), dtype='i') self.locationAcrossOffset = np.frombuffer(mp.Array('f', numlen).get_obj(), dtype='f') self.snr = np.frombuffer(mp.Array('f', numlen).get_obj(), dtype='f') + self.cov1 = np.frombuffer(mp.Array('f', numlen).get_obj(), dtype='f') + self.cov2 = np.frombuffer(mp.Array('f', numlen).get_obj(), dtype='f') + self.cov3 = np.frombuffer(mp.Array('f', numlen).get_obj(), dtype='f') self.locationDownOffset[:] = -10000.0 self.locationAcrossOffset[:] = -10000.0 self.snr[:] = 0.0 + self.cov1[:] = 999.0 + self.cov2[:] = 999.0 + self.cov3[:] = 999.0 ###run ampcor on parallel processes threads = [] @@ -421,7 +435,6 @@ class DenseAmpcor(Component): #### Scale images (default is 1.0 to keep as pixel) self.locationDownOffset *= self.scaling_factor self.locationAcrossOffset *= self.scaling_factor - self.snr *= self.scaling_factor self.write_slantrange_images() @@ -487,6 +500,9 @@ class DenseAmpcor(Component): self.locationAcross[firstind+i] = objAmpcor.locationAcross[j] self.locationAcrossOffset[firstind+i] = objAmpcor.locationAcrossOffset[j] self.snr[firstind+i] = objAmpcor.snrRet[j] + self.cov1[firstind+i] = objAmpcor.cov1Ret[j] + self.cov2[firstind+i] = objAmpcor.cov2Ret[j] + self.cov3[firstind+i] = objAmpcor.cov3Ret[j] j += 1 else: self.locationDown[firstind+i] = downInd @@ -494,6 +510,9 @@ class DenseAmpcor(Component): self.locationAcross[firstind+i] = acInd self.locationAcrossOffset[firstind+i] = -10000. self.snr[firstind+i] = 0. + self.cov1[firstind+i] = 999. + self.cov2[firstind+i] = 999. + self.cov3[firstind+i] = 999. return @@ -518,6 +537,14 @@ class DenseAmpcor(Component): if self.locationAcross.ndim == 1: self.locationAcross = self.locationAcross.reshape(-1,self.offsetCols) + if self.cov1.ndim == 1: + self.cov1 = self.cov1.reshape(-1,self.offsetCols) + + if self.cov2.ndim == 1: + self.cov2 = self.cov2.reshape(-1,self.offsetCols) + + if self.cov3.ndim == 1: + self.cov3 = self.cov3.reshape(-1,self.offsetCols) outdata = np.empty((2*self.offsetLines, self.offsetCols), dtype=np.float32) outdata[::2,:] = self.locationDownOffset @@ -544,7 +571,23 @@ class DenseAmpcor(Component): snrImg.setLength(self.offsetLines) snrImg.setAccessMode('read') snrImg.renderHdr() - + + ####Create covariance image + covdata = np.empty((3*self.offsetLines, self.offsetCols), dtype=np.float32) + covdata[::3,:] = self.cov1 + covdata[1::3,:] = self.cov2 + covdata[2::3,:] = self.cov3 + covdata.tofile(self.covImageName) + del covdata + covImg = isceobj.createImage() + covImg.setDataType('FLOAT') + covImg.setFilename(self.covImageName) + covImg.setBands(3) + covImg.scheme = 'BIL' + covImg.setWidth(self.offsetCols) + covImg.setLength(self.offsetLines) + covImg.setAccessMode('read') + covImg.renderHdr() def checkTypes(self): '''Check if the image datatypes are set.''' diff --git a/examples/input_files/master_SICD_RGZERO.xml b/examples/input_files/master_SICD_RGZERO.xml new file mode 100644 index 0000000..acc5951 --- /dev/null +++ b/examples/input_files/master_SICD_RGZERO.xml @@ -0,0 +1,19 @@ + + + data/RS2_RS2_OK64213_PK591209_DK527307_U14_20150612_151051_HH_SLC.nitf + + + 20110311 + + + +