概要

国土地理院からダウンロードできる数値標高モデル(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に変換することはできる。

簡単な準備

  1. python 3.x系をインストールする
  2. numpyをインストールする
  3. 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)

Next Post Previous Post