土木・建設分野に関連する仕事や研究を行われている方は、
基盤地図情報ダウンロードサービス（https://service.gsi.go.jp/kiban/app/）
から数値標高モデルをダウンロードし、地形図（DEMデータ）を作成することが多いと思います。
今回私の好きな福岡県朝倉市を例に、ダウンロードデータを処理し、geotiffデータに変換するPythonコードを紹介します
QGIS等で活用できるデータですので、ご参考になれば幸いです。
1.データダウンロード
まずは、上記のサイトからデータをダウンロードします。
今回は航空レーザ測量で得られた5mメッシュの地形データを入手します。
基盤地図情報でデータをダウンロードするためには、会員登録（ログイン）が必須ですので、登録の上、進めて下さい。
ダウンロードしたデータはすべて展開し、一つのフォルダに保存してください
（中々面倒な作業ですので、こちらもPythonで自動化できるかもしれません）
2.geotiffデータ作成
さっそくpythonで以下のコードに作業ディレクトリを追加し、実行してください。
（pythonの環境構築については、今後ご紹介します。）
実行すると、出力フォルダに画像のような結果が出力されます！
＊＊注意＊＊ 試行中のコードですので、あくまで自己責任のもと参考にしていただければと思います。
# -*- coding: utf-8 -*-
import re
import numpy as np
from osgeo import gdal
gdal.UseExceptions() # 例外を有効にする
from os.path import join,relpath
from glob import glob
def main():
#XMLを格納するフォルダ
#path = "C:\\●●\\●●\\●●" #ここにディレクトリを入れてください!
#path = "C:\\●●\\●●\\●●" #ここにディレクトリを入れてください!
path = "C:\\●●\\●●\\●●" #ここにディレクトリを入れてください!
#GeoTiffを出力するフォルダ
geopath = "C:\\●●\\●●\\●●" #ここにディレクトリを入れてください!
#ファイル名取得
files = [relpath(x,path) for x in glob(join(path,'*'))]
for fl in files:
xmlFile = join(path,fl)
#XMLを開く
with open(xmlFile, "r", encoding = "utf-8") as f:
r = re.compile("<gml:lowerCorner>(.+) (.+)</gml:lowerCorner>")
for ln in f:
m = r.search(ln)
#検索パターンとマッチした場合、スタートポジションを格納
if m != None:
lry = float2(m.group(1))
ulx = float2(m.group(2))
break
r = re.compile("<gml:upperCorner>(.+) (.+)</gml:upperCorner>")
for ln in f:
m = r.search(ln)
#検索パターンとマッチした場合、スタートポジションを格納
if m != None:
uly = float2(m.group(1))
lrx = float2(m.group(2))
break
# 検索パターンをコンパイル
r = re.compile("<gml:high>(.+) (.+)</gml:high>")
for ln in f:
m = r.search(ln)
#検索パターンとマッチした場合、縦横の領域を格納
if m != None:
xlen = int(m.group(1)) + 1
ylen = int(m.group(2)) + 1
break
startx = starty = 0
# 検索パターンをコンパイル
r = re.compile("<gml:startPoint>(.+) (.+)</gml:startPoint>")
for ln in f:
m = r.search(ln)
#検索パターンとマッチした場合、スタートポジションを格納
if m != None:
startx = int(m.group(1))
starty = int(m.group(2))
break
#numpy用にデータを格納しておく
with open(xmlFile, "r", encoding = "utf-8") as f:
src_document = f.read()
lines = src_document.split("\n")
num_lines = len(lines)
l1 = None
l2 = None
for i in range(num_lines):
if lines[i].find("<gml:tupleList>") != -1:
l1 = i + 1
break
for i in range(num_lines - 1, -1, -1):
if lines[i].find("</gml:tupleList>") != -1:
l2 = i - 1
break
#セルのサイズを算出
psize_x = (lrx - ulx) / xlen
psize_y = (lry - uly) / ylen
geotransform = [ulx, psize_x, 0, uly, 0, psize_y]
create_options = [] #空のままでOKっぽい
driver = gdal.GetDriverByName("GTiff")
#拡張子を変更する（ファイル名はそのまま）
dst = fl.replace('.xml', '.tif')
tiffFile = join(geopath,dst)
dst_ds = driver.Create(tiffFile, xlen, ylen, 1, gdal.GDT_Float32, create_options)
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"]]')
dst_ds.SetGeoTransform(geotransform)
rband = dst_ds.GetRasterBand(1)
narray = np.empty((ylen, xlen), np.float32)
narray.fill(0)
num_tuples = l2 - l1 + 1
#スタートポジションを算出
start_pos = starty*xlen + startx
i = 0
sx = startx
#標高を格納
for y in range(starty, ylen):
for x in range(sx, xlen):
if i < num_tuples:
vals = lines[i + l1].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
rband.WriteRaster(0, 0, xlen, ylen, narray.tobytes())
dst_ds.FlushCache()
def float2(str):
lc = ""
for i in range(len(str)):
c = str[i]
if c == lc:
renzoku += 1
if renzoku == 6:
return float(str[:i+1] + c * 10)
else:
lc = c
renzoku = 1
return float(str)
if __name__ == "__main__":
main()
3.結合
次に結合させます。
こちらもpythonで実行しましたが、今回は結果のみご紹介
QGISのツールを使ってもできるので、やってみてください。
結合したgeotiffをQGISに出力すると画像のような感じに！
4.まとめ
このような形で、少し気楽にgeotiffデータを作成できます。
今回はあまり詳しくご紹介できていませんが、geotiffを作成するコードの紹介がメインということでご了承ください。（不備があったら申し訳ございません。）
少しずつ、gis関連で紹介していけたらと思います。