地形の高さ (DEM, DSM, nDSM) を点群データ(オープンナガサキ)から生成する
高さ情報とは？
建物の高さを知りたいな〜と思っていても、、、
測量や衛星データを扱っている人はご存知の人も多いですが、皆さんの想像上以上に考慮することが多いです。
まず大切なことはどこからの高さなのかということです。
この常識は分野や前提ごとに変わるので意識しなければなりません。
以下より引用
ざっくりと以下まで考慮できれば十分だと思います。
- 地球楕円体
- ジオイド
- 標高
- 地表面測量
全て考えていたら正確かもしてませんがそれだけやることは増えます。
そろぞれの目的にあった高さ情報を利用しましょう。
DEM: Digital Terrain Model DSM: Digital Surface Model
-
DEM (Digital Terrain Model)とは、
海面からの地面の高さで、標高のことです。 - DSM (Digital Surface Model)とは、地表面の物体を含む高さのことです。
以下より引用
nDSM: Normalized Digital Surface Model
今回は nDSM と言って建物や樹木などの直感的な地表面の地物の高さを示します。
英語では、 Normalized Digital Surface Model のことです。
こちらの図が理解しやすいと思います。
以下より引用
点群データ
利用する点群データの紹介です。
オープンナガサキ
長崎県が点群データを C.C ライセンスで無償公開してくれています。
今回はこの点群データを利用して DSM を計算していきます。
この点群は RIDAR というレーザー測量によって計測されています。
場所は長崎県長崎市の長崎駅です。
The ナガサキ！
可視化
可視化には CloudCompare を利用しました。
あらゆる角度や距離から可視化が可能です。
3D 可視化は楽しいですね！
位置情報
オープンナガサキのデータには位置情報も付与されています。
QGIS を使って可視化します。
しかしながら、オープンナガサキの位置情報は座標情報が長崎になっていないので修正します。
こちらは点群の地理空間処理の PDAL を使用します。
Point データの GDAL ですね！
PATH_LAS = '****_org.las'
PATH_LAS_T = PATH_LAS.replace('.las', '_trandlate.las')
# 長崎は EPSG:6669
!pdal translate -i {PATH_LAS} -o {PATH_LAS_T} --writers.las.a_srs="EPSG:6669"
では可視化します。
長崎駅も場所も正しく細部まで見えています。
Z軸の取得
この点群データで XY軸 は地理空間の座標で、 Z軸 は地表面の高さ情報になっています。
なので Z軸 から DSM を取得します。
cfg.json を以下で保存します。
{
"pipeline":[
{
"type":"readers.las",
"filename":"****_org_trandlate.las"
},
{
"type":"writers.gdal",
"filename":"****_org_trandlate_dsm.tif",
"dimension":"Z",
"data_type":"float32",
"output_type":"mean",
"resolution": 1.0
}
]
}
この設定情報を利用して PDAL のパイプライン処理を行います。
pdal pipeline ./cfg.json -v 1
Python で可視化します。
import matplotlib.pyplot as plt
import tifffile
PATH_DSM = '****_org_trandlate_dsm.tif'
img = tifffile.imread(PATH_DSM)
img[img < 0.5] = 0.5
plt.figure(figsize=(20,12), facecolor='w')
plt.imshow(img)
plt.colorbar()
plt.savefig('dsm.png')
plt.show();
建物が 20 ~ 30m でそれっぽい値が取得できていますね！
これが DSM になります。
地物の高さ
nDSM を計算しましょう。
DEM の取得
さて、ここで nDSM を取得するために DEM を求めます。
今回は解像度が高く、 C.C 4.0 ライセンスの 国土地理院 の DEM を使用します。
場所は運悪く、２面に重複しているの２つ取得します。
DEM を可視化します。
DEM( 青色 ) と点群データを可視化します。
次に、 DEM を切り抜きます。
点群の位置座標と画像のサイズを QGIS で見てみます。
gdalwarp -te 34000.0000000000000000 -27750.0000000000000000 34999.9899999999979627 -27000.0000000000000000 \
-te_srs EPSG:6669 -t_srs EPSG:6669 -ts 1001 751 \
****_org_trandlate_dsm.tif dem_crop.tif
切り抜いた DEM( 赤色 ) を可視化します。
nDSM の計算
を実装します。
dem = tifffile.imread(PATH_DEM)
dsm = tifffile.imread(PATH_DSM)
n_dem = dsm - dem
n_dem[n_dem < 0] = 0
tifffile.imwrite(PATH_nDSM, n_dem)
では、nDSM( 緑色 ) を可視化します。
地表面の地物だけの高さを取得できました！
nDSM による 3D 可視化
以下の以前の記事の方法で 3D 可視化をして確認してみます。
Google Map
Google Map での可視化を行います。建物の光学情報が確認できます。
動画だとよく分かりますね。県庁？かながかなり高いですね。
DEM( 赤色 ) も同時に可視化してみると、補正された部分がわかります。
衛星データ
Sentinel-2 という光学衛星で可視化します。
aws s3 cp s3://sentinel-s2-l2a/tiles/52/S/EB/2023/8/27/0/ . --recursive --request-payer requester
上記で AWS の S3 から取得可能です。
長崎駅あたりはこんな感じです。
バンド8 (近赤外)を 3D 可視化します。
人工物などの金属が白色で明るくなっていますね。
衛星データと点群データでこのような解析ができます。
まとめ
nDSM の求め方とその利用例、他のデータとの解析・可視化を行いました。
点群データや衛星データの 近未来のデータ を活用して新しい解析をしていきたいですね！
これらを利用することで、建物の数がわかると都市の発展具合、開発状況がよくわかります。
他にも地形変化や自然災害、カーナビへの利用などワクワクしますね！
最後まで読んでいただきありがとうございます！*
皆さんのお役に立てられると光栄です。
おまけ
こちら以外にも記事執筆をしているのでご参考になれば幸いです
衛星データ解析として、宙畑のライターもしています。
SAR 解析をよくやっていますが、地理空間や衛星データ、点群データに関心があります。
勉強している人は好きなので楽しく絡んでくれると嬉しいです。
Discussion