ISCE_INSAR/components/isceobj/Sensor/GRD/Sentinel1.py

#!/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.getchildren():
            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.getchildren():
            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.getchildren():
            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.getchildren()):
            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.getchildren()):
            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.getchildren()):
                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.getchildren()):
            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