ISCE_INSAR/contrib/geo_autoRIFT/geogrid/GeogridOptical.py

#!/usr/bin/env python3

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Copyright 2019 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.
#
# Authors: Piyush Agram, Yang Lei
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~






import pdb
import subprocess
import re
import string

class GeogridOptical():
    '''
    Class for mapping regular geographic grid on radar imagery.
    '''

    def runGeogrid(self):
        '''
        Do the actual processing.
        '''
        
        from . import geogridOptical

        ##Determine appropriate EPSG system
        self.epsgDem = self.getProjectionSystem(self.demname)
        self.epsgDat = self.getProjectionSystem(self.dat1name)

        ###Determine extent of data needed
        bbox = self.determineBbox()
        
        ##Create and set parameters
        self.setState()
        
        ##check parameters
        self.checkState()

        ##Run
        geogridOptical.geogridOptical_Py(self._geogridOptical)
        self.get_center_latlon()
    
        ##Get parameters
        self.getState()
    
        ##Clean up
        self.finalize()
    
    def get_center_latlon(self):
        '''
        Get center lat/lon of the image.
        '''
        from osgeo import gdal
        gdal.AllRegister()
        self.epsgDem = 4326
        self.epsgDat = self.getProjectionSystem(self.dat1name)
        self.determineBbox()
        if gdal.__version__[0] == '2':
            self.cen_lat = (self._ylim[0] + self._ylim[1]) / 2
            self.cen_lon = (self._xlim[0] + self._xlim[1]) / 2
        else:
            self.cen_lon = (self._ylim[0] + self._ylim[1]) / 2
            self.cen_lat = (self._xlim[0] + self._xlim[1]) / 2
        print("Scene-center lat/lon: " + str(self.cen_lat) + "  " + str(self.cen_lon))
    

    def getProjectionSystem(self, filename):
        '''
        Testing with Greenland.
        '''
        if not filename:
            raise Exception('File {0} does not exist'.format(filename))

        from osgeo import gdal, osr
        ds = gdal.Open(filename, gdal.GA_ReadOnly)
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection())
        srs.AutoIdentifyEPSG()
        ds = None
#        pdb.set_trace()

        if srs.IsGeographic():
            epsgstr = srs.GetAuthorityCode('GEOGCS')
        elif srs.IsProjected():
            epsgstr = srs.GetAuthorityCode('PROJCS')
        elif srs.IsLocal():
            raise Exception('Local coordinate system encountered')
        else:
            raise Exception('Non-standard coordinate system encountered')
        if not epsgstr:  #Empty string->use shell command gdalsrsinfo for last trial
            cmd = 'gdalsrsinfo -o epsg {0}'.format(filename)
            epsgstr = subprocess.check_output(cmd, shell=True)
#            pdb.set_trace()
            epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0]
#            pdb.set_trace()
        if not epsgstr:  #Empty string
            raise Exception('Could not auto-identify epsg code')
#        pdb.set_trace()
        epsgcode = int(epsgstr)
#        pdb.set_trace()
        return epsgcode

    def determineBbox(self, zrange=[-200,4000]):
        '''
        Dummy.
        '''
        import numpy as np
        import datetime
        from osgeo import osr

