translates DEM of National Land Numerical Information provided by MLIT into GeoTIFF.

ksjdem.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# ksjdem.py
# translates DEM of National Land Numerical Information provided by MLIT into GeoTIFF
# Copyright (c) 2012, Minoru Akagi.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

__version__ = 'beta 2012/05/30'

import sys,os
import glob
import zipfile
import codecs
from xml.dom import minidom
import numpy

try:
    from osgeo import gdal
except ImportError:
    import gdal

try:
    progress = gdal.TermProgress_nocb
except:
    progress = gdal.TermProgress

GDT2nodata_value = {gdal.GDT_Byte:255, gdal.GDT_Int32:-9999, gdal.GDT_Float32:-9999, gdal.GDT_Float64:-9999}
GDT2numpy = {gdal.GDT_Byte:numpy.int8, gdal.GDT_Int32:numpy.int32, gdal.GDT_Float32:numpy.float32, gdal.GDT_Float64:numpy.float64}
# http://gdal.org/python/osgeo.gdalconst-module.html and http://docs.scipy.org/doc/numpy/user/basics.types.html

#### data
mesh_info = {}
mesh_info['pid'] = {'G04':2}    # 1: key=id, 2: key=pid
mesh_info['level'] = {'a':[45./3600,-30./3600,80,80],'b':[4.5/3600,-3./3600,800,800],'c':[22.5/3600,-15./3600,160,160],'d':[11.25/3600,-7.5/3600,320,320]}
mesh_info['pid2level'] = {}
mesh_info['type'] = {'G04':[gdal.GDT_Float64] + [gdal.GDT_Float32] * 3 + [gdal.GDT_Byte,gdal.GDT_Float32] * 3}
mesh_info['suffix'] = {'G04':['_meshcode','_ave_elev','_max_elev','_min_elev','_min_elev_code','_max_slope','_max_slope_dir','_min_slope','_min_slope_dir','_ave_slope']}
mesh_info['color_table'] = {'G04':{0:(0,0,0),1:(255,255,255),2:(191,191,191),3:(127,127,127),4:(63,63,63),5:(0,0,0),6:(63,63,63),7:(127,127,127),8:(191,191,191)}}
####

