PythonでDEMのXMLファイルをgeotiffファイルに変換する

PythonでDEMのXMLファイルをgeotiffファイルに変換します．

ライブラリのインストール
geotiffファイルの変換
geotiffファイルの変換実行
まとめ
参考サイト

# ライブラリのインストール

下記のコードでライブラリをインストールします．

!sudo apt-get install python-gdal python3-gdal
!sudo apt-get update
!sudo apt-get install -y build-essential
!sudo apt-get install -y gdal-bin
!sudo apt-get install -y libgdal-dev

!gdal-config --version
!pip install setuptools==57.4.0
!pip install GDAL==$(gdal-config --version) --global-option=build_ext --global-option="-I/usr/include/gdal"

# geotiffファイルの変換

下記のコードでgeotiffファイルへ変換します.

%%writefile zdem2tif.py
import os
import sys
import zipfile
import numpy as np
import xml.etree.ElementTree as et
import osgeo.gdal
import osgeo.osr


class XMLContents(object):
    def __init__(self):
        self.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'
                   }

    def read_xml(self, zf, filename, dem_mode):
        if dem_mode == 'DEM10':
            self.posix = 0
            self.posiy = 0
        elif dem_mode == 'DEM5':
            _split = os.path.basename(filename).split('-')
            self.posix = int(_split[4][1])
            self.posiy = int(_split[4][0])

        with zf.open(filename, 'r') as file:
            self.text = file.read().decode('utf_8')
            self.root = et.fromstring(self.text)
            self.dem = self.root.find('ns:DEM', self.ns)
            self.coverage = self.dem.find('ns:coverage', self.ns)
            self.envelope = self.coverage.find(
                'gml:boundedBy//gml:Envelope', self.ns)
            self.lower = self.envelope.find('gml:lowerCorner', self.ns).text
            self.upper = self.envelope.find('gml:upperCorner', self.ns).text
            self.grid = self.coverage.find(
                'gml:gridDomain//gml:Grid//gml:limits//gml:GridEnvelope', self.ns)
            self.low = self.grid.find('gml:low', self.ns).text
            self.high = self.grid.find('gml:high', self.ns).text
            self.tuplelist = self.coverage.find(
                'gml:rangeSet//gml:DataBlock//gml:tupleList', self.ns).text
            self.gridfunc = self.coverage.find(
                'gml:coverageFunction//gml:GridFunction', self.ns)
            self.rule = self.gridfunc.find('gml:sequenceRule', self.ns)
            self.start = self.gridfunc.find('gml:startPoint', self.ns).text

            if self.rule.attrib['order'] != '+x-y':
                print('warning sequence order not +x-y')
            if self.rule.text != 'Linear':
                print('warning sequence rule not Linear')

            file.close()


class DEMInArray(object):
    def __init__(self, contents):
        self._noValue = -1.
        self.llat, self.llon = np.array(
            contents.lower.split(' '), dtype=np.float64)
        self.ulat, self.ulon = np.array(
            contents.upper.split(' '), dtype=np.float64)
        self.lowx, self.lowy = np.array(
            contents.low.split(' '), dtype=np.int)
        self.highx, self.highy = np.array(
            contents.high.split(' '), dtype=np.int)
        self.sizex, self.sizey = self.get_size()
        self.x_init, self.y_init = np.array(
            contents.start.split(' '), dtype=np.int)
        self.datapoints = contents.tuplelist.splitlines()
        self.raster = self.get_raster()

    def get_size(self):
        sizex = self.highx - self.lowx + 1
        sizey = self.highy - self.lowy + 1
        return sizex, sizey

    def get_raster(self):
        _x, _y = self.x_init, self.y_init
        raster = np.zeros([self.sizey, self.sizex])
        raster.fill(self._noValue)
        for dp in self.datapoints:
            if _y > self.highy:
                print('DEM datapoints are unexpectedly too long')
                break
            s = dp.split(',')
            if len(s) == 1:
                continue
            desc, value = s[0], np.float64(s[1])
            # if desc == '地表面':
            raster[_y, _x] = value if value > -9998 else self._noValue
            _x += 1
            if _x > self.highx:
                _x = 0
                _y += 1
        return raster


