217 lines
8.7 KiB
Python
217 lines
8.7 KiB
Python
|
#!/usr/bin/env python3
|
|||
|
|
|
|||
|
|
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|||
|
|
# LAMP License
|
|||
|
|
#
|
|||
|
|
# Author: chenzenghui
|
|||
|
|
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|||
|
|
|
|||
|
|
# 自定义dem管理
|
|||
|
|
# 1. 创建二进制文件 ENVI hdr
|
|||
|
|
# 2. 读取文件构建.vrt .xml
|
|||
|
|
|
|||
|
|
import argparse
|
|||
|
|
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
|
|||
|
|
import numpy as np
|
|||
|
|
from isce import logging
|
|||
|
|
from iscesys.Component.Component import Component
|
|||
|
|
from isceobj.Image import createDemImage
|
|||
|
|
from osgeo import gdal,osr,ogr
|
|||
|
|
import xml.etree.ElementTree as ET
|
|||
|
|
from html.parser import HTMLParser
|
|||
|
|
import time
|
|||
|
|
env_str = os.path.dirname(os.path.abspath(sys.argv[0]))
|
|||
|
|
os.environ['PROJ_LIB'] = env_str
|
|||
|
|
|
|||
|
|
|
|||
|
|
class DEM2ISCE(Component):
|
|||
|
|
|
|||
|
|
def dem_merged(self,in_dem_path, out_dem_path):
|
|||
|
|
'''
|
|||
|
|
DEM重采样函数，默认坐标系为WGS84
|
|||
|
|
agrs:
|
|||
|
|
in_dem_path: 输入的DEM文件夹路径
|
|||
|
|
meta_file_path: 输入的xml元文件路径
|
|||
|
|
out_dem_path: 输出的DEM文件夹路径
|
|||
|
|
'''
|
|||
|
|
# 读取文件夹中所有的DEM
|
|||
|
|
# dem_file_paths = in_dem_path
|
|||
|
|
lons_min = []
|
|||
|
|
lons_max = []
|
|||
|
|
lats_min = []
|
|||
|
|
lats_max = []
|
|||
|
|
dem_file_paths=[os.path.join(in_dem_path,dem_name) for dem_name in os.listdir(in_dem_path) if (dem_name.find(".tif")>=0 or dem_name.find(".tiff")>=0) and dem_name.find(".tif.")==-1]
|
|||
|
|
spatialreference=osr.SpatialReference()
|
|||
|
|
spatialreference.SetWellKnownGeogCS("WGS84") # 设置地理坐标，单位为度 degree # 设置投影坐标，单位为度 degree
|
|||
|
|
spatialproj=spatialreference.ExportToWkt() # 导出投影结果
|
|||
|
|
# 将DEM拼接成一张大图
|
|||
|
|
for file in dem_file_paths:
|
|||
|
|
dataset = gdal.Open(file)
|
|||
|
|
transform = dataset.GetGeoTransform()
|
|||
|
|
im_width = dataset.RasterXSize
|
|||
|
|
im_height = dataset.RasterYSize
|
|||
|
|
lon_right = transform[0] + transform[1] * im_width
|
|||
|
|
lat_down = transform[3] + transform[5] * im_height
|
|||
|
|
lons_min.append(float(transform[0]))
|
|||
|
|
lons_max.append(float(lon_right))
|
|||
|
|
lats_min.append(float(lat_down))
|
|||
|
|
lats_max.append(float(transform[3]))
|
|||
|
|
lon_min = round(np.min(lons_min))
|
|||
|
|
lon_max = round(np.max(lons_max))
|
|||
|
|
lat_min = round(np.min(lats_min))
|
|||
|
|
lat_max = round(np.max(lats_max))
|
|||
|
|
|
|||
|
|
|
|||
|
|
all_file_name = 'demLat_N' + str(lat_min) + '_N' + str(lat_max) + '_Lon_E' + str(lon_min) + '_E' + str(lon_max) + '.tiff'
|
|||
|
|
# print("lon_min :{} lon_max :{}, lat_min :{}, lat_max: {}".format(lon_min, lon_max, lat_min, lat_max))
|
|||
|
|
# print(file_name)
|
|||
|
|
out_DEM_path=os.path.join(out_dem_path, all_file_name)
|
|||
|
|
# out_DEM=out_dem_path
|
|||
|
|
gdal.Warp(out_DEM_path,
|
|||
|
|
dem_file_paths,
|
|||
|
|
format="GTiff",
|
|||
|
|
dstSRS=spatialproj,
|
|||
|
|
dstNodata=self._NoDataValue,
|
|||
|
|
outputType=gdal.GDT_Float32)
|
|||
|
|
|
|||
|
|
file_name_vrt = 'demLat_N' + str(lat_min) + '_N' + str(lat_max) + '_Lon_E' + str(lon_min) + '_E' + str(lon_max) + '.dem.wgs84.vrt'
|
|||
|
|
file_name = 'demLat_N' + str(lat_min) + '_N' + str(lat_max) + '_Lon_E' + str(lon_min) + '_E' + str(lon_max) + '.dem.wgs84'
|
|||
|
|
# print("lon_min :{} lon_max :{}, lat_min :{}, lat_max: {}".format(lon_min, lon_max, lat_min, lat_max))
|
|||
|
|
# print(file_name)
|
|||
|
|
|
|||
|
|
mergeFile =gdal.BuildVRT(os.path.join(out_dem_path,file_name_vrt),out_DEM_path)
|
|||
|
|
out_DEM=os.path.join(out_dem_path, file_name)
|
|||
|
|
# out_DEM=out_dem_path
|
|||
|
|
gdal.Warp(out_DEM,
|
|||
|
|
mergeFile,
|
|||
|
|
format="ENVI",
|
|||
|
|
dstSRS=spatialproj,
|
|||
|
|
dstNodata=self._NoDataValue,
|
|||
|
|
outputType=gdal.GDT_Float32)
|
|||
|
|
|
|||
|
|
time.sleep(3)
|
|||
|
|
return out_DEM
|
|||
|
|
|
|||
|
|
#this method also create an actual DeimImage object that is returned by the getImage() method
|
|||
|
|
def createXmlMetadata(self,outname):
|
|||
|
|
demImage = self.createImage(outname)
|
|||
|
|
demImage.renderHdr()
|
|||
|
|
|
|||
|
|
def getDemWidth(self,outname):
|
|||
|
|
gdal.AllRegister()
|
|||
|
|
dataset=gdal.Open(outname)
|
|||
|
|
width=dataset.RasterXSize
|
|||
|
|
del dataset
|
|||
|
|
return width
|
|||
|
|
|
|||
|
|
def getDemHeight(self,outname):
|
|||
|
|
gdal.AllRegister()
|
|||
|
|
dataset=gdal.Open(outname)
|
|||
|
|
height=dataset.RasterYSize
|
|||
|
|
del dataset
|
|||
|
|
return height
|
|||
|
|
|
|||
|
|
def getGeotransform(self,outname):
|
|||
|
|
gdal.AllRegister()
|
|||
|
|
dataset=gdal.Open(outname)
|
|||
|
|
geotransform = dataset.GetGeoTransform()
|
|||
|
|
del dataset
|
|||
|
|
return geotransform
|
|||
|
|
|
|||
|
|
def createImage(self,outname):
|
|||
|
|
demImage = createDemImage()
|
|||
|
|
width = self.getDemWidth(outname)
|
|||
|
|
height=self.getDemHeight(outname)
|
|||
|
|
demImage.initImage(outname,'read',width,type="float")
|
|||
|
|
length = demImage.getLength()
|
|||
|
|
# 获取分辨率
|
|||
|
|
geotransform=self.getGeotransform(outname)
|
|||
|
|
|
|||
|
|
dictProp = {'METADATA_LOCATION':outname+'.xml','REFERENCE':self._reference,'Coordinate1':{'size':width,'startingValue':geotransform[0],'delta':geotransform[1]},'Coordinate2':{'size':length,'startingValue':geotransform[3],'delta':geotransform[5]},'FILE_NAME':outname}
|
|||
|
|
#no need to pass the dictionaryOfFacilities since init will use the default one
|
|||
|
|
demImage.init(dictProp)
|
|||
|
|
self._image = demImage
|
|||
|
|
return demImage
|
|||
|
|
|
|||
|
|
def setFillingValue(self,val):
|
|||
|
|
self._fillingValue = val
|
|||
|
|
|
|||
|
|
def setNoDataValue(self,val):
|
|||
|
|
self._NoDataValue = val
|
|||
|
|
|
|||
|
|
|
|||
|
|
def stitchDems(self,source, outname):
|
|||
|
|
import glob
|
|||
|
|
# 合并数据
|
|||
|
|
out_dem = self.dem_merged(source, outname)
|
|||
|
|
self.createXmlMetadata(out_dem)
|
|||
|
|
family = 'DEM2ISCE'
|
|||
|
|
def __init__(self,family = '', name = ''):
|
|||
|
|
self._extension = '.tif'
|
|||
|
|
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(DEM2ISCE, 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.DEM2ISCE')
|
|||
|
|
|
|||
|
|
def getImage(self):
|
|||
|
|
return self._image
|
|||
|
|
|
|||
|
|
|
|||
|
|
# DEM转换主流程
|
|||
|
|
def processDEM2ISCE(name,source_path,target_path,fillvalue,noDataValue):
|
|||
|
|
ds = DEM2ISCE(name=name)
|
|||
|
|
# 构建
|
|||
|
|
ds.setFillingValue(fillvalue)
|
|||
|
|
ds.setNoDataValue(noDataValue)
|
|||
|
|
ds.stitchDems(source_path,target_path)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def main():
|
|||
|
|
#if not argument provided force the --help flag
|
|||
|
|
if(len(sys.argv) == 1):
|
|||
|
|
sys.argv.append('-h')
|
|||
|
|
# Use the epilog to add usage examples
|
|||
|
|
epilog = '将格式为tif 的DEM 转换为ISCE 支持的DEM格式:\n\n'
|
|||
|
|
epilog += 'Usage examples:\n\n'
|
|||
|
|
epilog += 'DEM2ISCE.py -s /mnt/d/codestorage/isce2/青海省.tif -o /mnt/d/codestorage/isce2/青海省_wgs84 -fillvalue -9999 -Nodata -9999\n\n'
|
|||
|
|
|
|||
|
|
#set the formatter_class=argparse.RawDescriptionHelpFormatter otherwise it splits the epilog lines with its own default format
|
|||
|
|
parser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter,epilog=epilog)
|
|||
|
|
parser.add_argument('-s', '--source', type = str, default ="/mnt/d/codestorage/isce2/青海省.tif", dest = 'source_path', help = '输入dem,格式为tif')
|
|||
|
|
parser.add_argument('-o', '--outpath', type = str, default = '/mnt/d/codestorage/isce2/青海省_wgs84', dest = 'outpath', help = '输出isce 支持的DEM ')
|
|||
|
|
parser.add_argument('-fillvalue', '--fillvalue', type = float, default = -9999, dest = 'fillvalue', help = '空值填充')
|
|||
|
|
parser.add_argument('-Nodata', '--Nodata', type = float, default = -9999, dest = 'Nodatavalue', help = '无效值填充')
|
|||
|
|
args = parser.parse_args()
|
|||
|
|
processDEM2ISCE("DEM2ISCE",args.source_path,args.outpath,args.fillvalue,args.Nodatavalue)
|
|||
|
|
print("DEM==>ISCE ok")
|
|||
|
|
return -1
|
|||
|
|
|
|||
|
|
if __name__ == '__main__':
|
|||
|
|
main()
|