概要
国土地理院からダウンロードできる数値標高モデル(Digital Elevation Model:DEM)はJPGIS(GML)形式で書かれている。JPGIS(GML)形式はXML文書なのでPythonでこれを読み取ってGISソフトウェアなどで取り扱いやすいGeotiff形式に変換する。
方法
結論から言えば、このページにあるPythonコードを適当にjpgis_convert.pyなどの名前で保存して
python jpgis_convert.py FG-GML-XXXX-XX-DEM5A.zip FG-GML-XXXX-XX-DEM5B.zip FG-GML-XXXX-XX-DEM5A.zipのような感じで実行すればFG-GML-XXXX-XX-DEM5A.tifファイルができる。コマンドラインツールとして基本的な機能が何も整っていないが最低限JPGIS(GML)をGeotiffに変換することはできる。
簡単な準備
- python 3.x系をインストールする
- numpyをインストールする
- GDALをインストールする
以上
Pythonコード
例によって自由に使っていただいて構わないがいつでも正しく動作する保証はしません。気が向いたら改良を加えるかもしれません。メタデータが欲しければ適当にprintしてください。
#!python3
# -*- coding:utf-8 -*-
import os
import sys
import codecs
import zipfile
import csv
import numpy as np
import xml.etree.ElementTree as et
import osgeo.gdal,osgeo.osr
def convert(input_filename, output_filename):
ns = {
'ns' : 'http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema',
'gml' : 'http://www.opengis.net/gml/3.2',
'xsi' : 'http://www.w3.org/2001/XMLSchema-instance',
'xlink' : 'http://www.w3.org/1999/xlink' }
nodata = -1
dataset = []
mllon, mllat, mulon, mulat = 180, 90, -180, -90 # most upper, lower
with zipfile.ZipFile(input_filename, 'r') as zf:
filelist = zf.namelist()
for filename in filelist:
print('loading '+filename)
with zf.open(filename, 'r') as file:
text = file.read().decode('utf_8')
root = et.fromstring(text)
dem = root.find('ns:DEM', ns)
coverage = dem.find('ns:coverage', ns)
envelope = coverage.find('gml:boundedBy//gml:Envelope', ns)
lower = envelope.find('gml:lowerCorner', ns).text
upper = envelope.find('gml:upperCorner', ns).text
grid = coverage.find('gml:gridDomain//gml:Grid//gml:limits//gml:GridEnvelope', ns)
low = grid.find('gml:low', ns).text
high = grid.find('gml:high', ns).text
tuplelist = coverage.find('gml:rangeSet//gml:DataBlock//gml:tupleList', ns).text
gridfunc = coverage.find('gml:coverageFunction//gml:GridFunction', ns)
rule = gridfunc.find('gml:sequenceRule', ns)
start = gridfunc.find('gml:startPoint', ns).text
if rule.attrib['order'] != '+x-y':
print('warning sequence order not +x-y')
if rule.text != 'Linear':
print('warning sequence rule not Linear')
s = np.array(lower.split(' '), dtype=np.float64)
llat, llon = s[0], s[1] #lower latitude, lower longitude
if llat < mllat: mllat = llat
if llon < mllon: mllon = llon
s = np.array(upper.split(' '), dtype=np.float64)
ulat, ulon = s[0], s[1] #upper latutude, upper longitude
if ulat > mulat: mulat = ulat
if ulon > mulon: mulon = ulon
s = low.split(' ')
lowx, lowy = int(s[0]), int(s[1])
s = high.split(' ')
highx, highy = int(s[0]), int(s[1])
datapoints = tuplelist.splitlines()
raster = np.zeros([highy - lowy + 1, highx - lowx + 1])
raster.fill(nodata)
s = start.split(' ')
x, y = int(s[0]), int(s[1])
for datapoint in datapoints:
s = datapoint.split(',')
if len(s) == 1: continue
desc,value = s[0],np.float64(s[1])
if desc == '地表面': raster[y][x] = value
else: raster[y][x] = nodata
x += 1
if x > highx:
x = 0
y += 1
if y > highy: break
data = {}
data['llat'] = llat
data['llon'] = llon
data['ulat'] = ulat
data['ulon'] = ulon
data['xsize'] = highx - lowx + 1
data['ysize'] = highy - lowy + 1
data['raster'] = raster
dataset.append(data)
file.close()
zf.close()
xsize = dataset[0]['xsize']
ysize = dataset[0]['ysize']
for data in dataset:
if data['xsize'] != xsize or data['ysize'] != ysize:
print('warning grid size not equal for dataset')
merged_raster = np.zeros([ysize*10, xsize*10])
merged_raster.fill(nodata)
xbandsize = mulon - mllon
ybandsize = mulat - mllat
for data in dataset:
xbase = int(round(10*(data['llon'] - mllon)/xbandsize))
ybase = int(round(10*(mulat - data['ulat'])/ybandsize))
raster = data['raster']
for i in range(ysize):
for j in range(xsize):
merged_raster[ybase*ysize + i][xbase*xsize + j] = raster[i][j]
pixelx = (mulon - mllon)/(xsize*10 - 1)
pixely = (mulat - mllat)/(ysize*10 - 1)
trans = [mllon, pixelx, 0, mulat, 0, -pixely]
srs = osgeo.osr.SpatialReference()
srs.ImportFromEPSG(4612)
driver = osgeo.gdal.GetDriverByName('GTiff')
output = driver.Create(output_filename, xsize*10, ysize*10, 1, osgeo.gdal.GDT_Float32)
output.GetRasterBand(1).WriteArray(merged_raster)
output.GetRasterBand(1).SetNoDataValue(nodata)
output.SetGeoTransform(trans)
output.SetProjection(srs.ExportToWkt())
output.FlushCache()
def main(argv):
for i in range(1,len(argv)):
convert(argv[i], os.path.splitext(os.path.basename(argv[i]))[0] + '.tif')
if __name__ == '__main__':
main(sys.argv)