Introduction

　このWEBサイトに設置したリンクは，新しいタブ・ウィンドウで開くことを推奨します．バグの影響で，新しいタブ・ウィンドウでないと開けないリンクが多いためです．予めご容赦ください．

Objectives

　このRPubsは，Rによる空間データの可視化の事例として，建物の3Dプロットを行う方法をメモしたものです．具体事例として，大阪市の御堂筋沿道の建物を3Dで可視化してみます．国土交通省のPLATEAU（プラトー）などの登場により，3次元の属性を持つGIS（オープン）データを用いた空間分析の可能性が広がりつつあります．筆者の知る限り，Rによる3Dデータの空間分析は，他のGISソフトほどは発展していません（処理能力の問題もある気がしますが）．一方，今回取り上げるrayshaderパッケージ等，少しずつ3Dデータを扱うための環境が整い始めています．
　rayshaderパッケージは，ベクタ・ラスタデータの両方の3Dプロットを可能にしたパッケージで，それ単体で用いることができる他，ggplot2パッケージのような既存の可視化パッケージで作成された2Dプロットを3Dにするといった用い方もできます．このチュートリアルでも，そちらの方法をお示ししたいと思います．
　このチュートリアルでは，rayshaderのインストールについては詳しく述べませんが，筆者が調べた限り，インストールが一発成功することは稀なようです．もしうまく行かない場合は，例えばMasahiko Satoさんによるこちらのページを参照し，インストールを行うことをお勧めします．また，インストールの際には，Rのバージョンは可能な限り最新にしておいてください．特に，RStudio上でインストールする際には，native pipeを利用可能な状態にしておく必要があります（詳細はUryu Shinya先生によるこちらのページを例えば参照）．
　なお，このRPubsは，Rを用いたGIS分析に関する基礎的な知識を持っている読者向けのものです．RによるGIS分析については，別途チュートリアルサイトを用意してあるので，そちらを参照ください．

Exercise

Load Packages

　まず，今回の実装例で用いるRパッケージを読み込みます．もし，インストールが済んでいないものがあれば，事前にインストールしておいてください．

library(sf)
library(dplyr)
library(ggplot2)
library(rayshader)

Import data

　今回用いる3Dデータは，国土交通省のPLATEAU（プラトー）のうち，大阪市のデータである3D都市モデル（Project PLATEAU）大阪市（2020年度）です．Rでも扱える形式は，このうちファイルジオデータベース形式ですので，こちらのページからダウンロードを行ってください．データのサイズが大きく，Rでダウンロードするには向かないので，ダウンロードのステップはここでは省略します．
　ダウンロードされたデータは「27100_osaka-shi_fgdb.7z」という圧縮ファイルなので，解凍ソフトで解凍した後得られる「27100_osaka-shi.gdb」を，作業ディレクトリに設置します（このチュートリアルは，Rコードと同じディレクトリにgdbファイル（gdbはGISデータの一種で，詳細はこちらのページを例えば参照）を設置しているので，作業ディレクトリをコード内で明示していません）．gdb内に入っているレイヤは，sfのst_layers関数で表示できます．今回用いるのはそのうち，建物レイヤの「lod0_Building」です．
　gdb内のレイヤの読み込みは，sfのread_sf関数（st_read関数の亜種）で行うことができます．読み込んだ建物データには，各建物ポリゴンのz軸座標が含まれますが，sfではz軸を持つデータを原則扱わないのと，高さを表す変数は別途用意されているので，st_zm関数を用いて座標を削除します．今回可視化するエリアは中央区内だけなので，建物住所の変数「gen_区名」が中央区であるレコードのみ残します．

#gdbファイルに含まれるレイヤの一覧
sf::st_layers(dsn="27100_osaka-shi.gdb")

## Driver: OpenFileGDB 
## Available layers:
##                layer_name    geometry_type features fields
## 1          lod1_TinRelief 3D Unknown (any)      277      3
## 2           lod2_Building 3D Multi Polygon     2967     86
## 3           lod1_Building 3D Multi Polygon   615938     86
## 4           lod0_Building 3D Multi Polygon   615938     88
## 5       lod2_BuildingPart 3D Multi Polygon      467     78
## 6       lod1_BuildingPart 3D Multi Polygon      467     78
## 7        lod2_RoofSurface 3D Multi Polygon     3464      2
## 8        lod2_WallSurface 3D Multi Polygon     3760      2
## 9      lod2_GroundSurface 3D Multi Polygon     3435      2
## 10 lod2_OuterFloorSurface 3D Unknown (any)        1      2
## 11           lod1_LandUse 3D Multi Polygon   775801     17
## 12              lod1_Road 3D Multi Polygon   116826     17