# ksjmesh2geotiff
def ksjmesh2geotiff(src_file,dst_dir,att_numbers,driver,create_options = [],init_value = None):
    if src_file.find('tky') != -1:
        print 'unsupport file.', src_file
        sys.exit(1)

    if quiet == 0:
        print 'parsing %s...' % src_file,

    xml = minidom.parse(src_file)

    if quiet == 0:
        print 'ok'
    if debug_mode:
        print 'schemaLocation : %s' % xml.childNodes[0].attributes['xsi:schemaLocation'].value
        print 'description : %s' % xml.getElementsByTagName('gml:description')[0].childNodes[0].data

    p = os.path.basename(src_file).split('_')[0].split('-')
    p.pop()
    if len(p) == 2:
        id = '-'.join(p)
        pid, level = id.split('-')
    else:
        id = pid = p[0]
        level = mesh_info['pid2level'][pid]

    if pid in mesh_info['pid'] and level in mesh_info['level']:
        if mesh_info['pid'][pid] == 1:
            key = id
        else:
            key = pid
        psize_x,psize_y,xsize,ysize = mesh_info['level'][level]
    else:
        print 'Cannot read %s. Unsupport file.' % src_file
        sys.exit(1)

    lry, ulx = map(float2, xml.getElementsByTagName('gml:lowerCorner')[0].childNodes[0].data.split(' '))
    uly, lrx = map(float2, xml.getElementsByTagName('gml:upperCorner')[0].childNodes[0].data.split(' '))
    geotransform = [ulx, psize_x, 0, uly, 0, psize_y]

    tupleList = xml.getElementsByTagName('gml:tupleList')

    if att_numbers is None:
        att_numbers = range(2,len(mesh_info['type'][key])+1)

    if debug_mode:
        print 'bounds : %f, %f - %f, %f' % (lry,ulx,uly,lrx)
        print 'cell size : %f, %f' % (psize_x,psize_y)
        print 'number of tupleList : %d' % len(tupleList)
        print 'attribute numbers :', att_numbers

    filetitle = os.path.splitext(os.path.basename(src_file))[0]
    ext = '.tif'
    bands = 1
    att_processed = 0
    for an in att_numbers:
        an_1 = an - 1
        band_type = mesh_info['type'][key][an_1]
        if debug_mode:
            print 'an:%d suffix:%s, data type:%d' % (an,mesh_info['suffix'][key][an_1],band_type)
        if band_type:
            filename = filetitle + mesh_info['suffix'][key][an_1] + ext
            filepath = os.path.join(dst_dir, filename)

            if quiet == 0:
                print filename + '...',

            dst_ds = driver.Create(filepath, xsize, ysize, bands, band_type, create_options)
            if dst_ds is None:
                print 'Creation failed.'
                sys.exit(1)

            dst_ds.SetGeoTransform(geotransform)
            # The following Well Known Text of CRS has been quoted from http://spatialreference.org/ref/epsg/4612/ogcwkt/
            dst_ds.SetProjection('GEOGCS["JGD2000",DATUM["Japanese_Geodetic_Datum_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6612"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4612"]]')

            if init_value is None:
                nv = GDT2nodata_value[band_type]
            else:
                nv = init_value

            narray = numpy.empty((ysize,xsize),GDT2numpy[band_type])
            narray.fill(nv)

            for tl in tupleList:
                rows = tl.childNodes[0].data.strip().split("\n")
                for row in rows:
                    fields = row.split(",")
                    meshcode = fields[0]
                    if level == 'd':
                        c8 = int(meshcode[8]) - 1
                        c9 = int(meshcode[9]) - 1
                        x = int(meshcode[5]) * 40 + int(meshcode[7]) * 4 + (c8 % 2) * 2 + (c9 % 2)
                        y = ysize - 1 - (int(meshcode[4]) * 40 + int(meshcode[6]) * 4 + (c8 // 2) * 2 + (c9 // 2))
                    elif level == 'c':
                        c8 = int(meshcode[8]) - 1
                        x = int(meshcode[5]) * 20 + int(meshcode[7]) * 2 + (c8 % 2)
                        y = ysize - 1 - (int(meshcode[4]) * 20 + int(meshcode[6]) * 2 + (c8 // 2))
                    elif level == 'b':
                        x = int(meshcode[5]) * 100 + int(meshcode[7]) * 10 + int(meshcode[9])
                        y = ysize - 1 - (int(meshcode[4]) * 100 + int(meshcode[6]) * 10 + int(meshcode[8]))
                    else:   # level == 'a'
                        x = int(meshcode[5]) * 10 + int(meshcode[7])
                        y = ysize - 1 - (int(meshcode[4]) * 10 + int(meshcode[6]))

                    if fields[an_1] != 'unknown':
                        if band_type == gdal.GDT_Byte:
                            if key == 'L03-b':
                                narray[y][x] = tochiriyo2int(fields[an_1])
                            else:
                                narray[y][x] = int(fields[an_1])
                        elif band_type == gdal.GDT_Int32:
                            narray[y][x] = int(fields[an_1])
                        else:
                            narray[y][x] = float(fields[an_1])

            rband = dst_ds.GetRasterBand(1)
            rband.WriteArray(narray,0,0)
            if init_value is None:
                rband.SetNoDataValue(nv)
            if band_type == gdal.GDT_Byte:
                rband.SetColorTable(create_color_table(key))

            dst_ds = None

            if quiet == 0:
                print 'ok'

        att_processed += 1
#        if quiet == 0 and verbose == 0:
#            progress(att_processed / float(len(att_numbers)))

def tochiriyo2int(val):
    i = ord(val)
    if i < 65:
        return i - 48
    return i - 55

def create_color_table(key):
    table = gdal.ColorTable()
    entries = mesh_info['color_table'][key]
    for i in entries:
        table.SetColorEntry(i,entries[i])
    return table

def unzip_xml(src_file, dest=None):
    if not os.path.isfile(src_file):
        return None

    if dest is None:
        dest = os.path.splitext(src_file)[0]
    zf = zipfile.ZipFile(src_file, mode='r')
    names = zf.namelist()

    if verbose:
        print 'files\n%s' % '\n'.join(names)

    for n in names:
        if os.path.splitext(n)[1] == '.xml' and n.find('META') == -1:
            zf.extract(n,dest)
            zf.close()

            if verbose:
                print('unzipped : '+names[0])

            return os.path.join(dest,n)
    zf.close()

def float2(str):
    lc = ''
    for i in range(len(str)):
        c = str[i]
        if c == lc:
            renzoku += 1
            if renzoku == 5:    #
                return float(str[:i+1] + c * 11)    #
        else:
            lc = c
            renzoku = 1
    return float(str)

def Usage():
    print '=== Usage ==='
    print 'python ksjdem.py [-out_dir output_directory][-init value][-an attribute_number(s)][-q] [-v] src_files*'
    print ''
    print 'src_files: The source file name(s). JPGIS 2.1(GML) zip/xml files.'
    print 'You can use -version option to display the version of this script.'
    return 0

def main(argv=None):
    global verbose, quiet, debug_mode
    verbose = 0
    quiet = 0
    debug_mode = 0
    format = 'GTiff'
    filenames = []
    out_dir = ''
    att_numbers = [2]
    init_value = None

    gdal.AllRegister()

    if argv is None:
        argv = sys.argv
    argv = gdal.GeneralCmdLineProcessor(argv)
    if argv is None:
        return 0

    # Parse command line arguments.
    i = 1
    while i < len(argv):
        arg = argv[i]

        if arg == '-init':
            i += 1
            init_value = int(argv[i])

        elif arg == '-an':
            i += 1
            att_numbers = map(int,argv[i].split(','))

        elif arg == '-out_dir':
            i += 1
            out_dir = argv[i]

        elif arg == '-v':
            verbose = 1

        elif arg == '-q':
            quiet = 1

        elif arg == '-debug':
            debug_mode = 1

        elif arg == '-help' or arg == '--help':
            return Usage()

        elif arg == '-version':
            print 'ksjdem.py Version %s' % __version__
            return 0

        elif arg[:1] == '-':
            print 'Unrecognised command option: %s' % arg
            return Usage()

        else:
            # Expand any possible wildcards
            f = glob.glob(arg)
            if len(f) == 0:
                print 'File not found: "%s"' % arg

            filenames += f

        i = i + 1

    if verbose or debug_mode:
        sys.stdout = codecs.getwriter('shift_jis')(sys.stdout)      # Error occurs in Python >= 3.0
    if debug_mode:
        print 'args: %s' % ' '.join(argv)
        print 'files: %s' % ','.join(filenames)

    if len(filenames) == 0:
        print 'No input files selected.'
        return Usage()

    if out_dir != '' and os.path.exists(out_dir) == False:
        os.mkdir(out_dir)
        if verbose:
            print '"%s" directory has been created.' % out_dir

    driver = gdal.GetDriverByName(format)
    if driver is None:
        print 'Format driver %s not found.' % format
        return 0

    fi_processed = 0
    for src_file in filenames:
        if quiet == 0 and len(filenames) > 1:
            print 'Processing files ( %d / %d )' % (fi_processed+1, len(filenames))

        src_root, ext = os.path.splitext(src_file)
        if out_dir == '':
            dst_root = src_root
        else:
            dst_root = out_dir

        ext = ext.lower()
        if ext == '.zip':
            if quiet == 0:
                print 'extracting %s...' % src_file,

            file_in = unzip_xml(src_file,dst_root)

            if quiet == 0:
                print 'ok'

            ksjmesh2geotiff(file_in,dst_root,att_numbers,driver,[],init_value)

            os.remove(file_in)

            if verbose:
                print 'Temporary files have been removed.'
            if quiet == 0:
                print ''

        elif ext == '.xml':
            if src_file.find('META') == -1:
                if os.path.exists(dst_root) == False:
                    os.mkdir(dst_root)

                ksjmesh2geotiff(src_file,dst_root,att_numbers,driver,[],init_value)
            else:
                print '%s is META file.' % src_file

        fi_processed += 1

    if quiet == 0:
        print 'completed'
    return 0

if __name__ == '__main__':
    sys.exit(main())

ksjdem.py.droptarget.bat

@echo off
SET OSGEO4W_ROOT=C:\OSGeo4W
echo OSGEO4W home is %OSGEO4W_ROOT%

SET PATH=%PATH%;%OSGEO4W_ROOT%\bin
for %%f in (%OSGEO4W_ROOT%\etc\ini\*.bat) do call %%f

echo on
python "%~dp0ksjdem.py" "%1"

pause

ksjdem.py.md

ksjdem.py

概要

国土交通省が公開している国土数値情報のJPGIS 2.1(GML)形式のメッシュデータのうち世界測地系の標高・傾斜度(3次・4次・5次)メッシュをラスタ(GeoTIFF)に変換するPythonスクリプト。

国土数値情報ダウンロードサービス: http://nlftp.mlit.go.jp/ksj/

動作環境

WindowsではOSGeo4Wでgdal-pythonパッケージをインストールしておく。 Linuxではpython-gdalパッケージをインストールしておく。 いずれもpythonのバージョン 2.6以上が必要。

インストール・操作方法

• バッチファイルへのドラッグ&ドロップによる方法(Windowsのみ)

1. ksjdem.py.droptarget.batをテキストエディタで開き，2行目のOSGEO4W_ROOTの値をOSGeo4Wがインストールされているフォルダパスに変更する。
2. 国土数値情報ダウンロードサイトからダウンロードしたZIPファイルをksjdem.py.droptarget.batにドラッグ&ドロップすると「平均標高」がGeoTIFFで出力される。

コマンドの説明

ksjdem.py [-out_dir output_directory][-init value][-an attribute_numbers(s)] src_file(s)

• src_file(s):

JPGIS 2.1(GML)形式のXMLファイルまたはそれを含むZIPファイル。タイトルがG04で始まりタイトルにjgdを含むファイル。

• -out_dir output_directory:

出力先フォルダ。指定がないときは入力ファイルがあるフォルダにフォルダが作成される。

• -init value:

出力ラスタバンドのセル初期値。データが存在しないセルはこの値となる。指定がない場合は既定値(整数型・実数型では-9999，コードリスト型では255)で初期化され既定値が「データなし」の値として設定される。このオプションを指定した場合，ラスタバンドに「データなし」の値は設定されない。

• -an attribute_numbers(s):

  変換する属性の番号(shp属性名の後ろ3桁の番号。【出力ファイル名】を参照)を指定する。複数指定する時はコンマで区切る。既定値は2(平均標高)。

出力ファイル名

xmlファイルのファイルタイトルに次の文字列を付加したGeoTIFFファイルが出力される。

属性番号	1	2	3	4	5	6	7	8	9	10
属性名	メッシュコード	平均標高	最高標高	最低標高	最低標高コード	最大傾斜角度	最大傾斜方向	最小傾斜角度	最小傾斜方向	平均傾斜角度
接尾文字列	_meshcode	_ave_elev	_max_elev	_min_elev	_min_elev_code	_max_slope	_max_slope_dir	_min_slope	_min_slope_dir	_ave_slope

留意点

• GISソフトの中にはメッシュの縦横のサイズが異なるラスタデータを正しく読み込めないものがあります。

更新履歴

Version Beta

• 2012/05/30 公開。

Copyright (c) 2012, Minoru Akagi