兵庫県土の1m精度3次元データをWeb地図上に3D可視化【DSM】
1. 概要
2020年1月、兵庫県が全国で初めて「全県土分の高精度３次元データ」をオープンデータとして公開しました。この中には地表面データ（DEM: Digital Elevation Model）だけでなく建物・樹木などの高さを含んだ地球表面データ（DSM: Digital Surface Model）も含まれており、兵庫県の実際の3次元形状が1メートル四方の精度で把握できることになります。
今回は、DSMデータ（建物・樹木を含む方）をダウンロードして加工し、Web地図上に点群として可視化してみました。マウス等を使ってグリグリと様々な方向から眺めたり、点群を標高に応じて色分けしたりすることも可能になります。
2. 可視化結果
細かい手法の説明は後回しにして、まずは結果を眺めてみましょう。今回は神戸市灘区の王子動物園・神戸高校付近のデータを取得して可視化してみました。
点群は3D地図上の実際の位置・高さに配置されており、標高に応じて色分けされています。
もう少しズームすると以下のようになります。
上半分は山なのですが、1m精度で地表面の凹凸がわかるため、木のもこもこ感が見て取れると思います。 画面中央には平坦な校庭がありますが、よく見ると粒々の人間？が点在していたり、左右に一対のサッカーゴール？があることすらも分かります。
ちなみにGoogle Earthで同じエリアを見ると以下のようになっています。
そして遠くから見ると、、
単純に地図で見るよりも地面の形状がよく分かりますね。
ブラウザでインタラクティブに操作することもできます。 動画↓
1m精度の「建物・樹木などの地物の高さを含む地球表面 [DSM]
— sw1227 (@_sw1227_) 2020年1月11日
データ」を兵庫県が昨日から公開してくれたとのことで、地図上に3次元表示してみた。神戸高校のあたり。
Point cloud layer of DSM(Digital Surface Model) data in Japan.#BuiltWithMapbox and @UberEng ‘s https://t.co/9rOy5jhgB5 pic.twitter.com/1lLfdTSryC
3. 可視化方法
3.1. データの取得
3.1.1. エリア別のダウンロード
データは「G空間情報センター」のWebサイトからダウンロードできます。
兵庫県_全域数値地形図_ポータル（2010年度～2018年度） - データセット
容量が大きいため、ブロックごとにダウンロードする形式となっています。エリアごとのブロック番号は以下の図で確認できます。
兵庫県_全域DSM（2010年度～2018年度） - Indexmap（PDFファイル）
3.1.2. フォーマット
平面直角座標y[m], 平面直角座標x[m], 標高[m]
が連なった形式になっています。
64000.500 -137999.500 81.94 64001.500 -137999.500 81.13 64002.500 -137999.500 81.59 64003.500 -137999.500 82.45
3.1.3. ライセンス
兵庫県オープンデータ利用規約に記載があります。
対象データの利用（複製、改変）については、原則として、クリエイティブ・コモンズ・ライセンス表示 4.0 国際に規定される著作権利用許諾条件によるものとします。
3.2. データの前処理
3.2.1. Pandasによる読み込み
先述のようなスペース区切りのフォーマットは、Python + Pandasで以下のように読み込むことができます。
import pandas as pd
df = pd.read_csv("DSM_05OG701_1g.txt", delimiter=" ", header=None)
df = df.rename(columns={0: "y", 1: "x", 2: "h"})
3.2.2. 緯度経度への変換
座標が平面直角座標系で表現されているため、今回はより一般的な緯度経度に変換します。 平面直角座標と緯度経度の相互変換については実装と考え方をQiita・当ブログに書いたことがあるため、詳しくはそちらを参照してください。
緯度経度と平面直角座標の相互変換をPythonで実装する - Qiita
今回は以下のようなConverterクラスを作成しました。
class Converter():
def __init__(self, phi0_deg, lambda0_deg):
# 定数 (a, F: 世界測地系-測地基準系1980（GRS80）楕円体)
self.m0 = 0.9999
self.a = 6378137.
self.F = 298.257222101
# 平面直角座標系原点をラジアンに直す
self.phi0_rad = np.deg2rad(phi0_deg)
self.lambda0_rad = np.deg2rad(lambda0_deg)
# 補助関数
def A_array(self, n):
A0 = 1 + (n**2)/4. + (n**4)/64.
A1 = - (3./2)*( n - (n**3)/8. - (n**5)/64. )
A2 = (15./16)*( n**2 - (n**4)/4. )
A3 = - (35./48)*( n**3 - (5./16)*(n**5) )
A4 = (315./512)*( n**4 )
A5 = -(693./1280)*( n**5 )
return np.array([A0, A1, A2, A3, A4, A5])
def beta_array(self, n):
b0 = np.nan # dummy
b1 = (1./2)*n - (2./3)*(n**2) + (37./96)*(n**3) - (1./360)*(n**4) - (81./512)*(n**5)
b2 = (1./48)*(n**2) + (1./15)*(n**3) - (437./1440)*(n**4) + (46./105)*(n**5)
b3 = (17./480)*(n**3) - (37./840)*(n**4) - (209./4480)*(n**5)
b4 = (4397./161280)*(n**4) - (11./504)*(n**5)
b5 = (4583./161280)*(n**5)
return np.array([b0, b1, b2, b3, b4, b5])
def delta_array(self, n):
d0 = np.nan # dummy
d1 = 2.*n - (2./3)*(n**2) - 2.*(n**3) + (116./45)*(n**4) + (26./45)*(n**5) - (2854./675)*(n**6)
d2 = (7./3)*(n**2) - (8./5)*(n**3) - (227./45)*(n**4) + (2704./315)*(n**5) + (2323./945)*(n**6)
d3 = (56./15)*(n**3) - (136./35)*(n**4) - (1262./105)*(n**5) + (73814./2835)*(n**6)
d4 = (4279./630)*(n**4) - (332./35)*(n**5) - (399572./14175)*(n**6)
d5 = (4174./315)*(n**5) - (144838./6237)*(n**6)
d6 = (601676./22275)*(n**6)
return np.array([d0, d1, d2, d3, d4, d5, d6])
def convert(self, x, y):
# (1) n, A_i, beta_i, delta_iの計算
n = 1. / (2*self.F - 1)
A_array = self.A_array(n)
beta_array = self.beta_array(n)
delta_array = self.delta_array(n)
# (2), S, Aの計算
A_ = ( (self.m0*self.a)/(1.+n) )*A_array[0]
S_ = ( (self.m0*self.a)/(1.+n) )*( A_array[0]*self.phi0_rad + np.dot(A_array[1:], np.sin(2*self.phi0_rad*np.arange(1,6))) )
# (3) xi, etaの計算
xi = (x + S_) / A_
eta = y / A_
# (4) xi', eta'の計算
xi2 = xi - np.sum(np.multiply(beta_array[1:],
np.multiply(np.sin(2*xi*np.arange(1,6)),
np.cosh(2*eta*np.arange(1,6)))))
eta2 = eta - np.sum(np.multiply(beta_array[1:],
np.multiply(np.cos(2*xi*np.arange(1,6)),
np.sinh(2*eta*np.arange(1,6)))))
# (5) chiの計算
chi = np.arcsin( np.sin(xi2)/np.cosh(eta2) ) # [rad]
latitude = chi + np.dot(delta_array[1:], np.sin(2*chi*np.arange(1, 7))) # [rad]
# (6) 緯度(latitude), 経度(longitude)の計算
longitude = self.lambda0_rad + np.arctan( np.sinh(eta2)/np.cos(xi2) ) # [rad]
# ラジアンを度になおしてreturn
return np.rad2deg(latitude), np.rad2deg(longitude) # [deg]
これを用いて以下のように変換処理を実行し、緯度経度の列を追加します。
ここで平面直角座標系の原点は地域ごとに異なることに注意が必要です。兵庫県はV系に属するため、134度20分0秒, 36度0分0秒が原点となります。以下のリンク（国土地理院）に原点の一覧が記載されています。
converter = Converter(36, 134+20/60)
df["lat"] = df.apply(lambda row: converter.convert(row.x, row.y)[0], axis=1)
df["lon"] = df.apply(lambda row: converter.convert(row.x, row.y)[1], axis=1)
データ量が多いため時間がかかります。何度も行うようなら並列化すると良いでしょう。
3.2.3. CSV出力
これをCSV形式で出力すれば可視化のための前処理は完了です。
df[["h", "lat", "lon"]].to_csv("nada_latlon.csv", index=None)
3.3. MapboxとDeck.glによる3次元地図
緯度経度表記に変換された点群データをブラウザで地図上に可視化するため、以下の技術を用います。
- React
- Mapbox GL JS: ベースマップの表示用
- Deck.gl の PointCloudLayer: 点群の表示用
JSXのシンタックスハイライトが上手くいかないので、コードはQiitaに別途記載しました↓
兵庫県の地表形状（建物等含む）を地図上に可視化【DSM】 - Qiita
4. 最後に
兵庫県全土で公開されている3次元データのうち、特定の1エリアに（データ量の関係で）限定して可視化を行いました。WebGLを活かしたDeck.glのおかげで、手元のMacbook Proでも100万オーダーの点群をリアルタイムに動かすことが可能です。
本来は、地図上で兵庫県内を移動した時に動的に点群をfetchできると嬉しいです。Point Cloudの3D Tilingを行えば良さそうということは分かっているのですが、軽く調べた感じではLAS（LIDAR点群データの形式）からタイリングを行う事例しか見当たらず、いったん断念しています。LASへの変換を自力で実装できれば良いのですが、GeoJSONみたいに軽いノリで生成できる感じでもないので…
300万点でもMacbookで普通に動かせるのすごい#javascript pic.twitter.com/CWmxMCA6hU
— sw1227 (@_sw1227_) 2020年1月11日
関連記事
地理情報プログラミング系の関連記事です。 sw1227.hatenablog.com