#        import pdb
#        pdb.set_trace()


        samples = self.startingX + np.array([0, self.numberOfSamples-1]) * self.XSize
        lines = self.startingY + np.array([0, self.numberOfLines-1]) * self.YSize

        coordDat = osr.SpatialReference()
        if self.epsgDat:
            coordDat.ImportFromEPSG(self.epsgDat)
        else:
            raise Exception('EPSG code does not exist for image data')


        coordDem = osr.SpatialReference()
        if self.epsgDem:
            coordDem.ImportFromEPSG(self.epsgDem)
        else:
            raise Exception('EPSG code does not exist for DEM')


        trans = osr.CoordinateTransformation(coordDat, coordDem)



        utms = []
        xyzs = []


        ### Four corner coordinates
        for ss in samples:
            for ll in lines:
                for zz in zrange:
                    utms.append([ss,ll,zz])
                    x,y,z = trans.TransformPoint(ss, ll, zz)
                    xyzs.append([x,y,z])

        utms = np.array(utms)
        xyzs = np.array(xyzs)

        self._xlim = [np.min(xyzs[:,0]), np.max(xyzs[:,0])]
        self._ylim = [np.min(xyzs[:,1]), np.max(xyzs[:,1])]


    def getState(self):

        from . import geogridOptical
        
        self.pOff = geogridOptical.getXOff_Py(self._geogridOptical)
        self.lOff = geogridOptical.getYOff_Py(self._geogridOptical)
        self.pCount = geogridOptical.getXCount_Py(self._geogridOptical)
        self.lCount = geogridOptical.getYCount_Py(self._geogridOptical)
        self.X_res = geogridOptical.getXPixelSize_Py(self._geogridOptical)
        self.Y_res = geogridOptical.getYPixelSize_Py(self._geogridOptical)
    
    def setState(self):
        '''
        Create C object and populate.
        '''

        from . import geogridOptical

        if self._geogridOptical is not None:
            geogridOptical.destroyGeoGridOptical_Py(self._geogridOptical)
    
        self._geogridOptical = geogridOptical.createGeoGridOptical_Py()
        geogridOptical.setOpticalImageDimensions_Py( self._geogridOptical, self.numberOfSamples, self.numberOfLines)
        geogridOptical.setXParameters_Py( self._geogridOptical, self.startingX, self.XSize)
        geogridOptical.setYParameters_Py( self._geogridOptical, self.startingY, self.YSize)
        geogridOptical.setRepeatTime_Py(self._geogridOptical, self.repeatTime)
        
        geogridOptical.setDtUnity_Py( self._geogridOptical, self.srs_dt_unity)
        geogridOptical.setMaxFactor_Py( self._geogridOptical, self.srs_max_scale)
        geogridOptical.setUpperThreshold_Py( self._geogridOptical, self.srs_max_search)
        geogridOptical.setLowerThreshold_Py(self._geogridOptical, self.srs_min_search)
        
        geogridOptical.setEPSG_Py(self._geogridOptical, self.epsgDem, self.epsgDat)
        geogridOptical.setChipSizeX0_Py(self._geogridOptical, self.chipSizeX0)
        geogridOptical.setGridSpacingX_Py(self._geogridOptical, self.gridSpacingX)
        
        geogridOptical.setXLimits_Py(self._geogridOptical, self._xlim[0], self._xlim[1])
        geogridOptical.setYLimits_Py(self._geogridOptical, self._ylim[0], self._ylim[1])
        if self.demname:
            geogridOptical.setDEM_Py(self._geogridOptical, self.demname)

        if (self.dhdxname is not None) and (self.dhdyname is not None):
            geogridOptical.setSlopes_Py(self._geogridOptical, self.dhdxname, self.dhdyname)
        
        if (self.vxname is not None) and (self.vyname is not None):
            geogridOptical.setVelocities_Py(self._geogridOptical, self.vxname, self.vyname)
        
        if (self.srxname is not None) and (self.sryname is not None):
            geogridOptical.setSearchRange_Py(self._geogridOptical, self.srxname, self.sryname)

        if (self.csminxname is not None) and (self.csminyname is not None):
            geogridOptical.setChipSizeMin_Py(self._geogridOptical, self.csminxname, self.csminyname)
        
        if (self.csmaxxname is not None) and (self.csmaxyname is not None):
            geogridOptical.setChipSizeMax_Py(self._geogridOptical, self.csmaxxname, self.csmaxyname)

        if (self.ssmname is not None):
            geogridOptical.setStableSurfaceMask_Py(self._geogridOptical, self.ssmname)
        
        geogridOptical.setWindowLocationsFilename_Py( self._geogridOptical, self.winlocname)
        geogridOptical.setWindowOffsetsFilename_Py( self._geogridOptical, self.winoffname)
        geogridOptical.setWindowSearchRangeFilename_Py( self._geogridOptical, self.winsrname)
        geogridOptical.setWindowChipSizeMinFilename_Py( self._geogridOptical, self.wincsminname)
        geogridOptical.setWindowChipSizeMaxFilename_Py( self._geogridOptical, self.wincsmaxname)
        geogridOptical.setWindowStableSurfaceMaskFilename_Py( self._geogridOptical, self.winssmname)
        geogridOptical.setRO2VXFilename_Py( self._geogridOptical, self.winro2vxname)
        geogridOptical.setRO2VYFilename_Py( self._geogridOptical, self.winro2vyname)
        geogridOptical.setNodataOut_Py(self._geogridOptical, self.nodata_out)
        
    
    def checkState(self):
        '''
        Create C object and populate.
        '''
        if self.repeatTime < 0:
            raise Exception('Input image 1 must be older than input image 2')


    def finalize(self):
        '''
        Clean up all the C pointers.
        '''

        from . import geogridOptical
        
        geogridOptical.destroyGeoGridOptical_Py(self._geogridOptical)
        self._geogridOptical = None


    






    def coregister(self,in1,in2):
        import os
        import numpy as np

        from osgeo import gdal, osr
        import struct

        DS1 = gdal.Open(in1, gdal.GA_ReadOnly)
        trans1 = DS1.GetGeoTransform()
        xsize1 = DS1.RasterXSize
        ysize1 = DS1.RasterYSize
        epsg1 = self.getProjectionSystem(in1)

        DS2 = gdal.Open(in2, gdal.GA_ReadOnly)
        trans2 = DS2.GetGeoTransform()
        xsize2 = DS2.RasterXSize
        ysize2 = DS2.RasterYSize
        epsg2 = self.getProjectionSystem(in2)
        
        if epsg1 != epsg2:
            raise Exception('The current version of geo_autoRIFT assumes the two images are in the same projection, i.e. it cannot handle two different projections; the users are thus recommended to do the tranformation themselves before running geo_autoRIFT.')
        
        

        W = np.max([trans1[0],trans2[0]])
        N = np.min([trans1[3],trans2[3]])
        E = np.min([trans1[0]+(xsize1-1)*trans1[1],trans2[0]+(xsize2-1)*trans2[1]])
        S = np.max([trans1[3]+(ysize1-1)*trans1[5],trans2[3]+(ysize2-1)*trans2[5]])

        x1a = int(np.round((W-trans1[0])/trans1[1]))
        x1b = int(np.round((E-trans1[0])/trans1[1]))
        y1a = int(np.round((N-trans1[3])/trans1[5]))
        y1b = int(np.round((S-trans1[3])/trans1[5]))

        x2a = int(np.round((W-trans2[0])/trans2[1]))
        x2b = int(np.round((E-trans2[0])/trans2[1]))
        y2a = int(np.round((N-trans2[3])/trans2[5]))
        y2b = int(np.round((S-trans2[3])/trans2[5]))
        
        if (x1a > (xsize1-1))|(x1b > (xsize1-1))|(x2a > (xsize2-1))|(x2b > (xsize2-1))|(y1a > (ysize1-1))|(y1b > (ysize1-1))|(y2a > (ysize2-1))|(y2b > (ysize2-1)):
            raise Exception('Uppper bound of coregistered image index should be <= size of image1 (and image2) minus 1')

        if (x1a < 0)|(x1b < 0)|(x2a < 0)|(x2b < 0)|(y1a < 0)|(y1b < 0)|(y2a < 0)|(y2b < 0):
            raise Exception('Lower bound of coregistered image index should be >= 0')

        if ((x1b-x1a) != (x2b-x2a))|((y1b-y1a) != (y2b-y2a)):
            raise Exception('Coregistered image size mismatch between image1 and image2')

        x1a = int(x1a)
        x1b = int(x1b)
        y1a = int(y1a)
        y1b = int(y1b)
        x2a = int(x2a)
        x2b = int(x2b)
        y2a = int(y2a)
        y2b = int(y2b)

        trans = (W, trans1[1], 0.0, N, 0.0, trans1[5])

        return x1a, y1a, x1b-x1a+1, y1b-y1a+1, x2a, y2a, x2b-x2a+1, y2b-y2a+1, trans








    def __init__(self):
        super(GeogridOptical, self).__init__()

        ##Optical image related parameters
        self.startingY = None
        self.startingX = None
        self.XSize = None
        self.YSize = None
        self.numberOfSamples = None
        self.numberOfLines = None
        self.repeatTime = None
        self.chipSizeX0 = None
        self.gridSpacingX = None

        ##Input related parameters
        self.dat1name = None
        self.demname = None
        self.dhdxname = None
        self.dhdyname = None
        self.vxname = None
        self.vyname = None
        self.srxname = None
        self.sryname = None
        self.csminxname = None
        self.csminyname = None
        self.csmaxxname = None
        self.csmaxyname = None
        self.ssmname = None

        ##Output related parameters
        self.winlocname = None
        self.winoffname = None
        self.winsrname = None
        self.wincsminname = None
        self.wincsmaxname = None
        self.winssmname = None
        self.winro2vxname = None
        self.winro2vyname = None
        
        ##dt-varying search range scale (srs) rountine parameters
        self.srs_dt_unity = 182
        self.srs_max_scale = 5
        self.srs_max_search = 20000
        self.srs_min_search = 0

        ##Coordinate system
        self.epsgDem = None
        self.epsgDat = None
        self._xlim = None
        self._ylim = None
        self.nodata_out = None
        
        ##Pointer to C
        self._geogridOptical = None

        ##parameters for autoRIFT
        self.pOff = None
        self.lOff = None
        self.pCount = None
        self.lCount = None
        self.X_res = None
        self.Y_res = None