今回は数値標高モデル(DEM)の解析を行ってみたいと思います。数値標高モデル(DEM)とは地形のデジタル表現のことで我らが国土地理院が公開しています。参考までにですが、これは火星のティトニウム・カズマの数値標高モデル(DEM)を三次元表示した画像です。
基盤地図情報（数値標高モデル）データのダウンロード
以下サイトでダウンロードできます。ユーザー登録が必要です。
基盤地図情報ダウンロードサービス
横須賀が含まれるデータをダウンロードします。
こんな感じで任意のフォルダにXMLのみを格納してください。承認申請のXMLは不要です。また、GeoTiff出力用のフォルダも作成してください。
基盤地図情報（数値標高モデル）データを解析する
ファイル名について
基盤地図情報（数値標高モデル）データには 5mメッシュと10mメッシュがありますが、今回はより細かい 5mメッシュを使用します。
ファイル名は以下の形式になっています。
FG-GML-aaaa-bb-cc-DEM5X-yyyymmdd.xml
・フォーマット
aaaaは1次メッシュ番号
bbは2次メッシュ番号
ccは3次メッシュ番号
Xは測量方法 - A:レーザー測量、B:写真測量
yyyymmddは年月日
「FG-GML-5439-13-24-DEM5A-20130702.xml」みたいな感じですね。
ファイルの中身について
重要な情報のみを記述します。
・緯度経度
これをもとにGeoTiffに変換します。
<gml:lowerCorner>35.175 139.6</gml:lowerCorner> <gml:upperCorner>35.183333333 139.6125</gml:upperCorner>
・GridEnvelope
メッシュを構成するグリッドセルの配列数を記載します。low属性は北西端を示す番号です。常に(0,0)です。
high属性は南東端を示す番号です。5mメッシュは(224,149)の固定です。つまり、5mメッシュのXMLは225×150＝33750行の標高データが格納されているということです。ただし、後述しますが、startPointの値によってはその限りではありません。
<gml:limits> <gml:GridEnvelope> <gml:low>0 0</gml:low> <gml:high>224 149</gml:high> </gml:GridEnvelope> </gml:limits>
・DataBlock
各グリッドセルの標高データが記載されます。tupleListタグで囲まれたデータが標高データになります。水面や標高データがない場合は、標高の値が"-9999."になります。
<gml:DataBlock> <gml:rangeParameters> <gml:QuantityList uom="DEM構成点"/> </gml:rangeParameters> <gml:tupleList> 地表面,23.16 地表面,23.16 地表面,23.61 ・ ・ ・ </gml:tupleList> </gml:DataBlock>
・coverageFunction
sequenceRuleでデータの並びを指定します。orderの"+x"はxが正の方向(西→東)であることを、"-y"はyは負の方向(北→南)であることを示します。北西端から南東端に向けて次のような並びで一列にしたデータがtupleListタグ内に配置されます。
gml:coverageFunction> <gml:GridFunction> <gml:sequenceRule order="+x-y">Linear</gml:sequenceRule> <gml:startPoint>0 0</gml:startPoint> </gml:GridFunction> </gml:coverageFunction>
startPointでグリッドセルの開始点を指定します。先頭部分にデータが存在しない場合、このようになります。
<gml:startPoint>5 1</gml:startPoint>
これの意味はデータのスタートがX方向が05列、Y方向が01行ということです。また、末尾部分にデータが存在しないグリッドは省略し、省略した旨を記載するタグは存在しません。
ですので、startPointを考慮してもデータ数が合わない場合は末尾部分のデータは省略されているということです。
NumpyとGDALのインストール
配列を扱う場合はNumpyですね。また、GDALも必要なので、インストールをお願いします。
サンプルソース
ちょっと無理矢理感は否めませんが、参考までに。
# -*- coding: utf-8 -*- import re import numpy as np import gdal from os.path import join,relpath from glob import glob def main(): #XMLを格納するフォルダ path = "D:\python\DEM" #GeoTiffを出力するフォルダ geopath = "D:\python\GeoTiff" #ファイル名取得 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.tostring()) 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()
結果を確認してみます。
Arcpyを使用してJ-SHISの表層地盤データを分析してみよう（横須賀編） ～ArcpyでShape読込、Shape出力、属性検索、空間検索、フィールド演算～ - GIS奮闘記 で使用した日本地図とGeoTiffを重ねあわせてみます。
おぉー。なんだかそれっぽいのが表示されました。しかし、まだ以下をする必要がありますね。
・取り込んだGeoTiffを一枚のラスターに変換（色の設定を一ファイルずつ行うのは非現実的なので。ちなみに、全部で239個のGeoTiffを取り込みました）。
・色の設定（白黒だと若干見づらいので）
ただ、ここからはArcpyを使用しますので、ArcGISのライセンスが無いという方はすみません。以下スクリプトでGeoTiffを一枚のラスタにします。
# -*- coding: utf-8 -*- import arcpy from glob import glob from os.path import join,relpath #GeoTiffを格納しているディレクトリ path = "D:\python\GeoTiff" #ファイル名取得 files = [relpath(x,path) for x in glob(join(path,'*'))] #MosaicToNewRaster_managementの引数を作成 for fl in files: raster = raster + fl raster = raster + ";" raster = raster[:-1] #以下のように任意のフォルダをワークスペースにする arcpy.env.workspace = "D:\python\Raster" #ラスタを一枚にする arcpy.MosaicToNewRaster_management(input_rasters = raster, output_location = "D:\python\Raster", \ raster_dataset_name_with_extension = "DEM.tif",\ pixel_type = "8_BIT_UNSIGNED", number_of_bands = "1", mosaic_method = "LAST",mosaic_colormap_mode = "FIRST")
こんな感じで一枚のラスタにすることができました。
ラスタのプロパティで色の設定をします。
それっぽい感じになりましたね。
あとは横須賀の範囲でラスタをクリップしたいのですが、Spatial Analystのライセンスを持っていないのでそこはあきらめます。単純なクリップ機能を使用すればある程度は抽出できますね。まずは以下スクリプトで横須賀のフィーチャを抽出します。
mxd = arcpy.mapping.MapDocument("CURRENT") df = arcpy.mapping.ListDataFrames(mxd)[0] #読み込んでいる日本地図を利用 addLayer = arcpy.mapping.Layer("japan_ver80") #レイヤを追加 arcpy.mapping.AddLayer(df, addLayer, "BOTTOM") #属性検索 zokuseiSearch = arcpy.SelectLayerByAttribute_management(addLayer.name,"NEW_SELECTION","SIKUCHOSON = '横須賀市'") #Shape出力 arcpy.FeatureClassToFeatureClass_conversion(zokuseiSearch, "D:\python\Raster", "yokosuka")
横須賀のみのレイヤの作成に成功しました。
次にラスタのクリップを行います。ラスタデータセット作成用のGDBを一つ作成しておく必要があります。
こんな感じでラスタを横須賀の範囲くらいまでクリップすることができました（座標の設定などもあったのでここはArcMapを使ってしまいました）。
色の設定をします。標高が低い順に水色→緑→茶色→白 と色が設定されました。
Arcpyを使用してJ-SHISの表層地盤データを分析してみよう（横須賀編） ～ArcpyでShape読込、Shape出力、属性検索、空間検索、フィールド演算～ - GIS奮闘記 で横須賀は丘陵が多いことがわかりましたが、そのせいか横須賀は全体的に標高が高いですね（たしかに山を切り開いてつくった町も多いですね）。
あとは街区データなどがあれば街ごとに標高の高低差を確認することができ、もっと面白くなりそうですね！ZMAP-TOWN2を買おうかしら。
以上、本日は「Pythonで基盤地図情報の数値標高モデルを解析してGeoTiffに変換する（横須賀編） ～Numpy、GDAL、Arcpyを駆使して横須賀の標高（地盤高）を視覚化する～」でした。