#gdbファイルから建物レイヤを読み込み
lod01 <- sf::read_sf(dsn="27100_osaka-shi.gdb",layer="lod0_Building") %>%
  #z軸を削除
  sf::st_zm() %>%
  #区名が中央区の建物のみ残す
  dplyr::filter(gen_区名=="中央区")

Filter data

　読み込んだ建物データを，更に御堂筋沿道の町丁目のものだけに絞ります．今回は簡単のため，愚直に該当する町丁目名ベクトルを作成し，それを用いてレコードを絞り込みます．また，建物の高さ変数「bldg_measuredHeight」以外の変数は，使わないのでドロップします．

#御堂筋沿道の町丁目のリスト
midosuji <- c("北浜４丁目",
              "北浜３丁目",
              "今橋４丁目",
              "今橋３丁目",
              "高麗橋４丁目",
              "高麗橋３丁目",
              "伏見町４丁目",
              "伏見町３丁目",
              "道修町４丁目",
              "道修町３丁目",
              "平野町４丁目",
              "平野町３丁目",
              "淡路町４丁目",
              "淡路町３丁目",
              "瓦町４丁目",
              "瓦町３丁目",
              "備後町４丁目",
              "備後町３丁目",
              "安土町３丁目",
              "本町３丁目",
              "本町４丁目",
              "南本町４丁目",
              "南本町３丁目",
              "久太郎町４丁目",
              "久太郎町３丁目",
              "北久宝寺町４丁目",
              "北久宝寺町３丁目",
              "南久宝寺町４丁目",
              "南久宝寺町３丁目",
              "博労町４丁目",
              "博労町３丁目",
              "南船場４丁目",
              "南船場３丁目")
lod01_mido <- lod01 %>%
  #上のリストに含まれる町丁目のレコードだけ残す
  dplyr::filter(gen_町丁目名称%in%midosuji) %>%
  #高さ変数だけ残す
  dplyr::select(bldg_measuredHeight)

3D plot

　3Dプロットの元になる2Dの建物高さのプロットを作成します．ggplot2のscale_fill_gradient関数を用いて，建物の高さが低い程黄色に，高い程オレンジ色に近くなるような色分けを行います．

plot_height_2d <- ggplot2::ggplot() +
  #sfオブジェクト可視化の為の設定（変数bldg_measuredHeightで色分け）
  ggplot2::geom_sf(data=lod01_mido,aes(fill=bldg_measuredHeight),lty=0) +
  #色分けの設定（高さが低いほど黄色，高いほどオレンジ色に近い）
  ggplot2::scale_fill_gradient(low="yellow",high="orange",na.value=NA) +
  #プロットのテーマの設定（最低限）
  ggplot2::theme_bw(base_size=5)
plot_height_2d

　ようやく（と行っても手順はかなり少ないですが），rayshaderを用いた3Dプロットを行います．ggplot2で作成した2Dプロットを元に3Dプロットを作成するには，plot_gg関数を用います．引数ggobjには元になる2Dプロットのオブジェクトを指定します．マルチコアが利用できる場合には，引数multicoreをTRUEに設定すると，描画が速くなります．引数widthには（原則）インチ単位のプロットの幅，heightには高さを指定します．引数zoomの値を変えることで，プロットの寄り・引きを調整できます．引数thetaではz軸周りの回転角度を，phiでは方位角（どのくらいの角度から見下ろすか）を指定します．3Dプロットの結果は，ggplot2のプロットとは異なり，（Macであれば）XQuartzのウィンドウが立ち上がって表示されます．作成したプロットは，render_snapshot関数を用いて画像として保存できます．
　ちなみに3Dプロットは，画像形式だけでなく動画形式でも保存できるようです．詳しくは，例えばこちらのWEBページを参照してください．

#上で作成した2Dプロットから3Dプロットを作成
rayshader::plot_gg(ggobj=plot_height_2d,multicore=TRUE,width=5,height=20,zoom=0.4,theta=15,phi=15)
#3Dプロットを画像として保存
rayshader::render_snapshot(filename="plot_height_3d.png",clear=TRUE)

　保存した3Dプロットの画像を表示してみます．百尺制限があった頃の名残か，淀屋橋側（図上側）のスカイラインは大体同じ高さですが，心斎橋側（下側）だと少々高さが不揃いになっていることが見て取れます．

作成された3Dプロット

もしもRで3Dプロットができたなら

Keisuke Takano

8/22/2022