GISコミュニティフォーラム開催報告 : 「ArcGISでデータサイエンスしよう -より高度で自由な地理空間分析へ-」①

1 Kudo

はじめに

2019年 5 月 22 日 (水) ～ 24 日 (金)、東京ミッドタウン六本木にて第15回GIS コミュニティフォーラム及びプレフォーラム・セミナーが開催されました！

公開可能なセッションの発表資料はESRIジャパンのサポートサイトで公開しております。また、GeoNetブログでは開発者向けのセッションで使用したコードなどを順次ご紹介していきます。発表資料や GeoNet記事、GitHub へのリンクなどの詳細は第15回GISコミュニティフォーラム開催報告記事をご覧ください。

テクニカルセッション「ArcGISでデータサイエンスしよう～より高度で自由な地理空間分析へ～」の中では、以下の 3 つのポイントをお伝えいたしました。

ArcGIS はデータサイエンスの様々な場面で活用できる
既に多くの統計及び機械学習の手法がジオプロセシングツールとして実装されている
既存のツールでは対応が難しい場合においても、オープンなプラットフォームである特徴を生かして R や Python の各種ライブラリ等と連携することで対応できる

特に ArcGIS の既存のツールでは対応が難しいケースとして、空間的自己相関を持つデータの回帰分析や、ディープラーニングを用いた物体検知を例に挙げ、それぞれ R およびArcGIS API for Python を用いてデモと共にご紹介しましたが、本ブログの中でセッション中ではご紹介しきれなかった具体的なコードの中身について詳しくご紹介します。

本稿では R との連携についてご紹介し、ArcGIS API for Python を用いたディープラーニングは別記事にてご紹介します。

概要

本稿は以下の流れに沿ってご紹介します。

R と ArcGIS Pro を連携させる R-ArcGIS Bridgeについて
環境構築の手順
テクニカルセッションでご紹介したデモの詳細
まとめ

R と ArcGIS Pro を連携させる R-ArcGIS Bridgeについて

セッション中でもご紹介しましたが、R-ArcGIS Bridge というツールを用いることで、R と ArcGIS Pro を連携させることが可能になります。R-ArcGIS Bridge は米国 Esri 社の GitHub 上で Apache License 2.0 に基づいて公開されているツールです。国内でのサポート等はしておりませんが、商用利用も可能な無償のツールとして提供されており、ArcGIS Pro 1.1 以降及び R 3.2.2 以降のバージョンに対応しています。

このツールが提供する機能は大きく分けて以下の 3 つです。

R と ArcGIS Pro 間でのデータを読み込み／書き込みする機能
R で用いられるデータフォーマットへ変換する機能
R をラッピングしてスクリプトツールにするための機能

本稿では 3 番目の機能については紹介しませんが、Python からスクリプトツールを作成できるように、R でも ArcGIS Pro 上で動作するスクリプトツールを作成することが可能です。詳細は GitHub 上でのチュートリアル (英語) を参照ください。

環境構築の手順

以下では実際に R-ArcGIS Bridge を使用する際の環境構築の手順を、予めArcGIS Pro 及び R がインストールされていることを想定してご紹介します。

ArcGIS Proのインストール方法については弊社サポートサイト内のインストールガイドを、R のインストール方法ついてはオンライン上の各種ドキュメントをご参照下さい。インストーラーは CRAN (Rや各種パッケージをダウンロードするための Web サイト) から入手可能で、国内では統計数理研究所様や山形大学様がミラーサイトを運営しています。

セッションで利用したバージョンは ArcGIS Pro 2.3.2 及び R 3.5.9 です。

※ オプションとして、統合開発環境 (IDE) を利用することで、R でのコーディングが容易になります。R で非常によく利用される IDE として、RStudio が挙げられます。他にVisual Studio や Jupyter Notebook 等でも実行可能です。

R-ArcGIS Bridge は、実際は arcgisbinding という名称の R のパッケージです。ArcGIS Pro 2.3.x では以下の手順でインストールすることができます。

ArcGIS Pro を開き、左下の [設定] をクリックします

[オプション] をクリックし、ダイアログボックスを開きます
[ジオプロセシング] をクリックし、[R-ArcGIS サポート] の項目の中の [検出された R ホームディレクトリ] にインストールした R のディレクトリが設定されていることを確認してください。ディレクトリが異なる場合は手動で変更してください。
[ArcGIS R 統合パッケージをインストールしてください] と記載されている右側のアイコンをクリックし、[インターネットからのパッケージのインストール] をクリックするとインストールが始まります。

別のインストールオプションとして、ファイルを指定してインストールすることも可能です。その場合は [最新バージョンのダウンロード...] をクリックし、リンク先から ZIP ファイルをダウンロード後、[ファイルからのインストール...] をクリックして、ダイアログに従ってダウンロードした ZIP ファイルを選択してください。

