106 lines
4.5 KiB
Python
106 lines
4.5 KiB
Python
#!/usr/bin/env python3
|
|
#
|
|
|
|
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
# Copyright 2017 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: David Bekaert
|
|
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
|
|
|
|
import os
|
|
import sys
|
|
import argparse
|
|
from osgeo import gdal
|
|
from isce.applications.gdal2isce_xml import gdal2isce_xml
|
|
|
|
|
|
# command line parsing of input file
|
|
def cmdLineParse():
|
|
'''
|
|
Command line parser.
|
|
'''
|
|
parser = argparse.ArgumentParser(description='Generate down-sample DEM from wgs84.vrt DEM file')
|
|
parser.add_argument('-i','--input', dest='input_dem_vrt', type=str, required=True, help='Input DEM vrt filename (GDAL supported)')
|
|
parser.add_argument('-rmeter','--res_meter', dest='res_meter', type=str, default='', required=False, help='DEM output resolution in m units')
|
|
parser.add_argument('-rsec','--res_seconds', dest='res_seconds', type=str, default ='', required=False, help='DEM output resolution in arc seconds units')
|
|
return parser.parse_args()
|
|
|
|
|
|
# main script
|
|
if __name__ == '__main__':
|
|
'''
|
|
Main driver.
|
|
'''
|
|
|
|
# Parse command line
|
|
inps = cmdLineParse()
|
|
|
|
if inps.res_meter == '' and inps.res_seconds == '':
|
|
raise Exception('Provide either rmeter or rsec argument for DEM resolution')
|
|
|
|
# check if the input file exist
|
|
if not os.path.isfile(inps.input_dem_vrt):
|
|
raise Exception('Input file is not found ....')
|
|
# check if the provided input file is a .vrt file and also get the envi filename
|
|
input_dem_envi, file_extension = os.path.splitext(inps.input_dem_vrt)
|
|
if file_extension != '.vrt':
|
|
raise Exception('Input file is not a vrt file ....')
|
|
# get the file path
|
|
input_path = os.path.dirname(os.path.abspath(inps.input_dem_vrt))
|
|
|
|
|
|
# convert the output resolution from m in degrees
|
|
# (this is approximate, could use instead exact expression)
|
|
if inps.res_meter != '':
|
|
gdal_opts = gdal.WarpOptions(format='ENVI',
|
|
outputType=gdal.GDT_Int16,
|
|
dstSRS='EPSG:4326',
|
|
xRes=float(inps.res_meter)/110/1000,
|
|
yRes=float(inps.res_meter)/110/1000,
|
|
targetAlignedPixels=True)
|
|
# res_degree = float(inps.res_meter)/110/1000
|
|
elif inps.res_seconds != '':
|
|
gdal_opts = gdal.WarpOptions(format='ENVI',
|
|
outputType=gdal.GDT_Int16,
|
|
dstSRS='EPSG:4326',
|
|
xRes=float(inps.res_seconds)*1/60*1/60,
|
|
yRes=float(inps.res_seconds)*1/60*1/60,
|
|
targetAlignedPixels=True)
|
|
# res_degree = float(1/60*1/60*float(inps.res_seconds))
|
|
|
|
# The ENVI filename of the coarse DEM to be generated
|
|
coarse_dem_envi = os.path.join(input_path, "Coarse_" + input_dem_envi)
|
|
|
|
# Using gdal to down-sample the WGS84 DEM
|
|
# cmd = "gdalwarp -t_srs EPSG:4326 -ot Int16 -of ENVI -tap -tr " + str(res_degree) + " " + str(res_degree) + " " + inps.input_dem_vrt + " " + coarse_dem_envi
|
|
# os.system(cmd)
|
|
ds = gdal.Warp(coarse_dem_envi,inps.input_dem_vrt,options=gdal_opts)
|
|
ds = None
|
|
|
|
# Generating the ISCE xml and vrt of this coarse DEM
|
|
gdal2isce_xml(coarse_dem_envi+'.vrt')
|
|
|