|
#!/usr/bin/env python3 |
|
# name : fgddem.py |
|
# purpose : translates digital elevation model of Fundamental Geospatial Data provided by GSI |
|
# into a GDAL supported format (which has write access). |
|
# license : The MIT License |
|
# copyright : (c) 2012, Minoru Akagi |
|
# |
|
# Permission is hereby granted, free of charge, to any person obtaining a |
|
# copy of this software and associated documentation files (the "Software"), |
|
# to deal in the Software without restriction, including without limitation |
|
# the rights to use, copy, modify, merge, publish, distribute, sublicense, |
|
# and/or sell copies of the Software, and to permit persons to whom the |
|
# Software is furnished to do so, subject to the following conditions: |
|
# |
|
# The above copyright notice and this permission notice shall be included |
|
# in all copies or substantial portions of the Software. |
|
# |
|
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
|
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
|
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
|
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
|
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
|
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
|
# DEALINGS IN THE SOFTWARE. |
|
|
|
script_version = "0.9" |
|
|
|
import sys, os |
|
import datetime |
|
import glob |
|
import numpy |
|
import re |
|
import shutil |
|
from xml.dom import minidom |
|
import zipfile |
|
from osgeo import gdal |
|
|
|
try: |
|
progress = gdal.TermProgress_nocb |
|
except: |
|
progress = gdal.TermProgress |
|
|
|
flush = sys.stdout.flush |
|
|
|
verbose = 0 |
|
quiet = 0 |
|
USE_GDALWARP = True |
|
|
|
|
|
def translate_jpgis_gml(text, dest_file, driver, create_options=None, replace_nodata_by_zero=False): |
|
"""translates JPGIS (GML) DEM file into a GDAL supported format (write access).""" |
|
|
|
create_options = create_options or [] |
|
|
|
# split document into 3 parts - header, body (tuple list) and footer |
|
header, body = text.split("<gml:tupleList>", 1) |
|
body, footer = body.rsplit("</gml:tupleList>", 1) |
|
hf = header + footer |
|
|
|
# parse xml |
|
doc = minidom.parseString(hf) |
|
lowerCorner = doc.getElementsByTagName("gml:lowerCorner")[0].childNodes[0].data.split(" ") |
|
upperCorner = doc.getElementsByTagName("gml:upperCorner")[0].childNodes[0].data.split(" ") |
|
lry = float2(lowerCorner[0]) |
|
ulx = float2(lowerCorner[1]) |
|
uly = float2(upperCorner[0]) |
|
lrx = float2(upperCorner[1]) |
|
|
|
high = doc.getElementsByTagName("gml:high")[0].childNodes[0].data.split(" ") |
|
xsize = int(high[0]) + 1 |
|
ysize = int(high[1]) + 1 |
|
|
|
startPoint = doc.getElementsByTagName("gml:startPoint")[0].childNodes[0].data.split(" ") |
|
startX = int(startPoint[0]) |
|
startY = int(startPoint[1]) |
|
|
|
psize_x = (lrx - ulx) / xsize |
|
psize_y = (lry - uly) / ysize |
|
geotransform = [ulx, psize_x, 0, uly, 0, psize_y] |
|
|
|
# create a new dataset |
|
dst_ds = driver.Create(dest_file, xsize, ysize, 1, gdal.GDT_Float32, create_options) |
|
if dst_ds is None: |
|
return "Cannot create file: " + dest_file |
|
|
|
# WKT is from: http://epsg.io/6668 |
|
dst_ds.SetProjection('GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["epsg","6668"]]') |
|
dst_ds.SetGeoTransform(geotransform) |
|
|
|
rband = dst_ds.GetRasterBand(1) |
|
if replace_nodata_by_zero: |
|
nodata_value = 0. |
|
else: |
|
nodata_value = -9999. |
|
rband.SetNoDataValue(nodata_value) |
|
|
|
# create an array initialized with nodata value |
|
narray = numpy.empty((ysize, xsize), numpy.float32) |
|
narray.fill(nodata_value) |
|
|
|
# read tuple list |
|
tuples = body.strip().split("\n") # no problem with \r character at line end: float("0.01\r") = 0.01 |
|
num_tuples = len(tuples) |
|
i = 0 |
|
sx = startX |
|
for y in range(startY, ysize): |
|
for x in range(sx, xsize): |
|
if i < num_tuples: |
|
vals = tuples[i].split(",") |
|
if len(vals) == 2 and vals[1].find("-99") == -1: |
|
narray[y][x] = float(vals[1]) |
|
i += 1 |
|
else: |
|
break |
|
if i == num_tuples: |
|
break |
|
sx = 0 |
|
|
|
# write array |
|
rband.WriteRaster(0, 0, xsize, ysize, narray.tostring()) |
|
|
|
# make sure that all data have been written |
|
dst_ds.FlushCache() |
|
|
|
if verbose: |
|
print("file: %s" % dest_file) |
|
print("name: %s" % doc.getElementsByTagName("gml:name")[0].childNodes[0].data) |
|
print("fid : %s" % doc.getElementsByTagName("fid")[0].childNodes[0].data) |
|
print("type: %s" % doc.getElementsByTagName("type")[0].childNodes[0].data) |
|
print("mesh: %s" % doc.getElementsByTagName("mesh")[0].childNodes[0].data) |
|
print("bounds : %f, %f - %f, %f" % (lry, ulx, uly, lrx)) |
|
print("cell size : %f, %f" % (psize_x, psize_y)) |
|
print("size : %d, %d" % (xsize, ysize)) |
|
print("start point : %d, %d\n" % (startX, startY)) |
|
return 0 |
|
|
|
|
|
def translate_zip(src_file, dst_file, driver, create_options=None, replace_nodata_by_zero=False): |
|
if not os.path.isfile(src_file): |
|
return "Source is not a file: " + src_file |
|
|
|
# create temporary directory |
|
temp_dir = os.path.splitext(dst_file)[0] + "_temp" + datetime.datetime.today().strftime("%Y%m%d%H%M%S") |
|
os.makedirs(temp_dir) |
|
|
|
# open zip file and translate xml files |
|
zf = zipfile.ZipFile(src_file, mode="r") |
|
namelist = zf.namelist() |
|
demlist = [] |
|
|
|
if not quiet and not verbose: |
|
progress(0.0) |
|
|
|
for i, name in enumerate(namelist): |
|
if name[-4:].lower() == ".xml" and not "meta" in name: |
|
tif_name = os.path.join(temp_dir, os.path.basename(name) + ".tif") |
|
with zf.open(name) as f: |
|
translate_jpgis_gml(f.read().decode("utf-8"), tif_name, driver, create_options, replace_nodata_by_zero) |
|
demlist.append(tif_name) |
|
if not quiet and not verbose: |
|
progress((i + 1.) / len(namelist)) |
|
zf.close() |
|
if len(demlist) == 0: |
|
return "Zip file includes no xml file: " + src_file |
|
|
|
if len(demlist) == 1: |
|
os.rename(demlist[0], dst_file) |
|
else: |
|
# create merge command |
|
gdal_merge_options = "" |
|
gdalwarp_options = "" |
|
gdal_merge_ext = "" |
|
if os.name != "nt": # nt: windows |
|
gdal_merge_ext = ".py" |
|
if quiet: |
|
gdal_merge_options += " -q" |
|
if not verbose: |
|
gdalwarp_options += " -q" |
|
if not replace_nodata_by_zero: |
|
gdal_merge_options += " -a_nodata -9999" |
|
gdalwarp_options += " -dstnodata -9999" |
|
re_non_ascii = re.compile(r"[^\x20-\x7E]") |
|
if not USE_GDALWARP and re_non_ascii.search(src_file + dst_file) is None: |
|
# write demlist into a file |
|
demlist_filename = os.path.join(temp_dir, 'demlist.txt') |
|
with open(demlist_filename, 'w') as f: |
|
f.write("\n".join(demlist)) |
|
f.write("\n") |
|
|
|
merge_command = 'gdal_merge%s%s -o "%s" --optfile "%s"' % (gdal_merge_ext, |
|
gdal_merge_options, |
|
dst_file, |
|
demlist_filename) |
|
# TODO: testing in Linux |
|
# Wildcards cannot be used for arguments now. See http://trac.osgeo.org/gdal/ticket/4542 (2012/04/08) |
|
else: |
|
merge_command = 'gdalwarp%s "%s" "%s"' % (gdalwarp_options, |
|
os.path.join(temp_dir, "*.tif"), |
|
dst_file) |
|
|
|
# do merge |
|
if not quiet: |
|
print("merging") |
|
flush() |
|
if verbose: |
|
print("execute %s" % merge_command) |
|
os.system(merge_command) |
|
|
|
# remove temporary directory |
|
shutil.rmtree(temp_dir) |
|
if not quiet: |
|
print("temporary files removed\n") |
|
flush() |
|
return 0 |
|
|
|
|
|
def unzip(src_file, dest=None): |
|
if not os.path.isfile(src_file): |
|
return False |
|
|
|
if dest is None: |
|
dest = os.path.splitext(src_file)[0] |
|
|
|
zf = zipfile.ZipFile(src_file, mode="r") |
|
zf.extractall(dest) |
|
zf.close() |
|
|
|
if verbose: |
|
print("unzipped : %s" % dest) |
|
|
|
return True |
|
|
|
def Usage(): |
|
print("=== Usage ===") |
|
print("python fgddem.py [-replace_nodata_by_zero] [-f format] [-out_dir output_directory] [-q] [-v] src_files*\n") |
|
print("src_files: The source file name(s). JPGIS(GML) DEM zip/xml files.") |
|
return 0 |
|
|
|
|
|
def main(argv=None): |
|
global verbose, quiet |
|
|
|
format = "GTiff" |
|
filenames = [] |
|
out_dir = "" |
|
replace_nodata_by_zero = False |
|
|
|
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO") |
|
os.putenv("GDAL_FILENAME_IS_UTF8", "NO") # for merging process |
|
|
|
# Parse command line arguments. |
|
i = 1 |
|
while i < len(argv): |
|
arg = argv[i] |
|
if arg == "-replace_nodata_by_zero": |
|
replace_nodata_by_zero = True |
|
elif arg == "-f": |
|
i += 1 |
|
format = argv[i] |
|
elif arg == "-out_dir": |
|
i += 1 |
|
out_dir = argv[i] |
|
elif arg == "-v": |
|
verbose = 1 |
|
elif arg == "-q": |
|
quiet = 1 |
|
elif arg == "-help" or arg == "--help": |
|
Usage() |
|
return 0 |
|
elif arg == "-version": |
|
print("fgddem.py version {}".format(script_version)) |
|
return 0 |
|
elif arg[:1] == "-": |
|
return "Unrecognised command option: %s" % arg |
|
else: |
|
# expand wildcards |
|
f = glob.glob(arg) |
|
if len(f) == 0: |
|
sys.stderr.write("File not found: %s\n" % arg) |
|
filenames += f |
|
i = i + 1 |
|
|
|
if len(filenames) == 0: |
|
return "No input files selected" |
|
|
|
# create output directory |
|
if out_dir and os.path.exists(out_dir) == False: |
|
os.makedirs(out_dir) |
|
if verbose: |
|
print("Directory has been created: %s" % out_dir) |
|
|
|
# get gdal driver |
|
driver = gdal.GetDriverByName(format) |
|
if driver is None: |
|
return "Driver %s not found" % format |
|
|
|
dst_ext = "." + (driver.GetMetadataItem(gdal.DMD_EXTENSIONS) or "unknown").split(" ")[0] |
|
err_count = 0 |
|
for i, src_file in enumerate(filenames): |
|
if not quiet: |
|
if len(filenames) > 1: |
|
print("(%d/%d): translating %s" % (i+1, len(filenames), src_file)) |
|
else: |
|
print("translating " + src_file) |
|
flush() |
|
|
|
if out_dir: |
|
filetitle, ext = os.path.splitext(os.path.basename(src_file)) |
|
dst_root = os.path.join(out_dir, filetitle) |
|
else: |
|
src_root, ext = os.path.splitext(src_file) |
|
dst_root = src_root |
|
|
|
ext = ext.lower() |
|
dst_file = dst_root + dst_ext |
|
|
|
# translate zip/xml file |
|
err = 0 |
|
if ext == ".zip": |
|
err = translate_zip(src_file, dst_file, driver, [], replace_nodata_by_zero) |
|
elif ext == ".xml" and not "meta" in src_file: |
|
with open(src_file, "r", encoding="utf-8") as f: |
|
err = translate_jpgis_gml(f.read(), dst_file, driver, [], replace_nodata_by_zero) |
|
else: |
|
err = "Not supported file: %s" % src_file |
|
|
|
if err: |
|
sys.stderr.write(err + "\n") |
|
err_count += 1 |
|
|
|
if not quiet and err_count == 0: |
|
print("completed") |
|
return 0 |
|
|
|
|
|
def float2(val, min_repeat=6): |
|
"""Increase number of decimal places of a repeating decimal. |
|
e.g. 34.111111 -> 34.1111111111111111""" |
|
repeat = 0 |
|
lc = "" |
|
for i in range(len(val)): |
|
c = val[i] |
|
if c == lc: |
|
repeat += 1 |
|
if repeat == min_repeat: |
|
return float(val[:i+1] + c * 10) |
|
else: |
|
lc = c |
|
repeat = 1 |
|
return float(val) |
|
|
|
|
|
if __name__ == "__main__": |
|
err = main(sys.argv) |
|
if err: |
|
sys.stderr.write(err + "\n") |
|
Usage() |
|
sys.exit(1) |
|
sys.exit(0) |