テクニカルセッションでご紹介したデモの詳細

本セッションにおいては、ボストンの各地区の住宅価格の中央値を、各地区の犯罪発生率や大気汚染の度合い等の変数を用いて予測する解析を、以下の 2 つのパターンで実行し結果を比較しました。

ArcGIS Pro の [一般化線形回帰分析] を使用した場合
R の spdep パッケージを用いて空間的自己相関を考慮した統計モデルを使用した場合

ここからは、セッションの中では詳しく紹介できなかった、R での解析について具体的なコードと共にご紹介します。

パッケージのインストール

R では多数のパッケージが公開されており、それらを利用することで機能拡張が可能です。デモでは以下のパッケージを利用しました。必要に応じてコメントアウトしてパッケージをインストールしてください。

# 必要に応じてパッケージをインストール=========================================================
# install.packages("spData", dependencies = TRUE)
# install.packages("sf", dependencies = TRUE)
# install.packages("tidyverse", dependencies = TRUE)
# install.packages("spdep", dependencies = TRUE)
# install.packages("spatialreg", dependencies = TRUE)
# install.packages("ggthemes", dependencies = TRUE)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

ArcGIS Pro へのデータの書き込み

上記でインストールした spData パッケージは多数の空間データを含んでいます。今回は spData の中にあるボストンの住宅価格のデータ (シェープファイル) を使用します。データの読み込みには sf パッケージの read_sf 関数を使用しています。

このデータを R から ArcGIS Pro へデータを書き込むために、R-ArcGIS Bridge (arcgisbinding) を読み込み、7 行目では arc.check_product 関数でライセンスのチェックを行っています。arc.write 関数を使用し、予め用意してあるジオデータベースへ bstn_house_price という名前のフィーチャクラスとして書き込んでいます。

# spDataパッケージのボストンの住宅価格のデータを ArcGIS へ渡す=================================
require(sf)
bstn_shp <- read_sf(system.file("shapes/boston_tracts.shp", package = "spData"))
head(bstn_shp)

# arcgisbindingをインポート
require(arcgisbinding)
arc.check_product()

# プロジェクトのFGBDへ書き込む
arc.write(path = "gcf2019_sem/gcf2019_sem.gdb/bstn_house_price", data = bstn_shp)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

ArcGIS Pro へのデータの読み込み

ArcGIS Pro に書き込んだフィーチャクラスのデータを、arc.open 関数を使って再度 R へ読み込みましょう。読み込んだデータは、arc.select 関数を用いて arc.data というデータフレーム形式のオブジェクトに変換できます。また、その際にフィールド名やSQL 形式の構文を引数に渡して必要なデータだけ抽出することができます。今回は住宅価格の中央値 (CMEDV)、犯罪率 (CRIM)、大気汚染の度合い (NOX) 等が格納されているフィールドだけ抽出しました。

# ArcGIS Proのデータを読み込む===============================================================
d <- arc.open("gcf2019_sem/gcf2019_sem.gdb/bstn_house_price")
bstn_fc <- arc.select(d, c("OBJECTID", "Shape", "CMEDV", "CRIM", "NOX", "RM", "DIS", "LSTAT", "BB", "POP"))
head(bstn_fc)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

R の spatilreg パッケージを使用した解析

データの準備ができたので、解析を実行します。今回は空間データの解析に係る関数を多数提供している spedep 及び spatialreg パッケージを利用します。

以下のコードでは、poly2nb 関数及び nb2listw 関数を使用して、フィーチャクラス内の各ポリゴンの隣接関係を重み付きの行列として定義しています。

その後、spatialreg パッケージの errorsarlm 関数を使用して空間的自己相関を加味した統計モデルの解析を実行します。引数の formula = CMEDV ~ CRIM + NOX + RM + DIS + LSTAT + BB + POP は、CMEDV を目的変数、その他のフィールドを説明変数としていることを意味します。

また、summary 関数で実行した結果のサマリーを確認することが可能です。

# Rのspdepパッケージを使って空間誤差モデルの解析をする=======================================
bstn_sp <- arc.data2sp(bstn_fc)

require(spdep)
W <- poly2nb(bstn_sp, queen = TRUE) %>% 
    nb2listw(style = "W", zero.policy = TRUE)

require(spatialreg)
SEM <- errorsarlm(
    formula = CMEDV ~ CRIM + NOX + RM + DIS + LSTAT + BB + POP,
    listw = W,
    data = bstn_sp
)
summary(SEM)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

解析結果を元のデータに付与

解析結果を代入したオブジェクト (SEM) の中には residual (残差)、fitted.value (推定結果) が含まれています。以下のコードで、それらを元のデータに付与しています。

