サンプルソース
ちょっと無理矢理感は否めませんが、参考までに。
import re
import numpy as np
import gdal
from os.path import join,relpath
from glob import glob
def main():
path = "D:\python\DEM"
geopath = "D:\python\GeoTiff"
files = [relpath(x,path) for x in glob(join(path,'*'))]
for fl in files:
xmlFile = join(path,fl)
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
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 = []
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を一枚のラスタにします。
import arcpy
from glob import glob
from os.path import join,relpath
path = "D:\python\GeoTiff"
files = [relpath(x,path) for x in glob(join(path,'*'))]
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 = '横須賀市'")
arcpy.FeatureClassToFeatureClass_conversion(zokuseiSearch,
"D:\python\Raster",
"yokosuka")
横須賀のみのレイヤの作成に成功しました。
次にラスタのクリップを行います。ラスタデータセット作成用のGDBを一つ作成しておく必要があります。
こんな感じでラスタを横須賀の範囲くらいまでクリップすることができました（座標の設定などもあったのでここはArcMapを使ってしまいました）。
色の設定をします。標高が低い順に水色→緑→茶色→白 と色が設定されました。
Arcpyを使用してJ-SHISの表層地盤データを分析してみよう（横須賀編） ～ArcpyでShape読込、Shape出力、属性検索、空間検索、フィールド演算～ - GIS奮闘記 で横須賀は丘陵が多いことがわかりましたが、そのせいか横須賀は全体的に標高が高いですね（たしかに山を切り開いてつくった町も多いですね）。
あとは街区データなどがあれば街ごとに標高の高低差を確認することができ、もっと面白くなりそうですね！ZMAP-TOWN2を買おうかしら。
以上、本日は「Pythonで基盤地図情報の数値標高モデルを解析してGeoTiffに変換する（横須賀編） ～Numpy、GDAL、Arcpyを駆使して横須賀の標高（地盤高）を視覚化する～」でした。
