1008 lines
42 KiB
Python
Executable File
1008 lines
42 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
|
|
|
|
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
# Copyright 2012 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: Giangi Sacco
|
|
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
|
|
|
|
|
|
|
|
|
|
from __future__ import print_function
|
|
import isce
|
|
from ctypes import cdll, c_char_p, c_int, byref
|
|
from array import array
|
|
import struct
|
|
import zipfile
|
|
import os
|
|
import sys
|
|
import math
|
|
import urllib.request, urllib.parse, urllib.error
|
|
from isce import logging
|
|
from iscesys.Component.Component import Component
|
|
|
|
import xml.etree.ElementTree as ET
|
|
from html.parser import HTMLParser
|
|
class DemDirParser(HTMLParser):
|
|
def __init__(self):
|
|
HTMLParser.__init__(self)
|
|
self._results = []
|
|
self._filterList = []
|
|
self._removeList = []
|
|
@property
|
|
def filterList(self):
|
|
return self._filterList
|
|
@filterList.setter
|
|
def filterList(self,filterList):
|
|
self._filterList = filterList
|
|
@property
|
|
def removeList(self):
|
|
return self._removeList
|
|
@removeList.setter
|
|
def removeList(self,removeList):
|
|
self._removeList = removeList
|
|
@property
|
|
def results(self):
|
|
return self._results
|
|
#implement the call back from data received
|
|
def handle_data(self,data):
|
|
#check that the data is one of the expected type
|
|
#based on filtesList
|
|
for filt in self.filterList:
|
|
isOk = True
|
|
#check that the data is not one that needs to be removed
|
|
for rm in self.removeList:
|
|
if data.count(rm):
|
|
isOk = False
|
|
break
|
|
if isOk and data.count(filt):
|
|
self._results.append(data.strip())
|
|
|
|
|
|
|
|
#Parameters definitions
|
|
URL = Component.Parameter('_url',
|
|
public_name = 'URL',default = 'http://dds.cr.usgs.gov',
|
|
type = str,
|
|
mandatory = False,
|
|
doc = "Top part of the url where the DEMs are stored. Used for SRTM version2")
|
|
USERNAME = Component.Parameter('_un',
|
|
public_name='username',
|
|
default = None,
|
|
type = str,
|
|
mandatory = False,
|
|
doc = "Username in case the url is password protected")
|
|
PASSWORD = Component.Parameter('_pw',
|
|
public_name='password',
|
|
default = None,
|
|
type = str,
|
|
mandatory = False,
|
|
doc = "Password in case the url is password protected")
|
|
KEEP_AFTER_FAILED = Component.Parameter('_keepAfterFailed',
|
|
public_name='keepAfterFailed',
|
|
default = False,
|
|
type = bool,
|
|
mandatory = False,
|
|
doc = "If the stitching for some reason fails, it keeps the downloaded files.\n" +\
|
|
"If 'useLocalDirectory' is set then this flag is forced to True to avoid \n" +\
|
|
"accidental deletion of files (default: False)")
|
|
DIRECTORY = Component.Parameter('_downloadDir',
|
|
public_name='directory',
|
|
default = './',
|
|
type = str,
|
|
mandatory = False,
|
|
doc = "If useLocalDirectory is False,it is used to download\n" + \
|
|
"the files and create the stitched file, otherwise it assumes that this is the\n" + \
|
|
"the local directory where the DEMs are (default: current working directory)")
|
|
ACTION = Component.Parameter('_action',
|
|
public_name='action',
|
|
default = 'stitch',
|
|
type = str,
|
|
mandatory = False,
|
|
doc = "Action to perform. Possible values are 'stitch' to stitch DEMs together\n" + \
|
|
"or 'download' to download the DEMs (default: 'stitch')")
|
|
CORRECT = Component.Parameter('_correct',
|
|
public_name='correct',
|
|
default = False,
|
|
type = bool,
|
|
mandatory = False,
|
|
doc = "Apply correction EGM96 -> WGS84 (default: True). The output metadata is in xml \n" +
|
|
"format only")
|
|
META = Component.Parameter('_meta',
|
|
public_name='meta',
|
|
default = 'xml',
|
|
type = str,
|
|
mandatory = False,
|
|
doc = "What type of metadata file is created. Possible values: xml or rsc (default: xml)")
|
|
SOURCE = Component.Parameter('_source',
|
|
public_name='source',
|
|
default = 1,
|
|
type = int,
|
|
mandatory = False,
|
|
doc = "DEM SRTM source. Possible values 1 or 3 (default: 1)")
|
|
NO_FILLING = Component.Parameter('_noFilling',
|
|
public_name='noFilling',
|
|
default = True,
|
|
type = bool,
|
|
mandatory = False,
|
|
doc = "If the flag is False the missing DEMs are filled with null values \n" + \
|
|
"(default: True, default null value -32768.")
|
|
FILLING_VALUE = Component.Parameter('_fillingValue',
|
|
public_name='fillingValue',
|
|
default = -32768,
|
|
type = int,
|
|
mandatory = False,
|
|
doc = "Value used to fill missing DEMs (default: -32768)")
|
|
BBOX = Component.Parameter('_bbox',
|
|
public_name='bbox',
|
|
default = None,
|
|
type = list,
|
|
mandatory = False,
|
|
doc = "Defines the spatial region in the format south north west east.\n" + \
|
|
"The values should be integers from (-90,90) for latitudes and (0,360) or " +\
|
|
"(-180,180) for longitudes.")
|
|
PAIRS = Component.Parameter('_pairs',
|
|
public_name='pairs',
|
|
default = None,
|
|
type = list,
|
|
mandatory = False,
|
|
doc = "Set of latitude and longitude pairs for which action = 'download' is performed.\n" +\
|
|
"The values should be integers from (-90,90)\n" + \
|
|
"for latitudes and (0,360) or (-180,180) for longitudes")
|
|
KEEP_DEMS = Component.Parameter('_keepDems',
|
|
public_name='keepDems',
|
|
default = False,
|
|
type = bool,
|
|
mandatory = False,
|
|
doc = "If the option is present then the single files used for stitching are kept.\n" + \
|
|
"If 'useLocalDirectory' is set then this flag is forced to True to avoid\n" + \
|
|
"accidental deletion of files (default: False)'")
|
|
REPORT = Component.Parameter('_report',
|
|
public_name='report',
|
|
default = False,
|
|
type = bool,
|
|
mandatory = False ,
|
|
doc = "If the option is present then failed and succeeded downloads are printed (default: False)")
|
|
USE_LOCAL_DIRECTORY = Component.Parameter('_useLocalDirectory',
|
|
public_name='useLocalDirectory',
|
|
default = False,
|
|
type = bool,
|
|
mandatory = False,
|
|
doc = "If the option is True then use the files that are in the location\n" + \
|
|
"specified by 'directory'. If not present 'directory' indicates\n" + \
|
|
"the directory where the files are downloaded (default: False)")
|
|
OUTPUT_FILE = Component.Parameter('_outputFile',
|
|
public_name='outputFile',
|
|
default = None,
|
|
type = str,
|
|
mandatory = False,
|
|
doc = "Name of the output file to be created in 'directory'.\n" + \
|
|
"If not provided the system generates one based on the bbox extremes")
|
|
|
|
REGIONS = Component.Parameter('_regions',
|
|
public_name='regions',
|
|
default = None,
|
|
type = list,
|
|
mandatory = False,
|
|
doc = "Regions where to look for the DEM files")
|
|
|
|
## This class provides a set of convenience method to retrieve and possibly combine different DEMs from the USGS server.
|
|
# \c NOTE: the latitudes and the longitudes that describe the DEMs refer to the bottom left corner of the image.
|
|
class DemStitcher(Component):
|
|
|
|
|
|
|
|
|
|
|
|
##
|
|
# Given a latitude and longitude in degrees it returns the expected filename.
|
|
# @param lat \c int latitude in the range (-90,90). Actual data are restricted to (-60,60) or so.
|
|
# @param lon \c int longitude in the range [-180,180) or [0,360).
|
|
# @return \c string the filename for that location
|
|
|
|
def createFilename(self,lat,lon,source = None):
|
|
|
|
if lon > 180:
|
|
lon = -(360 - lon)
|
|
else:
|
|
lon = lon
|
|
ns,ew = self.convertCoordinateToString(lat,lon)
|
|
return ns + ew + self._extension + self._zip
|
|
##
|
|
# Given a rectangle (in latitude,longitude space) defined by a maximum and minimum latitude and by a maximum and minimum longitude (in degrees) it returns
|
|
# an ordered list of the filenames defining the rectangle. The list is ordered first in ascending longitudes and teh ascending latitudes.
|
|
# @param lats \c list \c int list containing the minimum and maximum latitudes in the range (-90,90). Actual data are restricted to (-60,60) or so.
|
|
# @param lons \c list \c int list containing the minimum and maximum longitudes in the range [-180,180) or [0,360).
|
|
# @return \c tuple (\list strings the list of filenames covering the specified area, \c int the number of frames found along the longitude direction,
|
|
# \c int the number of frames found along the latitude direction)
|
|
|
|
def createNameList(self,lats,lons,source = None):
|
|
self._inputFileList = []
|
|
if lons[0] > 180:
|
|
lons[0] = -(360 - lons[0])
|
|
else:
|
|
lons[0] = lons[0]
|
|
if lons[1] > 180:
|
|
lons[1] = -(360 - lons[1])
|
|
else:
|
|
lons[1] = lons[1]
|
|
|
|
lonMin = min(lons[0],lons[1])
|
|
lons[1] = int(math.ceil(max(lons[0],lons[1])))
|
|
lons[0] = int(math.floor(lonMin))
|
|
#sanity check for lat
|
|
latMin = min(lats[0],lats[1])
|
|
lats[1] = int(math.ceil(max(lats[0],lats[1])))
|
|
lats[0] = int(math.floor(latMin))
|
|
# give error if crossing 180 and -180.
|
|
latList = []
|
|
lonList = []
|
|
for i in range(lats[0],lats[1]): # this leave out lats[1], but is ok because the last frame will go up to that point
|
|
latList.append(i)
|
|
#for lat go north to south
|
|
latList.reverse()
|
|
# create the list starting from the min to the max
|
|
if(lons[1] - lons[0] < 180):
|
|
for i in range(lons[0],lons[1]): # this leave out lons[1], but is ok because the last frame will go up to that point
|
|
lonList.append(i)
|
|
else:
|
|
print("Error. The crossing of E180 and W180 is not handled.")
|
|
raise Exception
|
|
self._latLonList = []
|
|
for lat in latList:
|
|
for lon in lonList:
|
|
name = self.createFilename(lat,lon,source)
|
|
self._inputFileList.append(name)
|
|
self._latLonList.append([lat,lon])
|
|
return self._inputFileList,len(latList),len(lonList)
|
|
|
|
##
|
|
# Given a rectangle (in latitude,longitude space) defined by a maximum and minimum
|
|
# latitude and by a maximum and minimum longitude (in degrees) it fetches
|
|
# the compressed (zip format) DEMs contained in that rectangle.
|
|
# @param lats \c list \c ints list containing the minimum and maximum latitudes
|
|
# in the range (-90,90). Actual data are restricted to (-60,60) or so.
|
|
# @param lons \c list \c ints list containing the minimum and maximum longitudes
|
|
# in the range [-180,180) or [0,360).
|
|
# @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data,
|
|
# source = 3 for 3 arcsec resolution data.
|
|
# @param downloadDir \c string the directory where the DEMs are downloaded.
|
|
# If the directory does not exists it will be created. If the argument is not provided
|
|
# then the files are downloaded in the location defined by the self._downloadDir
|
|
# that is defaulted to the current directory.
|
|
# @param region \c string region where to look for the files. If not provided the files
|
|
# are searched by scanning the content of each region. Use method getRagionList to get
|
|
# the list of possible region for a given source. Set region only if sure that all the
|
|
# requested file are contained in it.
|
|
def getDemsInBox(self,lats,lons,source,downloadDir = None,region = None):
|
|
nameList,numLat,numLon, = self.createNameList(lats,lons,source)
|
|
|
|
if downloadDir is None:
|
|
downloadDir = self._downloadDir
|
|
else:
|
|
self._downloadDir = downloadDir
|
|
#hackish. needs major refactoring. If self._useLocalDirectory is set we
|
|
#need only the nameList, no need to download
|
|
if not self._useLocalDirectory:
|
|
if region:
|
|
regionList = [region]*len(nameList)
|
|
else:
|
|
regionList = None
|
|
|
|
self.getDems(source,nameList,downloadDir,regionList)
|
|
|
|
else:
|
|
#create a fake download report from the nameList
|
|
files = os.listdir(downloadDir)
|
|
for fileNow in nameList:
|
|
#if file present then report success, failure otherwise
|
|
if files.count(fileNow):
|
|
self._downloadReport[fileNow] = self._succeded
|
|
else:
|
|
self._downloadReport[fileNow] = self._failed
|
|
|
|
return nameList,numLat,numLon
|
|
|
|
|
|
##
|
|
# Given a list of filenames it fetches the corresponding
|
|
# compressed (zip format) DEMs.
|
|
# @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data,
|
|
# source = 3 for 3 arcsec resolution data.
|
|
# @param listFile \c list of the filenames to be retrieved.
|
|
# @param downloadDir \c string the directory where the DEMs are downloaded.
|
|
# If the directory does not exists it will be created. If the argument is not
|
|
# provided then the files are downloaded in the location defined by the
|
|
# self._downloadDir that is defaulted to the current directory.
|
|
# @param region \c list \c strings regions where to look for the files. It must
|
|
# have the same length of \c listFile. If not provided the files are searched by
|
|
# scanning the content of each region. Use method getRegionList to get the list of
|
|
# possible regions for a given source. Set region only if sure that all the requested
|
|
# file are contained in it.
|
|
|
|
def getDems(self,source,listFile,downloadDir = None,region = None):
|
|
if downloadDir is None:
|
|
downloadDir = self._downloadDir
|
|
else:
|
|
self._downloadDir = downloadDir
|
|
|
|
if not (downloadDir) is None:
|
|
os.makedirs(downloadDir, exist_ok=True)
|
|
if region:
|
|
regionList = region
|
|
#region unknown, so try all of them
|
|
else:
|
|
# the scanning of the regions is time comsuming. get all the files in all region and create a big list
|
|
regionList = self.getRegionList(source)
|
|
|
|
regionMapping = []
|
|
fullList = []
|
|
for regionNow in regionList:
|
|
fileListUrl = self.getFileListPerRegion(source,regionNow)
|
|
if fileListUrl:
|
|
listNow = [file for file in fileListUrl]
|
|
fullList.extend(listNow)
|
|
regionNowMap = [regionNow]*len(fileListUrl)
|
|
regionMapping.extend(regionNowMap)
|
|
|
|
for fileNow in listFile:
|
|
url = ''
|
|
for i in range(len(fullList)):
|
|
if fileNow == fullList[i]:
|
|
regionNow = regionMapping[i]
|
|
url = self.getFullHttp(source,regionNow)
|
|
break
|
|
if not (url == ''):
|
|
try:
|
|
if not os.path.exists(os.path.join(downloadDir,fileNow)):
|
|
if(self._un is None or self._pw is None):
|
|
if os.path.exists(os.path.join(os.environ['HOME'],'.netrc')):
|
|
command = 'curl -n -L -c $HOME/.earthdatacookie -b $HOME/.earthdatacookie -k -f -O ' + os.path.join(url,fileNow)
|
|
else:
|
|
self.logger.error('Please create a .netrc file in your home directory containing\nmachine urs.earthdata.nasa.gov\n\tlogin yourusername\n\tpassword yourpassword')
|
|
sys.exit(1)
|
|
else:
|
|
command = 'curl -k -f -u ' + self._un + ':' + self._pw + ' -O ' + os.path.join(url,fileNow)
|
|
# curl with -O download in working dir, so save current, move to donwloadDir
|
|
# nd get back once download is finished
|
|
cwd = os.getcwd()
|
|
os.chdir(downloadDir)
|
|
if os.system(command):
|
|
os.chdir(cwd)
|
|
raise Exception
|
|
os.chdir(cwd)
|
|
self._downloadReport[fileNow] = self._succeded
|
|
except Exception as e:
|
|
self.logger.warning('There was a problem in retrieving the file %s. Exception %s'%(os.path.join(url,fileNow),str(e)))
|
|
self._downloadReport[fileNow] = self._failed
|
|
|
|
else:
|
|
self._downloadReport[fileNow] = self._failed
|
|
##
|
|
# After retriving DEMs this funtion prints the status of the download for each file, which could be 'succeded' or 'failed'
|
|
|
|
def printDownloadReport(self):
|
|
for k,v in self._downloadReport.items():
|
|
print('Download of file',k,v,'.')
|
|
##
|
|
# This function returns a dictionary whose keys are the attemped downloaded files and
|
|
# the values are the status of teh download, 'succeed' or 'failed'.
|
|
# @return \c dictionary whose keys are the attemped downloaded files and the values are
|
|
# the status of teh download, 'succeed' or 'failed'.
|
|
|
|
def getDownloadReport(self):
|
|
return self._downloadReport
|
|
|
|
##
|
|
# Given a list of latitudes and longitudes it fetches the corresponding
|
|
# compressed (zip format) DEMs.
|
|
# @param lats \c list \c int list containing set of latitudes in the range (-90,90).
|
|
# Actual data are restricted to (-60,60) or so.
|
|
# @param lons \c list \c int list containing set of longitudes in the range [-180,180)
|
|
# or [0,360).
|
|
# @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data,
|
|
# source = 3 for 3 arcsec resolution data.
|
|
# @param downloadDir \c string the directory where the DEMs are downloaded. If the
|
|
# directory does not exists it will be created. If the argument is not provided then
|
|
# the files are downloaded in the location defined by the self._downloadDir that is
|
|
# defaulted to the current directory.
|
|
# @param region \c list \c strings regions where to look for the files. It must have
|
|
# the same length of \c listFile. If not provided the files are searched by scanning
|
|
# the content of each region. Use method getRagionList to get the list of possible
|
|
# regions for a given source. Set region only if sure that all the requested file are
|
|
# contained in it.
|
|
|
|
def downloadFilesFromList(self,lats,lons,source,downloadDir = None,region = None):
|
|
|
|
inputFileList = []
|
|
for lat,lon in zip(lats,lons):
|
|
name = self.createFilename(lat,lon,source)
|
|
inputFileList.append(name)
|
|
self.getDems(source,inputFileList,downloadDir,region)
|
|
##
|
|
# Given a latitude and longitude it fetches the corresponding
|
|
# compressed (zip format) DEM.
|
|
# @param lat \c list \c int latitude in the range (-90,90). Actual data are restricted to (-60,60) or so.
|
|
# @param lons \c list \c int longitude in the range [-180,180) or [0,360).
|
|
# @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data, source = 3 for 3 arcsec resolution data.
|
|
# @param downloadDir \c string the directory where the DEMs are downloaded. If the directory does not exists it will be created. If the argument is not provided then the files are downloaded in the location defined by the self._downloadDir that is defaulted to the current directory.
|
|
# @param region \c list \c strings regions where to look for the files. It must have the same length of \c listFile. If not provided the files are searched by scanning the content of each region. Use method getRagionList to get the list of possible regions for a given source. Set region only if sure that all the requested file are contained in it.
|
|
def downloadFile(self,lat,lon,source,downloadDir = None,region = None):
|
|
name = self.createFilename(lat,lon,source)
|
|
inputFileList = [name]
|
|
regionList = [region]
|
|
self.getDems(source,inputFileList,downloadDir,regionList)
|
|
##
|
|
# It returns the list of DEMs for a give source and region (if provided). If the region is not provided the full list of files for that source type is returned.
|
|
# @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data, source = 3 for 3 arcsec resolution data.
|
|
# @param region \c list \c strings regions where to look for the files. If the region is not provided the full list of files for that source type is returned.
|
|
# @return \c list \c string list containing the the filenames found for the specific source and (if specified) region.
|
|
|
|
def getFileList(self,source,region = None):
|
|
retList = []
|
|
if region:
|
|
regionList = self.getRegionList(source)
|
|
foundRegion = False
|
|
for el in regionList:
|
|
if el == region:
|
|
foundRegion = True
|
|
if foundRegion:
|
|
retList = self.getFileListPerRegion(source,region)
|
|
else:
|
|
regionList = self.getRegionList(source)
|
|
for el in regionList:
|
|
retList.extend(self.getFileListPerRegion(source,el))
|
|
return retList
|
|
|
|
##
|
|
# It returns the list of DEMs for a given source and region.
|
|
# @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data, source = 3 for 3 arcsec resolution data.
|
|
# @param region \c list \c strings regions where to look for the files.
|
|
# @return \c list \c string list containing the the filenames found for the specific source and region.
|
|
def getFileListPerRegion(self,source,region):
|
|
url = self.getFullHttp(source,region)
|
|
return self.getUrlList(url,self._filters['fileExtension'], self._remove)
|
|
|
|
|
|
##
|
|
# It returns the list of regions for a given source.
|
|
# @param source \c int the type of DEM. source = 1 for 1 arcsec resolution data, source = 3 for 3 arcsec resolution data.
|
|
# @return \c list \c string list of region for the specified source.
|
|
def getRegionList(self,source):
|
|
# check first if it has been computed before
|
|
if self._regionList[str(source)] == []:
|
|
url = self.http + str(source)
|
|
self._regionList[str(source)] = self.getUrlList(url,self._filters['region'+str(source)], self._remove)
|
|
return self._regionList[str(source)]
|
|
|
|
def getUrlList(self,url,filterList = None, removeList = None):
|
|
if filterList is None:
|
|
filterList = []
|
|
if removeList is None:
|
|
removeList = []
|
|
if self._un is None or self._pw is None:
|
|
fp = urllib.request.urlopen(url)
|
|
allUrl = fp.read()
|
|
fp.close()
|
|
else:
|
|
# create a password manager
|
|
password_mgr = urllib.request.HTTPPasswordMgrWithDefaultRealm()
|
|
# Add the username and password.
|
|
password_mgr.add_password(None,url,self._un,self._pw)
|
|
handler = urllib.request.HTTPBasicAuthHandler(password_mgr)
|
|
# create "opener" (OpenerDirector instance)
|
|
opener = urllib.request.build_opener(handler)
|
|
# use the opener to fetch a URL
|
|
allUrl = opener.open(url).read()
|
|
|
|
ddp = DemDirParser()
|
|
# feed the data from the read() of the url to the parser. It will call the DemDirParser.handle_data everytime t
|
|
# a data type is parsed
|
|
ddp.filterList = filterList
|
|
ddp.removeList = removeList
|
|
ddp.feed(allUrl.decode('utf-8', 'replace'))
|
|
return ddp.results
|
|
|
|
##
|
|
# Setter function for the download directory.
|
|
# @param ddir \c string directory where the DEMs are downloaded. In self.stitchDem defines also the directory where the output stiched file is saved.
|
|
|
|
def setDownloadDirectory(self,ddir):
|
|
self._downloadDir = ddir
|
|
|
|
##
|
|
# Fuction that decompress the given file in zip format.
|
|
# @param filename \c strig the name of the file to decompress.
|
|
# @param downloadDir \c string the directory where the DEMs are downloaded. If the directory does not exists it will be created. If the argument is not provided then the files are downloaded in the location defined by the self._downloadDir that is defaulted to the current directory.
|
|
def decompress(self,filename,downloadDir = None,keep = None):
|
|
|
|
# keep .zip by default
|
|
if keep == None:
|
|
keep = True
|
|
if downloadDir is None:
|
|
downloadDir = self._downloadDir
|
|
else:
|
|
self._downloadDir = downloadDir
|
|
|
|
filen = os.path.join(downloadDir,filename)
|
|
try:
|
|
#some system might not have zlib so you a system call to unzip
|
|
zip = zipfile.ZipFile(filen,'r')
|
|
import zlib
|
|
zip.extractall(downloadDir)
|
|
except:
|
|
self.extract(downloadDir,filen)
|
|
|
|
if not keep:
|
|
os.remove(filen)
|
|
|
|
def extract(self,downloadDir,filen):
|
|
os.system('unzip -o -qq ' + os.path.join(filen) + ' -d ' + downloadDir)
|
|
|
|
|
|
def defaultName(self,snwe):
|
|
latMin = math.floor(snwe[0])
|
|
latMax = math.ceil(snwe[1])
|
|
lonMin = math.floor(snwe[2])
|
|
lonMax = math.ceil(snwe[3])
|
|
nsMin,ewMin = self.convertCoordinateToString(latMin, lonMin)
|
|
nsMax,ewMax = self.convertCoordinateToString(latMax, lonMax)
|
|
demName = (
|
|
'demLat_' + nsMin + '_' +nsMax +
|
|
'_Lon_' + ewMin +
|
|
'_' + ewMax + '.dem'
|
|
)
|
|
|
|
return demName
|
|
|
|
def convertCoordinateToString(self,lat,lon):
|
|
|
|
if(lon > 180):
|
|
lon = -(360 - lon)
|
|
if(lon < 0):
|
|
ew = 'W'
|
|
else:
|
|
ew = 'E'
|
|
lonAbs = int(math.fabs(lon))
|
|
if(lonAbs >= 100):
|
|
ew += str(lonAbs)
|
|
elif(lonAbs < 10):
|
|
ew += '00' + str(lonAbs)
|
|
else:
|
|
ew += '0' + str(lonAbs)
|
|
|
|
if(int(lat) >= 0):
|
|
ns = 'N'
|
|
else:
|
|
ns = 'S'
|
|
latAbs = int(math.fabs(lat))
|
|
if(latAbs >= 10):
|
|
ns += str(latAbs)
|
|
else:
|
|
ns += '0' +str(latAbs)
|
|
|
|
return ns,ew
|
|
|
|
|
|
#based on the source predict the width of the dem
|
|
def getDemWidth(self,lon,source):
|
|
if source == 3:
|
|
factor = 1200
|
|
else:
|
|
factor = 3600
|
|
return int(math.fabs((lon[1] - lon[0]))*factor)
|
|
|
|
#this method also create an actual DeimImage object that is returned by the getImage() method
|
|
def createXmlMetadata(self,lat,lon,source,outname):
|
|
|
|
demImage = self.createImage(lat,lon,source,outname)
|
|
demImage.renderHdr()
|
|
|
|
def createImage(self,lat,lon,source,outname):
|
|
from isceobj.Image import createDemImage
|
|
|
|
demImage = createDemImage()
|
|
if source == 3:
|
|
delta = 1/1200.0
|
|
else:
|
|
delta = 1/3600.0
|
|
|
|
os.makedirs(self._downloadDir, exist_ok=True)
|
|
|
|
width = self.getDemWidth(lon,source)
|
|
demImage.initImage(outname,'read',width)
|
|
length = demImage.getLength()
|
|
dictProp = {'METADATA_LOCATION':outname+'.xml','REFERENCE':self._reference,'Coordinate1':{'size':width,'startingValue':min(lon[0],lon[1]),'delta':delta},'Coordinate2':{'size':length,'startingValue':max(lat[0],lat[1]),'delta':-delta},'FILE_NAME':outname}
|
|
#no need to pass the dictionaryOfFacilities since init will use the default one
|
|
demImage.init(dictProp)
|
|
self._image = demImage
|
|
return demImage
|
|
|
|
##
|
|
#Function to indent an element of an ElementTree object. If the element passed is the root element, then all the ElementTree object is indented.
|
|
#@param elem element of an ElementTree object.
|
|
|
|
def indent(self,elem, depth = None,last = None):
|
|
if depth == None:
|
|
depth = [0]
|
|
if last == None:
|
|
last = False
|
|
tab = ' '*4
|
|
if(len(elem)):
|
|
depth[0] += 1
|
|
elem.text = '\n' + (depth[0])*tab
|
|
lenEl = len(elem)
|
|
lastCp = False
|
|
for i in range(lenEl):
|
|
if(i == lenEl - 1):
|
|
lastCp = True
|
|
self.indent(elem[i],depth,lastCp)
|
|
if(not last):
|
|
elem.tail = '\n' + (depth[0])*tab
|
|
else:
|
|
depth[0] -= 1
|
|
elem.tail = '\n' + (depth[0])*tab
|
|
else:
|
|
if(not last):
|
|
elem.tail = '\n' + (depth[0])*tab
|
|
else:
|
|
depth[0] -= 1
|
|
elem.tail = '\n' + (depth[0])*tab
|
|
|
|
def writeFileFromDictionary(self,file,dict, name = None):
|
|
if not name:
|
|
name = ''
|
|
root = ET.Element('component')
|
|
nameSubEl = ET.SubElement(root,'name')
|
|
nameSubEl.text = name
|
|
for key, val in dict.items():
|
|
propSubEl = ET.SubElement(root,'property')
|
|
ET.SubElement(propSubEl, 'name').text = key
|
|
ET.SubElement(propSubEl, 'value').text = str(val)
|
|
|
|
|
|
self.indent(root)
|
|
etObj = ET.ElementTree(root)
|
|
etObj.write(file)
|
|
def createRscMetadata(self,lat,lon,source,outname):
|
|
|
|
demImage = self.createImage(lat,lon,source,outname)
|
|
|
|
dict = {'WIDTH':demImage.width,'LENGTH':demImage.length,'X_FIRST':demImage.coord1.coordStart,'Y_FIRST':demImage.coord2.coordStart,'X_STEP':demImage.coord1.coordDelta,'Y_STEP':-demImage.coord2.coordDelta,'X_UNIT':'degrees','Y_UNIT':'degrees'}
|
|
os.makedirs(self._downloadDir, exist_ok=True)
|
|
extension = '.rsc'
|
|
outfile = outname + extension
|
|
fp = open(outfile,'w')
|
|
for k,v in dict.items():
|
|
fp.write(str(k) + '\t' + str(v) + '\n')
|
|
fp.close()
|
|
|
|
def setKeepDems(self,val):
|
|
self._keepDems = val
|
|
|
|
def setCreateXmlMetadata(self,val):
|
|
self._createXmlMetadata = val
|
|
|
|
def setCreateRscMetadata(self,val):
|
|
self._createRscMetadata = val
|
|
|
|
def setMetadataFilename(self,demName):
|
|
self._metadataFilename = demName
|
|
|
|
def setFillingFilename(self,name):
|
|
self._fillingFilename = name
|
|
|
|
def setFillingValue(self,val):
|
|
self._fillingValue = val
|
|
|
|
def setFilling(self):
|
|
self._noFilling = False
|
|
|
|
def setNoFilling(self):
|
|
self._noFilling = True
|
|
|
|
def setUseLocalDirectory(self,val):
|
|
self._useLocalDirectory = val
|
|
def getUrl(self):
|
|
return self._url
|
|
def setUrl(self,url):
|
|
self._url = url
|
|
#after the url has been set generate the full path
|
|
self._http = self._url + '/srtm/version2_1/SRTM'
|
|
|
|
def setUsername(self,un):
|
|
self._un = un
|
|
|
|
def setPassword(self,pw):
|
|
self._pw = pw
|
|
|
|
def createFillingTile(self,source,swap,filename):
|
|
fp = open(filename,'wb')
|
|
numSamples = 1201
|
|
if (source == 1):
|
|
numSamples = 3601
|
|
|
|
if swap:
|
|
# pack it as a big endian and unpack it, and get the swapped number
|
|
fillingValue = struct.unpack('h',struct.pack('>h',self._fillingValue))[0]
|
|
else:
|
|
fillingValue = self._fillingValue
|
|
fullTile = [fillingValue]*numSamples*numSamples
|
|
tile = array('h')
|
|
tile.fromlist(fullTile)
|
|
tile.tofile(fp)
|
|
fp.close()
|
|
|
|
|
|
#allow to overwrite from subclasses the nameing convention of the unzipped
|
|
def getUnzippedName(self,name,source = None):
|
|
return name.replace(self._zip,'')
|
|
def stitchDems(self,lat,lon,source, outname, downloadDir = None,region = None, keep = None, swap = None):
|
|
import glob
|
|
if downloadDir is None:
|
|
downloadDir = self._downloadDir
|
|
else:
|
|
self._downloadDir = downloadDir
|
|
|
|
swapFlag = 0
|
|
if swap:
|
|
if swap == True:#might be true or false
|
|
swapFlag = 1
|
|
else: # do it by default
|
|
swapFlag = 1
|
|
|
|
|
|
listNames,nLat,nLon = self.getDemsInBox(lat,lon,source,downloadDir,region)
|
|
unzip = True
|
|
#keep track of the synthetic ones since they don't need to be unzipped
|
|
syntheticTiles = []
|
|
if self._noFilling:
|
|
#make sure that we have all the file to cover the region. check if some download failed
|
|
for k,v in self._downloadReport.items():
|
|
if v == self._failed:
|
|
unzip = False
|
|
#clean up the dowloaded files if it failed since when trying a second source it might endup
|
|
#stitching them together beacaiuse it does not re-download the ones present and unfortunately
|
|
#the dems with different resolution have the same name convention
|
|
if not self._keepAfterFailed:
|
|
os.system("rm -rf " + downloadDir + "/*.hgt*")
|
|
break
|
|
else:
|
|
syntTileCreated = False
|
|
#check and send a warning if the full region is not available
|
|
if not self._succeded in self._downloadReport.values():
|
|
self.logger.warning('The full region of interested is not available. A DEM with all null values will be created.')
|
|
for k,v in self._downloadReport.items():
|
|
if v == self._failed:#symlink each missing file to the reference one created in createFillingFile
|
|
if not syntTileCreated:#create the synthetic Tile the first time around
|
|
#get the abs path otherwise the symlink doesn't work
|
|
tileName = os.path.abspath(os.path.join(downloadDir,self._fillingFilename))
|
|
self.createFillingTile(source,swapFlag,tileName)
|
|
syntTileCreated = True
|
|
|
|
syntheticTiles.append(k)
|
|
demName = os.path.join(downloadDir,self.getUnzippedName(k,source))
|
|
#check for lexists so it returns also broken links, just in case something went wrong before
|
|
if os.path.lexists(demName):#clean up to make sure that old names are not there. will cause problem if use old one and the resolution od the dem is changed
|
|
os.remove(demName)
|
|
os.symlink(tileName,demName)
|
|
|
|
if unzip:
|
|
decompressedList = []
|
|
for name in listNames:
|
|
if not name in syntheticTiles:#synthetic tiles don't need to be decompressed
|
|
self.decompress(name,downloadDir,keep)
|
|
|
|
newName = self.getUnzippedName(name,source)
|
|
if downloadDir:
|
|
newName = os.path.join(downloadDir,newName)
|
|
|
|
decompressedList.append(bytes(newName, 'utf-8'))
|
|
numSamples = 1201
|
|
if (source == 1):
|
|
numSamples = 3601
|
|
|
|
outname = os.path.join(downloadDir,outname)
|
|
numFiles = [nLat,nLon]
|
|
fileListIn_c = (c_char_p * len(decompressedList))()
|
|
fileListIn_c[:] = decompressedList
|
|
numFiles_c = (c_int * len(numFiles))()
|
|
numFiles_c[:] = numFiles
|
|
fileOut_c = c_char_p(bytes(outname, 'utf-8'))
|
|
numSamples_c = c_int(numSamples)
|
|
swapFlag_c = c_int(swapFlag)
|
|
self._lib.concatenateDem(fileListIn_c,numFiles_c,fileOut_c,byref(numSamples_c),byref(swapFlag_c))
|
|
|
|
if not self._keepDems:
|
|
for dem in decompressedList:
|
|
for d in glob.glob('*'+os.path.basename(dem.decode('UTF-8'))[:7]+'*'):
|
|
os.remove(d)
|
|
if self._createXmlMetadata:
|
|
self.createXmlMetadata(lat,lon,source,outname)
|
|
if self._createRscMetadata:
|
|
self.createRscMetadata(lat,lon,source,outname)
|
|
|
|
return unzip #if False it means that failed
|
|
|
|
## Corrects the self._image from EGM96 to WGS84 and viceversa.
|
|
#@param image \c Image if provided is used instead of the instance attribute self._image
|
|
#@param conversionType \c int -1 converts from EGM96 to WGS84, 1 converts from WGS84 to EGM96
|
|
#@return \c Image instance the converted Image
|
|
def correct(self,image = None,conversionType=-1):
|
|
'''Corrects the self._image from EGM96 to WGS84 and viceversa.'''
|
|
from contrib.demUtils.Correct_geoid_i2_srtm import (
|
|
Correct_geoid_i2_srtm
|
|
)
|
|
cg = Correct_geoid_i2_srtm()
|
|
return cg(image,conversionType) if image else cg(self._image,conversionType)
|
|
|
|
#still need to call it since the initialization calls the _url so the setter of
|
|
#url does not get called
|
|
def _configure(self):
|
|
#after the url has been set generate the full path
|
|
self._http = self._url + '/srtm/version2_1/SRTM'
|
|
|
|
def getImage(self):
|
|
return self._image
|
|
|
|
def __getstate__(self):
|
|
d = dict(self.__dict__)
|
|
del d['logger']
|
|
del d['_lib']
|
|
return d
|
|
|
|
def __setstate__(self,d):
|
|
self.__dict__.update(d)
|
|
self.logger = logging.getLogger('isce.contrib.demUtils.DemStitcher')
|
|
libName = os.path.join(os.path.dirname(__file__),self._loadLibName)
|
|
##self._keepAfterFailed = False #if True keeps the downloaded files even if the stitching failed.
|
|
self._lib = cdll.LoadLibrary(libName)
|
|
return
|
|
|
|
def main(self):
|
|
# prevent from deliting local files
|
|
if(self._useLocalDirectory):
|
|
self._keepAfterFailed = True
|
|
self._keepDems = True
|
|
# is a metadata file is created set the right type
|
|
if(self._meta == 'xml'):
|
|
self.setCreateXmlMetadata(True)
|
|
elif(self._meta == 'rsc'):
|
|
self.setCreateRscMetadata(True)
|
|
# check for the action to be performed
|
|
if(self._action == 'stitch'):
|
|
if(self._bbox):
|
|
lat = self._bbox[0:2]
|
|
lon = self._bbox[2:4]
|
|
if (self._outputFile is None):
|
|
self._outputFile = self.defaultName(self._bbox)
|
|
|
|
if not(self.stitchDems(lat,lon,self._source,self._outputFile,self._downloadDir, \
|
|
keep=self._keepDems)):
|
|
print('Could not create a stitched DEM. Some tiles are missing')
|
|
else:
|
|
if(self._correct):
|
|
width = self.getDemWidth(lon,self._source)
|
|
self.correct()
|
|
#self.correct(self._output,self._source,width,min(lat[0],lat[1]),min(lon[0],lon[1]))
|
|
else:
|
|
print('Error. The "bbox" attribute must be specified when the action is "stitch"')
|
|
raise ValueError
|
|
elif(self._action == 'download'):
|
|
if(self._bbox):
|
|
lat = self._bbox[0:2]
|
|
lon = self._bbox[2:4]
|
|
self.getDemsInBox(lat,lon,self._source,self._downloadDir)
|
|
#can make the bbox and pairs mutually esclusive if replace the if below with elif
|
|
if(self._pairs):
|
|
self.downloadFilesFromList(self._pairs[::2],self._pairs[1::2],self._source,self._downloadDir)
|
|
if(not (self._bbox or self._pairs)):
|
|
print('Error. Either the "bbox" attribute or the "pairs" attribute must be specified when --action download is used')
|
|
raise ValueError
|
|
|
|
else:
|
|
print('Unrecognized action ',self._action)
|
|
return
|
|
|
|
if(self._report):
|
|
for k,v in self._downloadReport.items():
|
|
print(k,'=',v)
|
|
|
|
|
|
def _facilities(self):
|
|
super(DemStitcher,self)._facilities()
|
|
def getFullHttp(self,source, region = None):
|
|
toAppend = ''
|
|
if region:
|
|
toAppend = ('/' + region + '/')
|
|
return self._http + str(source) + toAppend
|
|
@property
|
|
def http(self):
|
|
return self._http
|
|
family = 'demstitcher'
|
|
parameter_list = (
|
|
URL,
|
|
USERNAME,
|
|
PASSWORD,
|
|
KEEP_AFTER_FAILED,
|
|
DIRECTORY,
|
|
ACTION,
|
|
CORRECT,
|
|
META,
|
|
SOURCE,
|
|
NO_FILLING,
|
|
FILLING_VALUE,
|
|
BBOX,
|
|
PAIRS,
|
|
KEEP_DEMS,
|
|
REPORT,
|
|
USE_LOCAL_DIRECTORY,
|
|
OUTPUT_FILE,
|
|
REGIONS
|
|
)
|
|
def __init__(self,family = '', name = ''):
|
|
|
|
self._loadLibName = "demStitch.so"
|
|
libName = os.path.join(os.path.dirname(__file__),self._loadLibName)
|
|
##self._keepAfterFailed = False #if True keeps the downloaded files even if the stitching failed.
|
|
self._lib = cdll.LoadLibrary(libName)
|
|
self._downloadReport = {}
|
|
# Note if _useLocalDirectory is True then the donwloadDir is the local directory
|
|
##self._downloadDir = os.getcwd()#default to the cwd
|
|
self._inputFileList = []
|
|
##self._useLocalDirectory = False
|
|
##self._outputFile = ''
|
|
##self._un = un
|
|
##self._pw = pw
|
|
self._extension = '.hgt'
|
|
self._zip = '.zip'
|
|
|
|
#to make it working with other urls, make sure that the second part of the url
|
|
#it's /srtm/version2_1/SRTM(1,3)
|
|
self._filters = {'region1':['Region'],'region3':['Africa','Australia','Eurasia','Islands','America'],'fileExtension':['.hgt.zip']}
|
|
self._remove = ['.jpg']
|
|
self._metadataFilename = 'fileDem.dem'
|
|
self._createXmlMetadata = None
|
|
self._createRscMetadata = None
|
|
self._regionList = {'1':[],'3':[]}
|
|
##self._keepDems = False
|
|
self._fillingFilename = 'filling.hgt' # synthetic tile to cover holes
|
|
##self._fillingValue = -32768 # fill the synthetic tile with this value
|
|
##self._noFilling = False
|
|
self._failed = 'failed'
|
|
self._succeded = 'succeded'
|
|
self._image = None
|
|
self._reference = 'EGM96'
|
|
super(DemStitcher, self).__init__(family if family else self.__class__.family, name=name)
|
|
# logger not defined until baseclass is called
|
|
|
|
if not self.logger:
|
|
self.logger = logging.getLogger('isce.contrib.demUtils.DemStitcher')
|
|
|
|
url = property(getUrl,setUrl)
|