また、ここではパイプ演算子 (%>%) を使用しています。この演算子は関数等の実行結果を次の関数の第一引数に渡してくれるため、複数の処理を連続的に記述できます。中間処理結果を代入するための一時的な変数を作成する必要もなくなるため、コードの可読性が向上します。

# 解析結果を元のデータに付与する=================================================================
sem_values <- bstn_sp@data %>% 
    cbind(residuals = SEM$residuals) %>% 
    cbind(fitted = SEM$fitted.values)
bstn_sp@data <- (sem_values)
head(bstn_sp@data)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

補足: R でのチャート作成

セッション中では詳しく紹介しませんでしたが、R の利用に慣れている場合は ggplot2 というパッケージを利用してチャートを作成しても良いでしょう。

以下のコードは 2 種類のチャートを作成しています。1 つ目は横軸に予測結果、縦軸に残差を取ったものです。プロットされた点が 0 を示す赤い点線から離れているほど残差が大きい (予測と実際の値の差が大きい）ことを示します。2 つ目は縦軸に予測結果を、横軸に実際の値を取っています。こちらも同じように赤い点線から離れたポイントは予測と実際の値の差が大きいことを示しています。

# 結果の可視化===============================================================================
require(ggplot2)
require(ggthemes)
g <- ggplot(bstn_sp@data, aes(x = fitted, y = residuals, ymin = -30, ymax = 30, xmin = -10, xmax = 50))
g <- g + theme_hc(base_size = 16)
g <- g + coord_fixed(ratio=1)
g <- g + geom_point(colour = "#003366", alpha = 0.4)
g <- g + geom_abline(aes(slope=0, intercept=0), color='red', linetype = "31")
g <- g + labs(x = "Predicted", y = "Residuals")
g

g <- ggplot(bstn_sp@data, aes(x = CMEDV, y = fitted, ymin = 0, ymax = 55, xmin = 0, xmax = 55))
g <- g + theme_hc(base_size = 16)
g <- g + coord_fixed(ratio=1)
g <- g + geom_point(colour = "#003366", alpha = 0.4)
g <- g + geom_abline(aes(slope=1, intercept=0), color='red', linetype = "31")
g <- g + labs(x = "Observed", y = "Predicted")
g‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

結果を ArcGIS Pro へ書き込む

最後に、冒頭と同じように arc.write 関数を使用して、解析結果を付与したデータを新たなフィーチャクラスとしてファイルジオデータベースに書き込みます。

# データをArcGIS Proへ戻す==================================================================
arc.write(bstn_sp, path = "gcf2019_sem/gcf2019_sem.gdb/sem_rslt", overwrite = TRUE)‍‍‍‍‍‍‍‍‍‍

地図上での可視化

ArcGIS Pro 上で残差を地図上にプロットします。以下の画像では、ArcGIS Pro の [一般化線形モデル] の結果と、R での解析結果を交互に表示しています。両者の結果を比較できるようにシンボルを調整していますが、どちらも残差を等間隔で分類し塗り分けています。R での解析結果の方が色の濃い地区が減少しており、残差のばらつきが小さい (予測と実際の値の差が少ない) ことがわかります。

また、補足として R を用いたチャート作成のコードをご紹介しましたが、ArcGIS Pro にはノンコーディングで見栄えの良いチャートを作成する機能が標準で実装されています。以下では ArcGIS Pro での解析結果と、R を用いた場合の解析結果、それぞれの残差をプロットして比較しています。下側に表示されている R での解析結果の方が、残差が 0 付近に集まっているのがわかります。

作成したチャートはマップと連動しているため、チャート上でポイントを選択するとマップ上でも同じように選択され、マップ上でもすぐに確認することが可能です。

まとめ

さて、本稿ではテクニカルセッション中にご紹介した ArcGIS Pro と R の連携のデモンストレーションついて主に R のコードを中心に紹介しました。ArcGIS の既存のツールでの解析が難しい場合などは R の様々なパッケージを利用することで対応できる範囲を拡張することが可能です。

一方で、R 側で空間的な処理を全て行うよりも、ArcGIS Pro のジオプロセシングツールを利用した方がずっと手軽なケースも多いでしょう。また、R 単体では解析結果を表示する Web アプリを作成し、組織内外に共有するためには別途開発が必要です。Web GIS、特に ArcGIS Online を活用すればプラットフォームの利点を活かしてすぐに解析結果の共有が可能です。

本記事で紹介したコードは ESRIジャパンの GitHub 上で公開しています。是非ご参考にしていただき、双方の強みを組み合わせて、さらに ArcGIS をご活用いただければ幸いです。

※ ArcGIS API for Python を用いたディープラーニングは、別途本ブログにてご紹介いたします。

関連リンク