#!/usr/bin/env python3 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Copyright 2014 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: Piyush Agram #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ import isce import xml.etree.ElementTree as ElementTree from collections import defaultdict import datetime import isceobj from isceobj.Util import Poly1D, Poly2D from isceobj.Planet.Planet import Planet from isceobj.Orbit.Orbit import StateVector, Orbit from isceobj.Orbit.OrbitExtender import OrbitExtender from isceobj.Planet.AstronomicalHandbook import Const from iscesys.Component.Component import Component from iscesys.DateTimeUtil.DateTimeUtil import DateTimeUtil as DTUtil from .GRDProduct import GRDProduct import os import glob import json import numpy as np import shelve SAFE = Component.Parameter('safe', public_name = 'safe', default = None, type = str, doc = 'SAFE folder with the GRD product') POLARIZATION = Component.Parameter('polarization', public_name = 'polarization', default = None, type = str, doc = 'Polarization to unpack') ORBIT_DIR = Component.Parameter('orbitDir', public_name = 'orbit directory', default = None, type = str, doc = 'Directory to search for orbit data') ORBIT_FILE = Component.Parameter('orbitFile', public_name = 'orbit file', default = None, type = str, doc = 'External orbit file with state vectors') OUTPUT = Component.Parameter('output', public_name = 'output directory', default = None, type = str, doc = 'Directory where the data gets unpacked') SLANT_RANGE_FILE = Component.Parameter('slantRangeFile', public_name = 'slant range file', default = None, type = str, doc = 'Slant range file at full resolution') ####List of facilities PRODUCT = Component.Facility('product', public_name='product', module = 'isceobj.Sensor.GRD', factory = 'createGRDProduct', args = (), mandatory = True, doc = 'GRD Product populated by the reader') class Sentinel1(Component): """ A Class representing Sentinel1 data """ family = 's1grd' logging = 'isce.sensor.grd.s1' parameter_list = (SAFE, POLARIZATION, ORBIT_DIR, ORBIT_FILE, OUTPUT, SLANT_RANGE_FILE) facility_list = (PRODUCT,) def __init__(self): Component.__init__(self) self.xml = None self.tiff = None self.calibrationXml = None self.noiseXml = None self.manifest = None self.noiseCorrectionApplied = False self.betaLUT = None self.gr2srLUT = None self.noiseRangeLUT = None self.noiseAzimuthLUT = None self._xml_root=None def validateUserInputs(self): ''' Validate inputs provided by user. Populat tiff and xml list using SAFE folder names. ''' import zipfile import fnmatch if self.safe is None: raise Exception('SAFE directory is not provided') if self.polarization in ['',None]: raise Exception('Polarization is not provided') ###Check if zip file / unpacked directory is provided. iszipfile = False if self.safe.endswith('.zip'): iszipfile = True with zipfile.ZipFile(self.safe, 'r') as zf: namelist = zf.namelist() ####First find VV annotation file swathid = 's1?-iw-grd-' + self.polarization.lower() if iszipfile: #Find XML file patt = os.path.join('*.SAFE', 'annotation', swathid+'*') match = fnmatch.filter(namelist, patt) if len(match) == 0 : raise Exception('Annotation file for {0} not found in {1}'.format(self.polarization, self.safe)) self.xml = os.path.join('/vsizip', self.safe, match[0]) #Find TIFF file patt = os.path.join('*.SAFE', 'measurement', swathid+'*') match = fnmatch.filter(namelist, patt) if len(match) == 0 : raise Exception('Annotation file found for {0} but no measurement in {1}'.format(self.polarization, self.safe)) self.tiff = os.path.join('/vsizip', self.safe, match[0]) #Find Calibration file patt = os.path.join('*.SAFE', 'annotation', 'calibration', 'calibration-'+swathid+'*') match = fnmatch.filter(namelist, patt) if len(match) == 0 : raise Exception('Annotation file found for {0} but no calibration in {1}'.format(self.polarization, self.safe)) self.calibrationXml = os.path.join('/vsizip', self.safe, match[0]) #Find Noise file patt = os.path.join('*.SAFE', 'annotation', 'calibration', 'noise-'+swathid+'*') match = fnmatch.filter(namelist, patt) if len(match) == 0 : raise Exception('Annotation file found for {0} but no noise in {1}'.format(self.polarization, self.safe)) self.noiseXml = os.path.join('/vsizip', self.safe, match[0]) patt = os.path.join('*.SAFE', 'manifest.safe') match = fnmatch.filter(namelist, patt) if len(match) == 0: raise Exception('No manifest file found in {0}'.format(self.safe)) self.manifest = os.path.join('/vsizip', self.safe, match[0]) else: ####Find annotation file patt = os.path.join( self.safe, 'annotation', swathid + '*') match = glob.glob(patt) if len(match) == 0: raise Exception('Annotation file for {0} not found in {1}'.format(self.polarization, self.safe)) self.xml = match[0] ####Find TIFF file patt = os.path.join( self.safe, 'measurement', swathid+'*') match = glob.glob(patt) if len(match) == 0: raise Exception('Annotation file found for {0} but not measurement in {1}'.format(self.polarization, self.safe)) self.tiff= match[0] ####Find calibration file patt = os.path.join( self.safe, 'annotation', 'calibration', 'calibration-' + swathid + '*') match = glob.glob(patt) if len(match) == 0 : raise Exception('Annotation file found for {0} but not calibration in {1}'.format(self.polarization, self.safe)) self.calibrationXml= match[0] ####Find noise file patt = os.path.join( self.safe, 'annotation', 'calibration', 'noise-' + swathid + '*') match = glob.glob(patt) if len(match) == 0 : raise Exception('Annotation file found for {0} but not noise in {1}'.format(self.polarization, self.safe)) self.noiseXml = match[0] ####Find manifest file self.manifest = os.path.join(self.safe, 'manifest.safe') print('XML: ', self.xml) print('TIFF: ', self.tiff) print('CALIB: ', self.calibrationXml) print('NOISE: ', self.noiseXml) print('MANIFEST: ', self.manifest) return def parse(self): import zipfile self.validateUserInputs() if '.zip' in self.xml: try: parts = self.xml.split(os.path.sep) zipname = os.path.join('/',*(parts[:-3])) fname = os.path.join(*(parts[-3:])) with zipfile.ZipFile(zipname, 'r') as zf: xmlstr = zf.read(fname) except: raise Exception('Could not read xml file {0}'.format(self.xml)) else: try: with open(self.xml, 'r') as fid: xmlstr = fid.read() except: raise Exception('Could not read xml file {0}'.format(self.xml)) self._xml_root = ElementTree.fromstring(xmlstr) self.populateMetadata() self.populateBbox() ####Try and locate an orbit file if self.orbitFile is None: if self.orbitDir is not None: self.orbitFile = self.findOrbitFile() print('Found this orbitfile: %s' %self.orbitFile) ####Read in the orbits if '_POEORB_' in self.orbitFile: orb = self.extractPreciseOrbit() elif '_RESORB_' in self.orbitFile: orb = self.extractOrbit() self.product.orbit.setOrbitSource('Header') for sv in orb: self.product.orbit.addStateVector(sv) self.populateIPFVersion() self.extractBetaLUT() self.extractNoiseLUT() def getxmlattr(self, path, key): try: res = self._xml_root.find(path).attrib[key] except: raise Exception('Cannot find attribute %s at %s'%(key, path)) return res def getxmlvalue(self, path): try: res = self._xml_root.find(path).text except: raise Exception('Tag= %s not found'%(path)) if res is None: raise Exception('Tag = %s not found'%(path)) return res def getxmlelement(self, path): try: res = self._xml_root.find(path) except: raise Exception('Cannot find path %s'%(path)) if res is None: raise Exception('Cannot find path %s'%(path)) return res def convertToDateTime(self, string): dt = datetime.datetime.strptime(string,"%Y-%m-%dT%H:%M:%S.%f") return dt def populateMetadata(self): """ Create metadata objects from the metadata files """ ####Set each parameter one - by - one mission = self.getxmlvalue('adsHeader/missionId') swath = self.getxmlvalue('adsHeader/swath') polarization = self.getxmlvalue('adsHeader/polarisation') orbitnumber = int(self.getxmlvalue('adsHeader/absoluteOrbitNumber')) frequency = float(self.getxmlvalue('generalAnnotation/productInformation/radarFrequency')) passDirection = self.getxmlvalue('generalAnnotation/productInformation/pass') groundRangePixelSize = float(self.getxmlvalue('imageAnnotation/imageInformation/rangePixelSpacing')) azimuthPixelSize = float(self.getxmlvalue('imageAnnotation/imageInformation/azimuthPixelSpacing')) azimuthTimeInterval = float(self.getxmlvalue('imageAnnotation/imageInformation/azimuthTimeInterval')) lines = int(self.getxmlvalue('imageAnnotation/imageInformation/numberOfLines')) samples = int(self.getxmlvalue('imageAnnotation/imageInformation/numberOfSamples')) slantRangeTime = float(self.getxmlvalue('imageAnnotation/imageInformation/slantRangeTime')) startingSlantRange = float(self.getxmlvalue('imageAnnotation/imageInformation/slantRangeTime'))*Const.c/2.0 incidenceAngle = float(self.getxmlvalue('imageAnnotation/imageInformation/incidenceAngleMidSwath')) sensingStart = self.convertToDateTime( self.getxmlvalue('imageAnnotation/imageInformation/productFirstLineUtcTime')) sensingStop = self.convertToDateTime( self.getxmlvalue('imageAnnotation/imageInformation/productLastLineUtcTime')) ####Sentinel is always right looking lookSide = -1 ###Read ascending node for phase calibration ascTime = self.convertToDateTime(self.getxmlvalue('imageAnnotation/imageInformation/ascendingNodeTime')) ###Noise correction correctionApplied = self.getxmlvalue('imageAnnotation/processingInformation/thermalNoiseCorrectionPerformed').upper() == 'TRUE' self.product.lookSide = 'RIGHT' self.product.numberOfSamples = samples self.product.numberOfLines = lines self.product.startingGroundRange = 0.0 self.product.startingSlantRange = startingSlantRange self.product.trackNumber = ((orbitnumber-73)%175) + 1 self.product.orbitNumber = orbitnumber self.product.frameNumber = 1 self.product.polarization = polarization self.product.passDirection = passDirection self.product.radarWavelength = Const.c / frequency self.product.groundRangePixelSize = groundRangePixelSize self.product.azimuthPixelSize = azimuthPixelSize self.product.azimuthTimeInterval = azimuthTimeInterval self.product.ascendingNodeTime = ascTime self.product.slantRangeTime = slantRangeTime self.product.sensingStart = sensingStart self.product.sensingStop = sensingStop self.noiseCorrectionApplied = correctionApplied def populateBbox(self, margin=0.1): ''' Populate the bounding box from metadata. ''' glist = self.getxmlelement('geolocationGrid/geolocationGridPointList') lat = [] lon = [] for child in glist: lat.append( float(child.find('latitude').text)) lon.append( float(child.find('longitude').text)) self.product.bbox = [min(lat) - margin, max(lat) + margin, min(lon) - margin, max(lon) + margin] print(self.product.bbox) return def populateIPFVersion(self): ''' Get IPF version from the manifest file. ''' if self.manifest is None: return nsp = "{http://www.esa.int/safe/sentinel-1.0}" if '.zip' in self.manifest: import zipfile parts = self.manifest.split(os.path.sep) zipname = os.path.join('/',*(parts[:-2])) fname = os.path.join(*(parts[-2:])) try: with zipfile.ZipFile(zipname, 'r') as zf: xmlstr = zf.read(fname) except: raise Exception('Could not read manifest file : {0}'.format(self.manifest)) else: try: with open(self.manifest, 'r') as fid: xmlstr = fid.read() except: raise Exception('Could not read manifest file: {0}'.format(self.manifest)) try: root = ElementTree.fromstring(xmlstr) elem = root.find('.//metadataObject[@ID="processing"]') rdict = elem.find('.//xmlData/' + nsp + 'processing/' + nsp + 'facility/' + nsp + 'software').attrib self.IPFversion = rdict['version'] print('Setting IPF version to : ', self.IPFversion) except: print('Could not read version number successfully from manifest file: ', self.manifest) pass return def findOrbitFile(self): ''' Find correct orbit file in the orbit directory. ''' datefmt = "%Y%m%dT%H%M%S" types = ['POEORB', 'RESORB'] filelist = [] match = [] timeStamp = self.product.sensingStart+(self.product.sensingStop - self.product.sensingStart)/2. for orbType in types: files = glob.glob( os.path.join(self.orbitDir, 'S1?_OPER_AUX_' + orbType + '_OPOD*')) filelist.extend(files) ###List all orbit files for result in filelist: fields = result.split('_') taft = datetime.datetime.strptime(fields[-1][0:15], datefmt) tbef = datetime.datetime.strptime(fields[-2][1:16], datefmt) print(taft, tbef) #####Get all files that span the acquisition if (tbef <= timeStamp) and (taft >= timeStamp): tmid = tbef + 0.5 * (taft - tbef) match.append((result, abs((timeStamp-tmid).total_seconds()))) #####Return the file with the image is aligned best to the middle of the file if len(match) != 0: bestmatch = min(match, key = lambda x: x[1]) return bestmatch[0] if len(match) == 0: raise Exception('No suitable orbit file found. If you want to process anyway - unset the orbitdir parameter') def extractOrbit(self): ''' Extract orbit information from xml node. ''' node = self._xml_root.find('generalAnnotation/orbitList') print('Extracting orbit from annotation XML file') frameOrbit = Orbit() frameOrbit.configure() for child in node: timestamp = self.convertToDateTime(child.find('time').text) pos = [] vel = [] posnode = child.find('position') velnode = child.find('velocity') for tag in ['x','y','z']: pos.append(float(posnode.find(tag).text)) for tag in ['x','y','z']: vel.append(float(velnode.find(tag).text)) vec = StateVector() vec.setTime(timestamp) vec.setPosition(pos) vec.setVelocity(vel) frameOrbit.addStateVector(vec) return frameOrbit def extractPreciseOrbit(self): ''' Extract precise orbit from given Orbit file. ''' try: fp = open(self.orbitFile,'r') except IOError as strerr: print("IOError: %s" % strerr) return _xml_root = ElementTree.ElementTree(file=fp).getroot() node = _xml_root.find('Data_Block/List_of_OSVs') orb = Orbit() orb.configure() margin = datetime.timedelta(seconds=40.0) tstart = self.product.sensingStart - margin tend = self.product.sensingStop + margin for child in node: timestamp = self.convertToDateTime(child.find('UTC').text[4:]) if (timestamp >= tstart) and (timestamp < tend): pos = [] vel = [] for tag in ['VX','VY','VZ']: vel.append(float(child.find(tag).text)) for tag in ['X','Y','Z']: pos.append(float(child.find(tag).text)) vec = StateVector() vec.setTime(timestamp) vec.setPosition(pos) vec.setVelocity(vel) # print(vec) orb.addStateVector(vec) fp.close() return orb def extractBetaLUT(self): ''' Extract Beta nought look up table from calibration file. ''' from scipy.interpolate import RectBivariateSpline if self.calibrationXml is None: raise Exception('No calibration file provided') if '.zip' in self.calibrationXml: import zipfile parts = self.calibrationXml.split(os.path.sep) zipname = os.path.join('/',*(parts[:-4])) fname = os.path.join(*(parts[-4:])) try: with zipfile.ZipFile(zipname, 'r') as zf: xmlstr = zf.read(fname) except: raise Exception('Could not read calibration file: {0}'.format(self.calibrationXml)) else: try: with open(self.calibrationXml, 'r') as fid: xmlstr = fid.read() except: raise Exception('Could not read calibration file: {0}'.format(self.calibrationXml)) _xml_root = ElementTree.fromstring(xmlstr) node = _xml_root.find('calibrationVectorList') num = int(node.attrib['count']) lines = [] pixels = [] data = None for ii, child in enumerate(node): pixnode = child.find('pixel') nump = int(pixnode.attrib['count']) if ii==0: data = np.zeros((num,nump)) pixels = [float(x) for x in pixnode.text.split()] lines.append( int(child.find('line').text)) signode = child.find('betaNought') data[ii,:] = [float(x) for x in signode.text.split()] lines = np.array(lines) pixels = np.array(pixels) self.betaLUT = RectBivariateSpline(lines,pixels, data, kx=1, ky=1) if False: import matplotlib.pyplot as plt plt.figure() plt.imshow(data) plt.colorbar() plt.show() return def extractNoiseLUT(self): ''' Extract Noise look up table from calibration file. ''' from scipy.interpolate import interp1d, InterpolatedUnivariateSpline if self.noiseXml is None: raise Exception('No calibration file provided') if '.zip' in self.noiseXml: import zipfile parts = self.noiseXml.split('.zip/') zipname = parts[0] + '.zip' fname = parts[1] try: with zipfile.ZipFile(zipname, 'r') as zf: xmlstr = zf.read(fname) except: raise Exception('Could not read noise file: {0}'.format(self.calibrationXml)) else: try: with open(self.noiseXml, 'r') as fid: xmlstr = fid.read() except: raise Exception('Could not read noise file: {0}'.format(self.calibrationXml)) if float(self.IPFversion) < 2.90: noise_range_vector_name = "noiseVectorList" noise_range_lut_name = "noiseLut" has_azimuth_noise_vectors = False self.noiseAzimuthLUT = None else: noise_range_vector_name = "noiseRangeVectorList" noise_range_lut_name = "noiseRangeLut" has_azimuth_noise_vectors = True print("Extracting noise LUT's...") _xml_root = ElementTree.fromstring(xmlstr) node = _xml_root.find(noise_range_vector_name) num_vectors = int(node.attrib['count']) print("File contains {} range noise vectors.".format(num_vectors)) full_samples_range = np.arange(self.product.numberOfSamples) noise_range_lut_indices = np.zeros((num_vectors,)) noise_range_lut_values = np.zeros((num_vectors, self.product.numberOfSamples)) for ii, child in enumerate(node): print("Processing range noise vector {}/{}".format(ii + 1, num_vectors)) pixnode = child.find('pixel') sample_pixels = [float(x) for x in pixnode.text.split()] signode = child.find(noise_range_lut_name) vector = np.asarray([float(x) for x in signode.text.split()]) vector_interpolator = InterpolatedUnivariateSpline(sample_pixels, vector, k=1) vector_interpolated = vector_interpolator(full_samples_range) noise_range_lut_indices[ii] = int(child.find('line').text) noise_range_lut_values[ii] = vector_interpolated self.noiseRangeLUT = interp1d(noise_range_lut_indices, noise_range_lut_values, kind='linear', axis=0, fill_value="extrapolate") if has_azimuth_noise_vectors: node = _xml_root.find("noiseAzimuthVectorList") num_vectors = int(node.attrib['count']) print("File contains {} azimuth noise blocks.".format(num_vectors)) noise_azimuth_lut_indices = defaultdict(list) noise_azimuth_lut_values = defaultdict(list) for block_i, child in enumerate(node): print("Processing azimuth noise vector {}/{}".format(block_i + 1, num_vectors)) linenode = child.find('line') signode = child.find("noiseAzimuthLut") block_range_start = int(child.find('firstRangeSample').text) block_range_end = int(child.find('lastRangeSample').text) block_azimuth_start = int(child.find('firstAzimuthLine').text) block_azimuth_end = int(child.find('lastAzimuthLine').text) block_line_index = [float(x) for x in linenode.text.split()] block_vector = [float(x) for x in signode.text.split()] block_line_range = np.arange(block_azimuth_start, block_azimuth_end + 1) if len(block_vector) > 1: block_vector_interpolator = InterpolatedUnivariateSpline(block_line_index, block_vector, k=1) for line in block_line_range: noise_azimuth_lut_indices[line].extend([block_range_start, block_range_end]) noise_azimuth_lut_values[line].extend([block_vector_interpolator(line)] * 2) else: for line in block_line_range: noise_azimuth_lut_indices[line].extend([block_range_start, block_range_end]) noise_azimuth_lut_values[line].extend([block_vector[0]] * 2) self.noiseAzimuthLUT = (noise_azimuth_lut_indices, noise_azimuth_lut_values) else: print("File contains no azimuth noise blocks.") if False: import matplotlib.pyplot as plt plt.figure() plt.imshow(data) plt.colorbar() plt.show() return def extractImage(self, parse=False, removeNoise=False): """ Use gdal python bindings to extract image """ try: from osgeo import gdal except ImportError: raise Exception('GDAL python bindings not found. Need this for RSAT2/ TandemX / Sentinel1.') from scipy.interpolate import interp1d if parse: self.parse() print('Extracting normalized image ....') src = gdal.Open('/vsizip//'+self.tiff.strip(), gdal.GA_ReadOnly) band = src.GetRasterBand(1) if self.product.numberOfSamples != src.RasterXSize: raise Exception('Width in metadata and image dont match') if self.product.numberOfLines != src.RasterYSize: raise Exception('Length in metadata and image dont match') noiseFactor = 0.0 if (not removeNoise) and self.noiseCorrectionApplied: print('User asked for data without noise corrections, but product appears to be corrected. Adding back noise from LUT ....') noiseFactor = 1.0 elif removeNoise and (not self.noiseCorrectionApplied): print('User asked for data with noise corrections, but the products appears to not be corrected. Applying noise corrections from LUT ....') noiseFactor = -1.0 elif (not removeNoise) and (not self.noiseCorrectionApplied): print('User asked for data without noise corrections. The product contains uncorrected data ... unpacking ....') else: print('User asked for noise corrected data and the product appears to be noise corrected .... unpacking ....') ###Write original SLC to file fid = open(self.output, 'wb') pix = np.arange(self.product.numberOfSamples) if self.noiseAzimuthLUT is not None: noise_azimuth_lut_indices, noise_azimuth_lut_values = self.noiseAzimuthLUT for ii in range(self.product.numberOfLines//100 + 1): ymin = int(ii*100) ymax = int(np.clip(ymin+100,0, self.product.numberOfLines)) if ymax == ymin: break lin = np.arange(ymin,ymax) ####Read in one block of data data = 1.0 * band.ReadAsArray(0, ymin, self.product.numberOfSamples, ymax-ymin) lut = self.betaLUT(lin,pix,grid=True) if noiseFactor != 0.0: noise = self.noiseRangeLUT(lin) if self.noiseAzimuthLUT is not None: block_azimuth_noise = np.zeros_like(noise) for l_i, line in enumerate(lin): interpolator = interp1d(noise_azimuth_lut_indices[line], noise_azimuth_lut_values[line], kind='previous', fill_value="extrapolate") block_azimuth_noise[l_i] = interpolator(np.arange(self.product.numberOfSamples)) noise *= block_azimuth_noise else: noise = 0.0 #outdata = data outdata = np.clip(data*data/(lut*lut) + noiseFactor * noise/(lut*lut), 0, None) #outdata = 10 * np.log10(outdata) outdata.astype(np.float32).tofile(fid) fid.close() ####Render ISCE XML slcImage = isceobj.createImage() slcImage.setByteOrder('l') slcImage.dataType = 'FLOAT' slcImage.setFilename(self.output) slcImage.setAccessMode('read') slcImage.setWidth(self.product.numberOfSamples) slcImage.setLength(self.product.numberOfLines) slcImage.renderHdr() slcImage.renderVRT() self.product.image = slcImage band = None src = None return def extractSlantRange(self, full=True): ''' Extract pixel-by-pixel slant range file for GRD images. ''' print('Extracing slant range ....') from scipy.interpolate import interp1d from isceobj.Util import Poly1D node = self._xml_root.find('coordinateConversion/coordinateConversionList') num = int(node.attrib['count']) lines = [] data = [] for ii, child in enumerate(node): t0 = self.convertToDateTime(child.find('azimuthTime').text) lines.append( (t0-self.product.sensingStart).total_seconds()/self.product.azimuthTimeInterval) pp = [float(x) for x in child.find('grsrCoefficients').text.split()] meangr = float( child.find('gr0').text) if meangr !=0 : raise Exception('Ground range image does not start at zero coordinate') data.append(pp[::-1]) ###If polynomial starts beyond first line if lines[0] > 0: lines.insert(0, 0.) pp = data[0] data.insert(0,pp) ####If polynomial ends before last line if lines[-1] < (self.product.numberOfLines-1): lines.append(self.product.numberOfLines-1.0) pp = data[-1] data.append(pp) lines = np.array(lines) data = np.array(data) LUT = [] for ii in range(data.shape[1]): col = data[:,ii] LUT.append(interp1d(lines, col, bounds_error=False, assume_sorted=True)) self.gr2srLUT = LUT ###Write original SLC to file fid = open(self.slantRangeFile, 'wb') pix = np.arange(self.product.numberOfSamples) * self.product.groundRangePixelSize lin = np.arange(self.product.numberOfLines) polys = np.zeros((self.product.numberOfLines, len(self.gr2srLUT))) for ii, intp in enumerate(self.gr2srLUT): polys[:,ii] = intp(lin) minrng = 1e11 maxrng = -1e11 for ii in range(self.product.numberOfLines): pp = polys[ii,:] outdata = np.polyval(pp, pix) minrng = min(minrng, outdata[0]) ###Can be made more efficient maxrng = max(maxrng, outdata[-1]) outdata.tofile(fid) fid.close() self.product.startingSlantRange = minrng self.product.endingSlantRange = maxrng ####Render ISCE XML slcImage = isceobj.createImage() slcImage.setByteOrder('l') slcImage.dataType = 'DOUBLE' slcImage.setFilename(self.slantRangeFile) slcImage.setAccessMode('read') slcImage.setWidth(self.product.numberOfSamples) slcImage.setLength(self.product.numberOfLines) slcImage.renderHdr() slcImage.renderVRT() self.product.slantRangeImage = slcImage