407 lines
14 KiB
Python
Executable File
407 lines
14 KiB
Python
Executable File
#!/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
|