271 lines
8.8 KiB
Python
Executable File
271 lines
8.8 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
########################
|
|
#Author: Heresh Fattahi
|
|
#Copyright 2016
|
|
######################
|
|
import argparse
|
|
import isce
|
|
import isceobj
|
|
import os
|
|
import gdal
|
|
import numpy as np
|
|
import xml.etree.ElementTree as ET
|
|
|
|
def createParser():
|
|
'''
|
|
Create command line parser.
|
|
'''
|
|
|
|
parser = argparse.ArgumentParser( description='Create DEM simulation for merged images')
|
|
parser.add_argument('-l','--lat', dest='latFile', type=str, required=True,
|
|
help = 'latitude file in radar coordinate')
|
|
parser.add_argument('-L','--lon', dest='lonFile', type=str, required=True,
|
|
help = 'longitude file in radar coordinate')
|
|
parser.add_argument('-f', '--filelist', dest='prodlist', type=str, required=True,
|
|
help='Input file to be geocoded')
|
|
parser.add_argument('-o', '--xoff', dest='xOff', type=int, default=0,
|
|
help='Offset from the begining of geometry files in x direction. Default 0.0')
|
|
parser.add_argument('-p', '--yoff', dest='yOff', type=int, default=0,
|
|
help='Offset from the begining of geometry files in y direction. Default 0.0')
|
|
parser.add_argument('-r', '--resampling_method', dest='resamplingMethod', type=str, default='near',
|
|
help='Resampling method (gdalwarp resamplin methods)')
|
|
|
|
|
|
parser.add_argument('-b', '--bbox', dest='bbox', type=str, default='',
|
|
help='Bounding box (SNWE)')
|
|
parser.add_argument('-x', '--lon_step', dest='lonStep', type=str, default=0.001,
|
|
help='output pixel size (longitude) in degrees. Default 0.001')
|
|
parser.add_argument('-y', '--lat_step', dest='latStep', type=str, default=0.001,
|
|
help='output pixel size (latitude) in degrees. Default 0.001')
|
|
|
|
|
|
parser.add_argument('-t', '--tiff', dest='istiff', action='store_true', default=False,
|
|
help='Create GEOTIFF instead of standard ENVI / ISCE files')
|
|
parser.add_argument('--alex', dest='isAlexGrid', default=False,
|
|
action='store_true', help='Geocode to the Antarctica grid for optical offsets used by Alex Gardner')
|
|
|
|
|
|
return parser
|
|
|
|
def cmdLineParse(iargs = None):
|
|
'''
|
|
Command line parser.
|
|
'''
|
|
|
|
parser = createParser()
|
|
inps = parser.parse_args(args = iargs)
|
|
|
|
|
|
if not inps.isAlexGrid:
|
|
inps.bbox = [val for val in inps.bbox.split()]
|
|
if len(inps.bbox) != 4:
|
|
raise Exception('Bbox should contain 4 floating point values')
|
|
inps.outproj = 'EPSG:4326'
|
|
else:
|
|
print('Ignoring bbox and spacing inputs. Using standard grid for Antarctica.')
|
|
inps.lonStep = 240.0
|
|
inps.latStep = 240.0
|
|
inps.outproj = 'EPSG:3031'
|
|
|
|
inps.prodlist = inps.prodlist.split()
|
|
return inps
|
|
|
|
def prepare_lat_lon(inps):
|
|
|
|
latFile = os.path.abspath(inps.latFile)
|
|
lonFile = os.path.abspath(inps.lonFile)
|
|
#cmd = 'isce2gis.py vrt -i ' + latFile
|
|
#os.system(cmd)
|
|
#cmd = 'isce2gis.py vrt -i ' + lonFile
|
|
#os.system(cmd)
|
|
|
|
|
|
width, length = getSize(latFile)
|
|
widthFile , lengthFile = getSize(inps.prodlist[0])
|
|
|
|
print("size of lat and lon files (width, length) ", width, length)
|
|
print("size of input file to be geocoded (width, length): ", widthFile , lengthFile)
|
|
|
|
xOff = inps.xOff
|
|
yOff = inps.yOff
|
|
|
|
cmd = 'gdal_translate -of VRT -srcwin ' + str(xOff) + ' ' + str(yOff) \
|
|
+' '+ str(width - xOff) +' '+ str(length - yOff) +' -outsize ' + str(widthFile) + \
|
|
' '+ str(lengthFile) + ' -a_nodata 0 ' + latFile +'.vrt' + ' tempLAT.vrt'
|
|
|
|
os.system(cmd)
|
|
|
|
cmd = 'gdal_translate -of VRT -srcwin ' + str(xOff) + ' ' + str(yOff) \
|
|
+' '+ str(int(width-xOff)) +' '+ str(int(length-yOff)) +' -outsize ' + str(widthFile) +\
|
|
' '+ str(lengthFile) + ' -a_nodata 0 ' + lonFile +'.vrt' + ' tempLON.vrt'
|
|
|
|
os.system(cmd)
|
|
|
|
return 'tempLAT.vrt', 'tempLON.vrt'
|
|
|
|
# gdal_translate -of VRT -srcwin 384 384 64889 12785 -outsize 1013 199 ../../COMBINED/GEOM_MASTER/LAT.rdr LAT_off.vrt
|
|
|
|
|
|
def writeVRT(infile, latFile, lonFile):
|
|
#This function is modified from isce2gis.py
|
|
|
|
#cmd = 'isce2gis.py vrt -i ' + infile
|
|
#os.system(cmd)
|
|
|
|
tree = ET.parse(infile + '.vrt')
|
|
root = tree.getroot()
|
|
|
|
meta = ET.SubElement(root, 'metadata')
|
|
meta.attrib['domain'] = "GEOLOCATION"
|
|
meta.tail = '\n'
|
|
meta.text = '\n '
|
|
|
|
|
|
rdict = { 'Y_DATASET' : os.path.relpath(latFile , os.path.dirname(infile)),
|
|
'X_DATASET' : os.path.relpath(lonFile , os.path.dirname(infile)),
|
|
'X_BAND' : "1",
|
|
'Y_BAND' : "1",
|
|
'PIXEL_OFFSET': "0",
|
|
'LINE_OFFSET' : "0",
|
|
'LINE_STEP' : "1",
|
|
'PIXEL_STEP' : "1" }
|
|
|
|
for key, val in rdict.items():
|
|
data = ET.SubElement(meta, 'mdi')
|
|
data.text = val
|
|
data.attrib['key'] = key
|
|
data.tail = '\n '
|
|
|
|
data.tail = '\n'
|
|
tree.write(infile + '.vrt')
|
|
|
|
|
|
def runGeo(inps):
|
|
|
|
#for file in inps.prodlist:
|
|
#cmd = 'isce2gis.py envi -i ' + file
|
|
#os.system(cmd)
|
|
|
|
|
|
if not inps.isAlexGrid:
|
|
WSEN = str(inps.bbox[2]) + ' ' + str(inps.bbox[0]) + ' ' + str(inps.bbox[3]) + ' ' + str(inps.bbox[1])
|
|
latFile, lonFile = prepare_lat_lon(inps)
|
|
|
|
getBound(latFile,float(inps.bbox[0]),float(inps.bbox[1]),'lat')
|
|
getBound(lonFile,float(inps.bbox[2]),float(inps.bbox[3]),'lon')
|
|
|
|
for infile in inps.prodlist:
|
|
infile = os.path.abspath(infile)
|
|
print ('geocoding ' + infile)
|
|
outFile = os.path.join(os.path.dirname(infile), "geo_" + os.path.basename(infile))
|
|
#cmd = 'isce2gis.py vrt -i '+ file + ' --lon ' + lonFile + ' --lat '+ latFile
|
|
#os.system(cmd)
|
|
writeVRT(infile, latFile, lonFile)
|
|
|
|
cmd = 'gdalwarp -of ENVI -geoloc -te '+ WSEN + ' -tr ' + str(inps.latStep) + ' ' + str(inps.lonStep) + ' -srcnodata 0 -dstnodata 0 ' + ' -r ' + inps.resamplingMethod + ' ' + infile +'.vrt '+ outFile
|
|
print (cmd)
|
|
os.system(cmd)
|
|
|
|
write_xml(outFile)
|
|
|
|
|
|
else:
|
|
from geo2ant import getGridLimits
|
|
ylims, xlims = getGridLimits(latfile=latFile, lonfile=lonFile)
|
|
|
|
WSEN = str(xlim[0]) + ' ' + str(ylim[0]) + ' ' + str(xlim[1]) + ' ' + str(ylim[1])
|
|
if inps.istiff:
|
|
ext = '.tif'
|
|
outformat = 'GTiff'
|
|
else:
|
|
ext = '.ant'
|
|
outformat = 'ENVI'
|
|
|
|
for infile in inps.prodlist:
|
|
print('geocoding: ' + infile)
|
|
|
|
writeVRT(infile, latFile, lonFile)
|
|
|
|
cmd = 'gdalwarp -of ' + outformat + ' -t_srs ' + inps.outproj + ' -geoloc -te ' + WSEN + ' -tr ' + str(inps.lonStep) + ' ' + str(inps.latStep) + ' -srcnodata 0 -dstnodata 0 -r ' + inps.resamplingMethod + ' ' + infile + '.vrt ' + infile+'ext'
|
|
status = os.system(cmd)
|
|
if status:
|
|
raise Exception('Command {0} Failed'.format(cmd))
|
|
|
|
if not inps.istiff:
|
|
write_xml(infile+ext)
|
|
|
|
|
|
def getSize(infile):
|
|
|
|
ds=gdal.Open(infile + ".vrt")
|
|
b=ds.GetRasterBand(1)
|
|
return b.XSize, b.YSize
|
|
|
|
def getBound(infile,minval,maxval,latlon): #added by Minyan Zhong
|
|
|
|
ds=gdal.Open(infile)
|
|
b=ds.GetRasterBand(1)
|
|
data=b.ReadAsArray()
|
|
|
|
idx=np.where((data>=minval) & (data<=maxval))
|
|
|
|
if latlon=='lat':
|
|
print('latitide bound in cliped area:')
|
|
else:
|
|
print('longitude bound in cliped area:')
|
|
print(np.min(data[idx]),np.max(data[idx]))
|
|
|
|
|
|
def get_lat_lon(infile):
|
|
|
|
ds=gdal.Open(infile)
|
|
b=ds.GetRasterBand(1)
|
|
width = b.XSize
|
|
length = b.YSize
|
|
minLon = ds.GetGeoTransform()[0]
|
|
deltaLon = ds.GetGeoTransform()[1]
|
|
maxLat = ds.GetGeoTransform()[3]
|
|
deltaLat = ds.GetGeoTransform()[5]
|
|
minLat = maxLat + (b.YSize)*deltaLat
|
|
|
|
return maxLat, deltaLat, minLon, deltaLon, width, length
|
|
|
|
def write_xml(outFile):
|
|
|
|
maxLat, deltaLat, minLon, deltaLon, width, length = get_lat_lon(outFile)
|
|
|
|
unwImage = isceobj.Image.createImage()
|
|
unwImage.setFilename(outFile)
|
|
unwImage.setWidth(width)
|
|
unwImage.setLength(length)
|
|
unwImage.bands = 1
|
|
unwImage.scheme = 'BIL'
|
|
unwImage.dataType = 'FLOAT'
|
|
unwImage.setAccessMode('read')
|
|
|
|
unwImage.coord2.coordDescription = 'Latitude'
|
|
unwImage.coord2.coordUnits = 'degree'
|
|
unwImage.coord2.coordStart = maxLat
|
|
unwImage.coord2.coordDelta = deltaLat
|
|
unwImage.coord1.coordDescription = 'Longitude'
|
|
unwImage.coord1.coordUnits = 'degree'
|
|
unwImage.coord1.coordStart = minLon
|
|
unwImage.coord1.coordDelta = deltaLon
|
|
|
|
# unwImage.createImage()
|
|
unwImage.renderHdr()
|
|
unwImage.renderVRT()
|
|
|
|
def main(iargs=None):
|
|
'''
|
|
Main driver.
|
|
'''
|
|
inps = cmdLineParse(iargs)
|
|
runGeo(inps)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|
|
|
|
|