兵庫県1m DEMをさわってみる

はじめに

https://web.pref.hyogo.lg.jp/press/20200110_4566.html によると、兵庫県内の1mメッシュ(!)のDSM, DEMm CS立体図がダウンロードできます。

これは、これまでに類を見ない、大変なことです。全国整備されている基盤地図情報は5mメッシュと10mメッシュです。

さらに、ライセンスはCC-BYです。

そこで、試しに DEM を落として、様子を見てみました。

データのダウンロード

https://www.geospatial.jp/ckan/dataset/2010-2018-hyogo-geo-potal から CS立体図, DSM, DEM がダウンロードできます。

DEMは https://www.geospatial.jp/ckan/dataset/2010-2018-hyogo-geo-dem にて取得できます。

ZIPファイルがおかれていますので、必要なものをダウンロードします。

求める地域とZIPファイル名との関係が不明な場合には、"Indexmap (GeoJSONファイル)"を使用します。

"Indexmap (GeoJSONファイル)"をクリックすると、地図が表示されます。さらに地図上部のタブで「GeoJSON」を選択します。「Map Viewer」ではGeoJSONの属性が表示されません。メッシュをクリックすると"meshID"という属性の値が表示されます。これがZIPファイル名のプリフィックスです。

今回は、旧日高町付近を狙って 05LF90 にしました。

GDALに食わせる

ZIPを展開すると、4つのtxtファイル(DEM_05LF901_1g.txt, DEM_05LF902_1g.txt, DEM_05LF903_1g.txt, DEM_05LF904_1g.txt)が現れました。

まずは様子を見る

うち一つの様子をgdalinfoで見てみます。

gdalinfo DEM_05LF901_1g.txt  

多少時間がかかります。たぶん全部読まないと画像サイズ等が計算できないからだろうと思います。

Driver: XYZ/ASCII Gridded XYZ  
Files: DEM_05LF901_1g.txt  
Size is 2000, 1500  
Origin = (40000.000000000000000,-58500.000000000000000)  
Pixel Size = (1.000000000000000,1.000000000000000)  
Corner Coordinates:  
Upper Left  (   40000.000,  -58500.000)  
Lower Left  (   40000.000,  -57000.000)  
Upper Right (   42000.000,  -58500.000)  
Lower Right (   42000.000,  -57000.000)  
Center      (   41000.000,  -57750.000)  
Band 1 Block=2000x1 Type=Float32, ColorInterp=Undefined  
  Min=5.090 Max=98.760  

ということで、XYZファイルで、GDALで読める(OGRではないです)ことが分かりました。

空間参照系を推定する

このデータで使われているのは、たぶん、平面直角座標系のV系です。

まず、単位が度でないのは明らかです。そうすると、わが国の陸域ではUTMか平面直角座標系のいずれかです。

Y値が負になっているので、UTMは外れます。UTMの原点が赤道上であるため、北半球ではY値が負になることはありえません。なお、普通は南半球でも正数となるようにオフセットを足します(北半球とは別の系となります)。

平面直角系の何系かは適用範囲から分かります。わかりやすい平面直角座標系から、兵庫県ならV系であることが分かります。

平面直角座標系V系の空間参照系コードはEPSG:6673です。

GeoTIFFにする

gdalwarpで複数のXYZファイルをまとめて一つのGeoTIFFにしてみます。

gdalwarp -s_srs EPSG:6673 DEM_05LF901_1g.txt DEM_05LF902_1g.txt DEM_05LF903_1g.txt DEM_05LF904_1g.txt DEM_05LF90.tif  

やはり、XYZファイルの読み込みに時間がかかるので、ちょっと待ってください。

今回は一つのファイルにまとめましたが、必ずしもまとめる必要はありませんが、XYZファイルだと時間がかかるので、GeoTIFF等の一般的なラスタへの変換はやっておいた方がいいかと思います。

gdalinfo DEM_05LF90.tif  
Driver: GTiff/GeoTIFF  
Files: DEM_05LF90.tif  
Size is 4000, 3000  
Coordinate System is:  
PROJCRS["JGD2011 / Japan Plane Rectangular CS V",  
(中略)  
Data axis to CRS axis mapping: 2,1  
Origin = (40000.000000000000000,-57000.000000000000000)  
Pixel Size = (1.000000000000000,-1.000000000000000)  
Metadata:  
  AREA_OR_POINT=Area  
Image Structure Metadata:  
  INTERLEAVE=BAND  
Corner Coordinates:  
Upper Left  (   40000.000,  -57000.000) (134d46'27.00"E, 35d29' 7.51"N)  
Lower Left  (   40000.000,  -60000.000) (134d46'26.47"E, 35d27'30.16"N)  
Upper Right (   44000.000,  -57000.000) (134d49' 5.70"E, 35d29' 6.90"N)  
Lower Right (   44000.000,  -60000.000) (134d49' 5.12"E, 35d27'29.55"N)  
Center      (   42000.000,  -58500.000) (134d47'46.07"E, 35d28'18.53"N)  
Band 1 Block=4000x1 Type=Float32, ColorInterp=Gray  

ピクセルサイズも変わらないですし、大丈夫ですね。

QGISで眺める

まずは、ラスタレイヤを追加して、グレースケールで見てみます。


出典:兵庫県_全域DEM(2010年度~2018年度)、兵庫県

次に陰影図を見ます。陰影図の作成は「ラスタ」→「解析」→「陰影図 (hillshade) ...」を選択します。


出典:兵庫県_全域DEM(2010年度~2018年度)、兵庫県

田んぼのあぜが見えます。

3D表示もやってみましょう。詳細は標準装備のQGISで三次元表示ができるよ!を参照して下さい。


出典:兵庫県_全域DEM(2010年度~2018年度)、兵庫県

今回は、標高のグレースケールだけで表示していますが、空中写真等をかぶせると、結構見ごたえのある絵になります。

おわりに

いかがだったでしょうか。

兵庫県さんが出されたのは1mメッシュという非常に細かいデータです。これを自由に扱えるのですから、とんでもない時代になったと思います。

このデータはGDALが難なく食べられるので、特に追加的なコード作成もなく使いまわせます。ただし、XYZファイルは読み込みに時間がかかるので、GeoTIFF等に変換しておいた方が良いと思います。

現時点の私は、兵庫県を舞台にして何かしているわけではないので、実際に使う場面は無いのですが、そんな者にも触らせてくれるのはありがたいです。

本記事のライセンス


この記事は クリエイティブ・コモンズ 表示 4.0 国際 ライセンスの下に提供されています。