class DEMOutArray(object):
    def __init__(self):
        self._fillValue = -9999.

    def initialize_raster(self):
        raster = np.zeros([self.sizey, self.sizex])
        raster.fill(self._fillValue)
        return raster

    def get_trans(self):
        # +x-y
        return [self.llon, self.resx, 0, self.ulat, 0, -self.resy]

    def update_raster(self, tile, posiy, posix):
        xmin = posix * tile.shape[1]
        ymin = posiy * tile.shape[0]
        xmax = xmin + tile.shape[1]
        ymax = ymin + tile.shape[0]
        self.raster[ymin:ymax, xmin:xmax] = tile[::-1]

    def get_size(self):
        sizex = (self.highx - self.lowx + 1) * self.tilex
        sizey = (self.highy - self.lowy + 1) * self.tiley
        return sizex, sizey

    def mesh_info_from_zipname(self, zipname):
        '''
        Ref: Table 4 of https://www.stat.go.jp/data/mesh/pdf/gaiyo1.pdf
        '''
        _split = os.path.basename(zipname).split('-')
        _lat1 = int(_split[2][:2])
        _lon1 = int(_split[2][2:])
        _lat2 = int(_split[3][0])
        _lon2 = int(_split[3][1])
        self.llat = _lat1 * 0.66666666666 + _lat2 * 0.08333333333
        self.llon = _lon1 + 100. + _lon2 * 0.125
        self.ulat = self.llat + 0.08333333333
        self.ulon = self.llon + 0.125

        if _split[4][:5] == 'DEM10':
            self.mode = 'DEM10'
            self.tilex, self.tiley = 1, 1
            self.lowx, self.lowy = 0, 0
            self.highx, self.highy = 1124, 749
            self.sizex, self.sizey = self.get_size()
            self.resx, self.resy = 1.1111111111e-04, 1.1111111111e-04
        elif _split[4][:4] == 'DEM5':
            self.mode = 'DEM5'
            self.tilex, self.tiley = 10, 10
            self.lowx, self.lowy = 0, 0
            self.highx, self.highy = 224, 149
            self.sizex, self.sizey = self.get_size()
            self.resx, self.resy = 5.5555555555e-05, 5.5555555555e-05
        else:
            raise ValueError('Unexpected definition of DEM:', _split[4])

        self.raster = self.initialize_raster()

    def save_raster(self, out_filename):
        trans = self.get_trans()
        srs = osgeo.osr.SpatialReference()
        srs.ImportFromEPSG(4612)
        driver = osgeo.gdal.GetDriverByName('GTiff')
        output = driver.Create(out_filename, self.sizex, self.sizey,
                               1, osgeo.gdal.GDT_Float32)
#         output.GetRasterBand(1).WriteArray(self.raster) <- 修正@2020.06.10
        output.GetRasterBand(1).WriteArray(self.raster[::-1])
        output.GetRasterBand(1).SetNoDataValue(self._fillValue)
        output.SetGeoTransform(trans)
        output.SetProjection(srs.ExportToWkt())
        output.FlushCache()


def convert(in_filename, out_filename):
    output_array = DEMOutArray()
    output_array.mesh_info_from_zipname(in_filename)
    dem_mode = output_array.mode

    with zipfile.ZipFile(in_filename, 'r') as zf:
        filelist = zf.namelist()
        for filename in filelist:
            print('loading', filename)
            contents = XMLContents()
            contents.read_xml(zf, filename, dem_mode)
            posix, posiy = contents.posix, contents.posiy
            dem_tile = DEMInArray(contents)
            output_array.update_raster(dem_tile.raster, posiy, posix)
        zf.close()
    print('save', out_filename)
    output_array.save_raster(out_filename)


def main(argv):
    for i in range(1, len(argv)):
        convert(argv[i], argv[i].replace('.zip', '.tif'))


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

# geotiffファイルの変換実行

下記のコードでgeotiffファイルへ変換を実行します.

!python zdem2tif.py ../kure_city_data/geograph/dem/5Mmesh_zip/FG-GML-5132-03-DEM5A.zip \

                    ../kure_city_data/geograph/dem/5Mmesh_zip/FG-GML-5132-36-DEM5B.zip