『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門：空間解析・可視化編
Keisuke Takano
はじめに
チュートリアルの概要
このチュートリアルは，『事例で学ぶ経済・政策分析のためのGIS入門』中で用いられているRによる空間解析・可視化をより深く理解する為のものです．本文や付録内では，空間解析の為のパッケージとしてsfパッケージとdplyrパッケージを，可視化のためのパッケージとしてggplot2パッケージとmapviewパッケージを頻繁に用いています．
sfパッケージはベクタデータ（ポイント・ポリゴン・ラインデータ）の操作を行う際に便利なパッケージで，dplyrパッケージは空間データに限らない様々なデータの操作に用いられるパッケージです．また，ggplot2パッケージはRで扱われる様々なデータの可視化に用いられるパッケージで，mapviewパッケージは空間データのインタラクティブな参照を行うことが可能な空間可視化のパッケージです．ここでは，それらの基礎的な使い方を，WEB上で入手可能な空間データを用いた演習を通じて解説していきます．各パッケージの詳細については，以下の公式WEBページを参照してください．
dplyrパッケージ及びデータ構築に必要なR関数の基礎については，同じく『事例で学ぶ経済・政策分析のためのGIS入門』の付録に含まれている「『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門：基礎編」で解説しているので，R言語自体に馴染みがない方は，まずそちらで演習を行うことをお勧めします．
また，GISの基礎知識については，例えば本書『事例で学ぶ経済・政策分析のためのGIS入門』の第1章や「GISオープン教材」の「GISの基本概念」セクション，矢野桂司先生による『GIS地理情報システム』のChapter2等で学んだ上で，演習に取り組まれることをお勧めします．
このチュートリアルの後半では，重回帰分析や固定効果法，差分の差分法といった，計量経済学的手法を用いた統計分析を行います．それら手法の基礎知識について，例えば田中隆一先生による『計量経済学の第一歩』であらかじめ学んでおくと，理解がより深まります．
さらなる学習のために（このチュートリアルの参考文献集）
RによるGIS
sfパッケージとdplyrパッケージを組み合わせたGISに関して，さらに詳細・高度なジオプロセシングを解説しているWEBサイトをいくつか紹介します（作成者名の敬称略）．日本語のサイトとしては，以下がおすすめです．
- Rを使った地理空間データの可視化と分析
(S. Uryu)
- Rで行う公示地価パネルデータの整備と分析 (S. Kuroda)
- もしもArcGISのジオプロセシングをRで実装したら (K. Takano)
このチュートリアルは「もしも」をベースとして同じ筆者によって作成されたものですが，筆者によるこれらチュートリアルの当初のコンセプトは，「Rで行う公示地価パネルデータの整備と分析」の内容を，扱う空間データ操作の種類を増やしたり，dplyrで処理を記述する等によって補完するというものでした．英語のサイトとしては，以下がおすすめです．
- Geocomputation with
R (R. Lovelace, J. Nowosad, J. Muenchow)
- Introduction to
Spatial Data Programming with R (M. Dorman)
- R as
GIS for Economists (T. Mieno)
- r-spatial (E. Pebesma, M.
Appel, D. Nüst)
- Spatial Data
Science (E. Pebesma, R. Bivand)
- Using Spatial Data with R (C.A. Engel)
GISを用いた経済分析
社会経済データの分析に特化したGISのチュートリアルとしては，いずれもArcGIS使用のものですが，以下がおすすめです．
- オンラインGIS教材（『経済・政策分析のためのGIS入門』シリーズの補助教材） (M. Kawabata)
- Causal Inference with Spatial Data: ArcGIS 10 for Economics Research (M. Kudamatsu)
また，経済学におけるGISの活用事例については，下松真之先生による論文GIS for Credible Identification Strategies in Economics Researchが包括的にまとめていますので，そちらを学ぶことをお勧めします．
計量経済分析
このチュートリアルの後半では，Rを用いたいくつかの統計分析を行います．その前提知識を学ぶ上では，以下の教科書が入門編として最適です．
- 計量経済学の第一歩 (R. Tanaka)
計量経済学に関する発展的な内容を学びたい場合，以下の教科書が有用です．
また，現状英語版のみが利用可能ですが，以下の教科書はWEB上で無料で読むことができます．
- Causal Inference: The
Mixtape (S. Cunningham)
- Using R for Introductory Econometrics (F. Heiss)
加えて，このチュートリアルの主題である交通（・都市）に関する計量経済分析一般についてさらに詳しく学びたい場合は，以下の教科書・論文がおすすめです．
- 交通インフラ整備効果の因果推論：論点整理と展望 (T. Otazawa, Y. Ohira)
- Transportation Costs
and the Spatial Organization of Economic Activity (S. Redding, M.
Turner)
- An Introduction to Geographical and Urban Economics (S. Brakman, H. Garretsen, C. Van Marrewijk)
特に，このチュートリアルのパネルデータ分析に関する解説は，An Introduction to Geographical and Urban EconomicsのChapter 3に依拠して書かれています．
sfパッケージの概要
このセクションでは，sfパッケージの特徴について，主に，University of LeedsのRobin Lovelace先生らによって書かれた，Rによる空間解析の包括的なブックであるGeocomputation with Rに依拠しつつ，概説します．現状ブックは英語版のみ利用可能ですが，sfパッケージに関してより詳細な事項をもし学びたい場合には，そちらを参照することをお勧めします．
sfパッケージの魅力
Geocomputation with Rでは，sfパッケージの主たる魅力として，以下の5つを挙げています．
- データの高速読み取り・書き込みが可能
- 強化されたプロット・パフォーマンス
- ほとんどの操作でデータフレームとして扱うことが可能
- 関数名は比較的一貫性があり，直感に合う
- dplyrパッケージを含むtidyverseパッケージとの親和性が高い（パイプ演算子%>%が使える）
このチュートリアルでも頻繁に，sfパッケージによる空間データの作成と，dplyrパッケージによるデータフレーム操作を行き来しながら分析を進めていきますので，dplyrパッケージと両輪でsfパッケージを学ぶことが極めて重要となります．
sfパッケージで扱われる空間データ
sfパッケージで扱われる空間データのオブジェクトの種類や構造について概説します．「sf」は「simple features」の略であり，Rのsfパッケージでは7種類の形式が扱われます．各々1つのポイント・ライン・ポリゴンデータでなる形式（POINT, LINESTRING, POLYGON），それぞれの「マルチ」バージョン（同じタイプのフィーチャを1つのフィーチャにグループ化したもので，それぞれMULTIPOINT, MULTILINESTRING, MULTIPOLYGON），それらを全てひとまとめにしたGEOMETRYCOLLECTIONの合計7種類です（Figure 1を参照）．
Figure 1: Simple
featuresの形式
sfオブジェクトの構造を示したのが，下のFigure
2です．まず，地物とその形状を，座標値を元にポイント・ライン・ポリゴンで表したものがsfgクラスと呼ばれるものです．この状態では，空間データは座標系に関する情報も，地物の属性に関する情報も含んでいません．既存のshapefileをsfオブジェクトに変換して扱う際には，分析者が自身でsfgを作成することはありませんが，時には座標値を元にsfオブジェクトを作成する場面もあります．
sfgクラスに，座標系の情報を付与した（あるいは付与可能な状態にした）のがsfcオブジェクトで，更にデータフレーム等の形式で格納された地物の属性を付与したのが，我々がこのチュートリアルで主に扱うsfオブジェクトです．
Figure 2: sfオブジェクトの構造
使用データ
今回は，国土数値情報などから取得可能な下に示すデータセットを用いて，Rでのジオプロセシングの実践方法を紹介していきます．なお，分析対象の地域は福岡県福岡市です．各データセット内の変数の定義については，データ名のリンクをクリックすると確認できます．これらデータは，『事例で学ぶ経済・政策分析のためのGIS入門』サポートページからダウンロード可能です．
空間データ
- 福岡県の地価データ（評価地点のポイントデータ）：L01-22_40
- 鉄道路線データ（駅のポイントデータ）：N05-20_Station2
- 福岡県の小地域データ（国勢調査の調査単位である小地域＝町丁目のポリゴンデータ）：h27ka40
- 平成27年国勢調査町丁・字等別境界データ
- 統計GIS（総務省）からダウンロード．データの概要はこちらのページにある定義書を参照．
- 平成27年国勢調査町丁・字等別境界データ
テーブル
- 店舗データ（OpenStreetMapから取得された緯度経度座標付き店舗リスト）：fuk_p_lonlat
- Rのosmdataパッケージを用いてダウンロードされたデータをエクセルファイルに変換したもの．
- ダウンロード・変換プロセスは，別添の「intro_spat_vis_app.Rmd」を参照．
- ダウンロード・変換プロセスは，別添の「intro_spat_vis_app.Rmd」を参照．
- データの仕様はOpenStreetMap
Wikiを参照．
- データに含まれる変数の定義は以下の通り．緯度・経度の座標系はJGD2011で定義．
- osm_id：OSM内でのオブジェクトID
- name：店舗名
- shop：店舗分類（分類の詳細は，OpenStreetMap
Wikiで確認可能）
- KEN_NAME：都道府県名
- GST_NAME：市区町村名
- CSS_NAME：区名
- MOJI：小地域名
- lon：経度
- lat：緯度
- osm_id：OSM内でのオブジェクトID
- 現時点では，福岡市周辺店舗の網羅率は未だ高くなく，網羅範囲もチェーン店に偏っている点に注意が必要．
- Rのosmdataパッケージを用いてダウンロードされたデータをエクセルファイルに変換したもの．
演習
パッケージのインストール
まず，上で挙げた4つのパッケージをインストールするコードは以下のようになります．なお，このチュートリアルの筆者のPCでは，これらは全てインストール済みなので，Rmdファイル実行の都合上プログラムとしては実行しません（故に，コメント化されています）．同様に，既にこれらパッケージのインストールが完了する場合は実行する必要はありません．インストールが済んでいない方は，「#」を消去した上でプログラムを実行し，インストールを完了してください．
dplyr・ggplot2パッケージは，tidyverseパッケージの中に含まれていますので，tidyverseひとつをインストールすることで，一括でインストールが可能です．
# #tidyverseパッケージをインストール
# install.packages("tidyverse")
# #sfパッケージをインストール
# install.packages("sf")
# #mapviewパッケージをインストール
# install.packages("mapview")続いて，以下に挙げるその他のパッケージをインストールします．
- openxlsx：Rでxlsx形式のファイルを読み込み・書き出しするためのパッケージ
- ggspatial：ggplot2を用いて地図を作成する際に補助的に用いるパッケージ
- RColorBrewer：様々なプロットで使えるカラーパレットを収録したパッケージ
- fixest：固定効果モデルを高速推定するためのパッケージ
- modelsummary：回帰分析の結果を綺麗に表示するためのパッケージ
- nngeo：k近傍（k番目に近い）の空間データを抽出するためのパッケージ
# #openxlsxパッケージをインストール
# install.packages("openxlsx")
# #ggspatialパッケージをインストール
# install.packages("ggspatial")
# #RColorBrewerパッケージをインストール
# install.packages("RColorBrewer")
# #fixestパッケージをインストール
# install.packages("fixest")
# #modelsummaryパッケージをインストール
# install.packages("modelsummary")
# #nngeoパッケージをインストール
# install.packages("nngeo")インストールの過程で「Do you want to install from sources the packages which need compilation? (Yes/no/cancel)」というメッセージが表示された場合には，とりあえず「no」を選択することをお勧めします．
パッケージの読み込み
上でインストールしたパッケージをRに読み込みます．インストールは一度行えば良いのですが，パッケージの読み込みは，Rを起動するたびに行う必要があります．
#tidyverseパッケージを読み込み
library(tidyverse)
#sfパッケージを読み込み
library(sf)
#mapviewパッケージを読み込み
library(mapview)
#環境によって，プロットがうまくいかない場合は以下のコードを実行
mapview::mapviewOptions(fgb=FALSE)
#openxlsxパッケージを読み込み
library(openxlsx)
#ggspatialパッケージを読み込み
library(ggspatial)
#RColorBrewerパッケージを読み込み
library(RColorBrewer)
#fixestパッケージを読み込み
library(fixest)
#modelsummaryパッケージを読み込み
library(modelsummary)
#nngeoパッケージを読み込み
library(nngeo)動作環境によっては，sfパッケージの読み込み時にエラーが出る場合があります．その原因は，sfパッケージが依存しているRcppのバージョンが古いことが多いようです．もしエラーが出る場合には，以下のコードを実行して見てください（筆者のPCではエラーが出なかったので，コメント化してあります）．
# #Rcppパッケージをインストールする
# install.packages("Rcpp")
# #Rcppパッケージを読み込み
# library(Rcpp)データの読み込み
先ほどリストした，空間データとテーブルを読み込みます．空間データを読み込む際の最も基本的な関数は，sfのst_read関数です．今回は差し当たり，読み込む空間データはこのRmdファイルと同じディレクトリ・階層に設置してください．引数dsnには読み込む空間データのディレクトリ（Rmdファイルと同じディレクトリに設置されていれば，ファイル名を指定するだけでOK）です．引数stringsAsFactorsは，読み込まれた空間データの地物属性（データフレームとして読み込まれます）に含まれる文字列変数を因子型として扱うかを指定するものですが，今回はFALSEを指定します．引数optionsは，データ読み込み時の諸々のオプションを指定するもので，今回は読み込む空間データの地物属性の文字コードを指定します．引数quietは，読み込み結果の詳細を表示しないようにするもので，今回はTRUEを指定します．
読み込んだ空間データの最後の列に「geometry」という変数があります．この変数は，データ内の各レコードが該当するsimple
featuresの形式や，その座標値を表すものです．形式の後に座標値が続くこのような表記形式はWell-known
text（略してWKT）形式と呼ばれます．このチュートリアルでは以後，変数geometryをgeometry属性とも表記します．
##地価（福岡県）
lp_fukuoka <- sf::st_read(dsn="L01-22_40.shp",
stringsAsFactors=F,
options="ENCODING=CP932",
quiet=TRUE)
#先頭行を表示する
head(lp_fukuoka)## Simple feature collection with 6 features and 139 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 130.9289 ymin: 33.87318 xmax: 131.0075 ymax: 33.94818
## Geodetic CRS: JGD2000
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## 1 000 001 000 001 2022 92300 1.7 1 false
## 2 000 002 000 002 2022 59500 -1.2 1 false
## 3 000 003 000 003 2022 50600 -0.8 1 false
## 4 000 004 000 004 2022 75000 0.4 1 false
## 5 000 005 000 005 2022 26400 -1.1 1 false
## 6 000 006 000 006 2022 33200 -0.6 1 false
## L01_010 L01_011 L01_012 L01_013 L01_014 L01_015 L01_016 L01_017 L01_018
## 1 false false false false false false false false false
## 2 false false false false false false false false false
## 3 false false false false false false false false false
## 4 false false false false false false false false false
## 5 false false false false false false false false false
## 6 false false false false false false false false false
## L01_019 L01_020 L01_021 L01_022 L01_023
## 1 false false false 40101 北九州門司
## 2 false false false 40101 北九州門司
## 3 false false false 40101 北九州門司
## 4 false false false 40101 北九州門司
## 5 false false false 40101 北九州門司
## 6 false false false 40101 北九州門司
## L01_024 L01_025 L01_026
## 1 福岡県 北九州市門司区上馬寄１丁目２番２５ 上馬寄１−２−１６ 223
## 2 福岡県 北九州市門司区東門司２丁目１９５４番３ 東門司２−３−１８ 109
## 3 福岡県 北九州市門司区清滝３丁目１２番１３ 清滝３−１−５ 118
## 4 福岡県 北九州市門司区黄金町３１番８ 黄金町７−２０ 293
## 5 福岡県 北九州市門司区白野江４丁目１９８１番２ 白野江４−２−１５ 219
## 6 福岡県 北九州市門司区大字畑字今在家２１４４番１ _ 231
## L01_027 L01_028 L01_029 L01_030 L01_031 L01_032 L01_033 L01_034 L01_035
## 1 住宅 _ 001 W true true true _ 1.0
## 2 住宅 _ 001 W true true true _ 1.5
## 3 住宅 _ 001 W true true true _ 1.0
## 4 住宅 _ 001 W true true true _ 1.2
## 5 住宅 _ 001 W true false true _ 1.0
## 6 住宅 _ 001 W true false true 台形 1.0
## L01_036 L01_037 L01_038 L01_039 L01_040 L01_041 L01_042 L01_043 L01_044
## 1 1.5 2 0 市道 西 5.0 _ _ _
## 2 1.0 2 0 市道 南 6.0 _ _ _
## 3 1.2 2 0 市道 北 4.5 _ _ _
## 4 1.0 2 0 市道 東 6.0 _ _ _
## 5 1.5 2 0 県道 東 11.0 _ _ _
## 6 1.5 2 0 市道 西 4.6 _ _ _
## L01_045 L01_046 L01_047 L01_048 L01_049
## 1 _ _ 中規模一般住宅が区画整然と建ち並ぶ住宅地域 門司 1900
## 2 _ _ 一般住宅の中に小売店舗も見られる住宅地域 門司港 1400
## 3 _ _ 小規模一般住宅が多い閑静な住宅地域 門司港 900
## 4 _ _ 一般住宅のほかにマンションも多い住宅地域 門司 1100
## 5 _ _ 中規模一般住宅、農家住宅が混在する住宅地域 門司港 7100
## 6 _ _ 中規模一般住宅が建ち並ぶ住宅地域 門司 7500
## L01_050 L01_051 L01_052 L01_053 L01_054 L01_055 L01_056 L01_057 L01_058
## 1 1中専 _ 市街化 0 _ _ 60 200 false
## 2 1住居 準防 市街化 0 _ _ 60 200 false
## 3 1住居 _ 市街化 0 _ _ 60 200 false
## 4 1住居 準防 市街化 0 _ _ 60 200 false
## 5 1低専 _ 市街化 0 _ _ 50 80 false
## 6 1中専 _ 市街化 0 _ _ 60 200 false
## L01_059 L01_060 L01_061 L01_062 L01_063
## 1 true 0000000000000000000000000000001111111111 0 0 0
## 2 false 0000000011111111111111111111111111111111 0 0 0
## 3 false 1111111111111111111111111111111111111111 70600 73500 75500
## 4 false 0000000000001111111111111111111111111111 0 0 0
## 5 false 1111111111111111111111111111111111111111 40000 41500 42500
## 6 false 1111111111111111111111111111111111111111 45000 46500 47600
## L01_064 L01_065 L01_066 L01_067 L01_068 L01_069 L01_070 L01_071 L01_072
## 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 83000 85500 85800 85800
## 3 76500 76500 76500 77300 78100 82500 85000 85500 85500
## 4 0 0 0 0 0 0 0 0 0
## 5 43000 43000 43000 43000 43500 44800 46500 46600 46600
## 6 48000 48000 48700 49700 50500 56000 58000 59000 59000
## L01_073 L01_074 L01_075 L01_076 L01_077 L01_078 L01_079 L01_080 L01_081
## 1 0 0 0 0 0 0 0 0 0
## 2 85800 85800 85800 85800 85800 85800 84700 82800 80500
## 3 85500 85500 85500 85500 85500 85500 84700 84300 82600
## 4 110000 110000 110000 110000 110000 108000 106000 104000 102000
## 5 46600 46600 46600 46600 46600 46600 46200 45800 45300
## 6 59500 60100 60600 61000 61500 61500 60700 60000 58900
## L01_082 L01_083 L01_084 L01_085 L01_086 L01_087 L01_088 L01_089 L01_090
## 1 0 0 0 0 0 0 0 0 0
## 2 78000 75300 72400 69600 67100 65100 62800 61400 60900
## 3 79200 75600 72000 68500 66500 64800 62800 61000 59100
## 4 98000 94300 91000 87200 84100 81200 78400 76900 76000
## 5 44500 42500 40600 38600 36900 35200 33400 31900 30800
## 6 53500 50200 47200 44500 42900 41500 40000 38700 37500
## L01_091 L01_092 L01_093 L01_094 L01_095 L01_096 L01_097 L01_098 L01_099
## 1 84500 84500 84500 84800 85200 85400 87100 89900 90800
## 2 60600 60300 60000 59900 59900 59900 59900 60200 60200
## 3 57300 55900 54700 53800 53100 52400 51800 51400 51000
## 4 75200 74800 74500 74300 74300 74300 74300 74700 74700
## 5 29800 29100 28600 28200 27900 27600 27300 27000 26700
## 6 36300 35500 34800 34300 34000 33800 33600 33500 33400
## L01_100 L01_101 L01_102 L01_103 L01_104
## 1 92300 00000000000000 00000000000000 00000000000000 00000000000000
## 2 59500 00000000000000 00000000000000 00000000000000 00000000000000
## 3 50600 10000000000000 10000000000000 20000000000000 10000000000000
## 4 75000 00000000000000 00000000000000 00000000000000 00000000000000
## 5 26400 10000000000000 10011000000000 10000000000000 10000000000000
## 6 33200 10000000000000 10000000000000 10000000000000 10000000000000
## L01_105 L01_106 L01_107 L01_108 L01_109
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 00000000000000 00000000000000 00000000000000 40000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 5 11000000000000 10000000000000 10000000000000 10000000000000 10000100000000
## 6 10000100000010 10000000000000 10000100000000 10000000000000 10000000000000
## L01_110 L01_111 L01_112 L01_113 L01_114
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000001000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000001000000
## 4 00000000000000 00000000000000 40000000000000 10000000000000 10000001000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000001000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000001000010
## L01_115 L01_116 L01_117 L01_118 L01_119
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000010000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000010000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000001000010 10000000000000 10000000000000
## L01_120 L01_121 L01_122 L01_123 L01_124
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000010 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000010 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000010 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000011 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_125 L01_126 L01_127 L01_128 L01_129
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 20000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_130 L01_131 L01_132 L01_133 L01_134
## 1 40000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_135 L01_136 L01_137 L01_138 L01_139
## 1 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## geometry
## 1 POINT (130.9289 33.89107)
## 2 POINT (130.9727 33.94818)
## 3 POINT (130.9612 33.94013)
## 4 POINT (130.9405 33.90936)
## 5 POINT (131.0075 33.94478)
## 6 POINT (130.9751 33.87318)##鉄道時系列（全国・駅データ）
rail_jp <- sf::st_read(dsn="N05-20_Station2.shp",
stringsAsFactors=F,
options="ENCODING=CP932",
quiet=TRUE)
head(rail_jp)## Simple feature collection with 6 features and 11 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 141.9359 ymin: 43.62763 xmax: 142.3579 ymax: 43.78933
## Geodetic CRS: JGD2011
## N05_001 N05_002 N05_003 N05_004 N05_005b N05_005e
## 1 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 2 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 3 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 4 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 5 2 函館線 北海道旅客鉄道（旧国鉄） 1911 1950 9999
## 6 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## N05_006 N05_007 N05_008 N05_009 N05_011 geometry
## 1 EB03_11218106 <NA> <NA> <NA> 旭川 POINT (142.3579 43.76344)
## 2 EB03_11218099 <NA> <NA> <NA> 妹背牛 POINT (141.9666 43.69096)
## 3 EB03_11218103 <NA> <NA> <NA> 伊納 POINT (142.2717 43.7632)
## 4 EB03_11218100 <NA> <NA> <NA> 深川 POINT (142.0414 43.72125)
## 5 EB03_11218105 <NA> <NA> <NA> 近文 POINT (142.3253 43.78933)
## 6 EB03_11218098 <NA> <NA> <NA> 江部乙 POINT (141.9359 43.62763)##小地域（福岡県）
ct_fukuoka <- sf::st_read(dsn="h27ka40.shp",
stringsAsFactors=F,
options="ENCODING=CP932",
quiet=TRUE)
head(ct_fukuoka)## Simple feature collection with 6 features and 35 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 130.4367 ymin: 33.00472 xmax: 130.9907 ymax: 33.92599
## Geodetic CRS: JGD2000
## KEY_CODE PREF CITY S_AREA PREF_NAME CITY_NAME S_NAME KIGO_E HCODE
## 1 40 40 202 000000 福岡県 大牟田市 <NA> <NA> 8101
## 2 40 40 202 000000 福岡県 大牟田市 <NA> <NA> 8101
## 3 40 40 202 000000 福岡県 大牟田市 <NA> <NA> 8101
## 4 401010010 40 101 001000 福岡県 門司区 青葉台 <NA> 8101
## 5 401010020 40 101 002000 福岡県 門司区 大字伊川 <NA> 8101
## 6 401010030 40 101 003000 福岡県 門司区 泉ヶ丘 <NA> 8101
## AREA PERIMETER H27KAxx_ H27KAxx_ID KEN KEN_NAME SITYO_NAME GST_NAME
## 1 1729.630 254.400 6725 6724 40 福岡県 <NA> 大牟田市
## 2 11782.225 433.670 6726 6725 40 福岡県 <NA> 大牟田市
## 3 2401.786 316.007 6727 6726 40 福岡県 <NA> 大牟田市
## 4 101577.058 1474.768 6411 6410 40 福岡県 <NA> 北九州市
## 5 4699014.636 11449.527 4686 4685 40 福岡県 <NA> 北九州市
## 6 82882.547 1262.830 6848 6847 40 福岡県 <NA> 北九州市
## CSS_NAME KIHON1 DUMMY1 KIHON2 KEYCODE1 KEYCODE2 AREA_MAX_F KIGO_D N_KEN
## 1 <NA> 0000 - 00 202000000 <NA> <NA> D1 43
## 2 <NA> 0000 - 00 202000000 <NA> <NA> D1 43
## 3 <NA> 0000 - 00 202000000 <NA> <NA> D1 43
## 4 門司区 0010 - 00 101001000 1010010 M <NA> <NA>
## 5 門司区 0020 - 00 101002000 1010020 M <NA> <NA>
## 6 門司区 0030 - 00 101003000 1010030 M <NA> <NA>
## N_CITY KIGO_I MOJI KBSUM JINKO SETAI X_CODE Y_CODE KCODE1
## 1 204 <NA> <NA> 0 0 0 130.4538 33.00534 0000-00
## 2 204 <NA> <NA> 0 0 0 130.4474 33.00582 0000-00
## 3 204 <NA> <NA> 0 0 0 130.4373 33.00599 0000-00
## 4 <NA> <NA> 青葉台 19 256 103 130.9207 33.88516 0010-00
## 5 <NA> <NA> 大字伊川 20 837 373 130.9726 33.91148 0020-00
## 6 <NA> <NA> 泉ヶ丘 14 851 351 130.9351 33.89115 0030-00
## geometry
## 1 POLYGON ((130.4537 33.00497...
## 2 POLYGON ((130.4471 33.00648...
## 3 POLYGON ((130.4367 33.00625...
## 4 POLYGON ((130.9233 33.88419...
## 5 POLYGON ((130.964 33.90387,...
## 6 POLYGON ((130.935 33.88987,...後の分析の際の便利のため，地価データに対して通し番号変数を追加します（鉄道時系列データについては後ほど追加）．また，鉄道時系列データについて，文字列型変数N05_004（供用開始年），N05_005b（設置開始年），N05_005e（設置終了年）を数値型に変換します．
##地価（福岡県）
lp_fukuoka <- lp_fukuoka %>%
dplyr::mutate(LID=dplyr::row_number())
##鉄道時系列（全国・駅データ）
rail_jp <- rail_jp %>%
#変数N05_004からN05_005eを一括で数値型に変換
dplyr::mutate(across(N05_004:N05_005e,~as.numeric(.)))空間データの読み込みには，st_read関数の他に，read_sf関数も使うことができます．read_sf関数ではデフォルトで，引数stringsAsFactorsはFALSE，quietはTRUEになっていますので，コーディングをより節約することができます．
続いて，openxlsxパッケージのread.xlsx関数を使って，テーブルを読み込みます．引数xlsxFileにはxlsxファイルのディレクトリ（Rmdファイルと同じディレクトリに設置されていれば，ファイル名を指定するだけでOK），引数sheetには読み込みたいシート名を指定します．
##店舗データ（福岡市周辺）
fuk_p_lonlat <- openxlsx::read.xlsx(xlsxFile="fuk_p_lonlat.xlsx",
sheet="Sheet 1")
head(fuk_p_lonlat)## osm_id name shop KEN_NAME GST_NAME CSS_NAME
## 1 312503423 マリノアシティ mall 福岡県 福岡市 西区
## 2 473575921 ローソン convenience 福岡県 福岡市 博多区
## 3 477366418 ミニストップ convenience 福岡県 福岡市 東区
## 4 485707416 BOX TOWN yes 福岡県 福岡市 東区
## 5 485707433 ゆめタウン博多 supermarket 福岡県 福岡市 東区
## 6 495981699 サニー supermarket 福岡県 福岡市 博多区
## MOJI lon lat
## 1 小戸２丁目 130.3226 33.59501
## 2 博多駅東２丁目 130.4237 33.58691
## 3 大岳１丁目 130.3446 33.65198
## 4 箱崎４丁目 130.4168 33.61967
## 5 東浜１丁目 130.4120 33.61187
## 6 博多駅南３丁目 130.4264 33.58202投影変換
空間データ読み込みプログラムの実行結果の中で「Geodetic CRS」の部分に注目してください．この項目では，読み込んだ空間データに定義されている投影法が表されます．地価・小地域の空間データの座標系はJGD2000で，鉄道時系列データの座標系はJGD2011で定義されています．座標系は扱う空間データ間で一致している必要があるので，ここでは例えば，鉄道時系列データの座標系をJGD2000に変換します．座標系の変換にはsfのst_transform関数を用います．引数crsには，EPSGコードを元に生成される座標系の情報を指定します．JGD2000を表すEPSGコードは「4612」ですが，コードから座標系の情報を呼び出す際にはsfのst_crs関数を用います．
##鉄道時系列（全国・駅データ）
rail_jp <- rail_jp %>%
#JGD2000に投影変換
sf::st_transform(crs=sf::st_crs(4612))ちなみに，st_crs関数の実行結果な中身は以下のようになっていて，原点や座標系がカバーしているエリアに関する情報を含むことが分かります．
##st_crs関数の実行結果の中身
sf::st_crs(4612)## Coordinate Reference System:
## User input: EPSG:4612
## wkt:
## GEOGCRS["JGD2000",
## DATUM["Japanese Geodetic Datum 2000",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Japan - onshore and offshore."],
## BBOX[17.09,122.38,46.05,157.65]],
## ID["EPSG",4612]]現状，空間データの座標は緯度経度（地理座標系）で定義された状態です．このままでは，様々なジオプロセシングを正確に適用できないので，m単位の座標（投影座標系）に変換を行う必要があります．そこで以下では，各空間データについて，福岡県＆JGD2000に対応する投影座標系である平面直角座標第2系（EPSG:2444）への投影法の変換を行います．
##地価（福岡県）
#座標系を変換
lp_fukuoka_xy <- lp_fukuoka %>%
sf::st_transform(crs=sf::st_crs(2444))
##鉄道時系列（全国）
rail_jp_xy <- rail_jp %>%
sf::st_transform(crs=sf::st_crs(2444))
##小地域（福岡県）
ct_fukuoka_xy <- ct_fukuoka %>%
sf::st_transform(crs=sf::st_crs(2444))実際に，座標値がどのように変わったかを見てみましょう．空間データに含まれる各地物の座標は，地物属性内のgeometryという変数に格納されています．例えば鉄道時系列データの変換前のgeometryは，
##鉄道時系列データ（駅データ）
rail_jp %>%
#geometry変数をdplyrのselect関数で選択
dplyr::select(geometry) %>%
#先頭行を表示する
head()## Simple feature collection with 6 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 141.9359 ymin: 43.62763 xmax: 142.3579 ymax: 43.78933
## Geodetic CRS: JGD2000
## geometry
## 1 POINT (142.3579 43.76344)
## 2 POINT (141.9666 43.69096)
## 3 POINT (142.2717 43.7632)
## 4 POINT (142.0414 43.72125)
## 5 POINT (142.3253 43.78933)
## 6 POINT (141.9359 43.62763)一方，変換後のgeometryは以下のようになり，座標系が変換されたことがわかります．
rail_jp_xy %>%
dplyr::select(geometry) %>%
head()## Simple feature collection with 6 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 882745 ymin: 1238074 xmax: 914730.2 ymax: 1260321
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## geometry
## 1 POINT (914730.2 1257807)
## 2 POINT (884284.8 1245446)
## 3 POINT (907786.9 1256819)
## 4 POINT (889864.8 1249625)
## 5 POINT (911706.7 1260321)
## 6 POINT (882745 1238074)ポリゴンの融合（ディゾルブ）
このセクションと次のセクションを通じて，福岡市域に含まれる地価ポイントのみを抽出する空間データ操作を行います．ここでは，福岡市域に含まれる小地域ポリゴンを抽出した上で，政令区単位でポリゴンに融合します．まずは，小地域データから福岡市域の小地域に対応するポリゴンを抽出します．小地域データの変数のうち，日本語で記された住所情報に該当する変数は以下の通りです．
- KEN_NAME：都道府県名
- GST_NAME：市区町村名
- CSS_NAME：政令区名
- MOJI：小地域名
故に，福岡市域のポリゴンを抽出するには，GST_NAMEが「福岡市」のレコードのみに絞り込めば良いです．また，福岡市域ポリゴンには，「水面」「博多湾」のといった陸地ではない住所や，小地域名がNAのレコードが含まれているので，それらも除外しておきます．
ct_fukuoka_c_xy <- ct_fukuoka_xy %>%
#変数GST_NAME（市区町村名）が「福岡市」に等しいレコードのみ残す
dplyr::filter(GST_NAME=="福岡市") %>%
#小地域名がNAでない陸地ポリゴンのみに絞る
dplyr::filter(!is.na(MOJI)&!(MOJI%in%c("水面","博多湾")))ここで，抽出したポリゴンをプロットしてみましょう．まずは，ggplot2パッケージを用いた方法を示します．空間データの可視化には，geom_sf関数を用います．引数dataには，プロットしたい空間データのオブジェクトを指定します．その他の引数を何も指定しない場合，緯度経度グリッドの上に，灰色で空間データがプロットされます．
#今からプロットしますよという宣言
ggplot2::ggplot() +
#小地域データをプロットする
ggplot2::geom_sf(data=ct_fukuoka_c_xy)
一方，mapviewパッケージのmapview関数を用いると，インタラクティブに地物属性を参照したり，ズームイン・アウトが可能な地図が作成できます．引数xには，プロットしたい空間データのオブジェクトを指定します．その他の引数を何も指定しない場合．WEBから読み込まれた地図（CartoDB.Positron）の上に，水色のポリゴンがプロットされます．基本的な操作方法を説明しておくと，左上のプラス・マイナスボタンで，地図のズームイン・アウトが行えます．その下にある，四角形が何枚か重なった図形で表されるボタンをクリックすると，背景地図や表示させる空間データを選択する画面が現れます．背景地図には，OSMの地図も使うことができます．なお，mapview関数で作成した地図の背景は，インターネットに接続しない状態では表示できないので，注意が必要です．
#小地域データをプロットする
mapview::mapview(x=ct_fukuoka_c_xy)次に，抽出された小地域ポリゴンを，政令区単位で融合します．少々回りくどい方法ですが，ポリゴンの融合は，dplyrパッケージのgroup_by関数とsummarise関数を組み合わせることで，以下のように実行できます．
ct_fukuoka_c_xy_ward <- ct_fukuoka_c_xy %>%
#変数CSS_NAME（政令区名）でポリゴンをグループ化
dplyr::group_by(CSS_NAME) %>%
#グループに基づく集計
dplyr::summarise() %>%
#グループ情報を削除
dplyr::ungroup()
head(ct_fukuoka_c_xy_ward)## Simple feature collection with 6 features and 1 field
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -67669.94 ymin: 47375.32 xmax: -46845.54 ymax: 79180.44
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## # A tibble: 6 × 2
## CSS_NAME geometry
## <chr> <GEOMETRY [m]>
## 1 中央区 POLYGON ((-57942.63 63269.97, -57942.29 63272.12, -57932.58 63335.14…
## 2 南区 POLYGON ((-55152.3 57952.55, -55159.03 57947.9, -55163.97 57931.66, …
## 3 博多区 POLYGON ((-54320.12 64164.49, -54377.79 64217.35, -54404.17 64263.34…
## 4 城南区 POLYGON ((-59145.47 59293.05, -59164.53 59289.48, -59179.74 59284.54…
## 5 早良区 POLYGON ((-62755.97 50693.41, -62820.23 50716.51, -62887.4 50750.17,…
## 6 東区 MULTIPOLYGON (((-53975.99 67588.37, -54012.41 67608.75, -54017.23 67…融合したポリゴンをプロットしてみます．まずは，geom_sf関数での結果です．
ggplot2::ggplot() +
ggplot2::geom_sf(data=ct_fukuoka_c_xy_ward)
続いて，mapview関数での結果です．mapview関数は，自動的に凡例を作成する機能が装備されていて，今回は政令区ごとの色分けが表示されています．引数xに，プロットしたい空間データを指定します．
mapview::mapview(x=ct_fukuoka_c_xy_ward)空間データ間の共通部分を取得（交差，クリップ）
上で融合した福岡市域ポリゴンを用いて，福岡県全体の地価ポイントから，福岡市域に含まれるポイントのみを抽出してみます．この抽出作業を言い換えるならば，福岡市域ポリゴンと，福岡県全体の地価ポイントとの共通部分を取る作業になります．
このジオプロセシングは，sfのst_intersection関数を用いて行います．引数x,yで，共通部分を取る空間データを指定します．
まず，福岡市域ポリゴンの属性を地価ポイントに付与する場合のプログラムは以下のようになります．福岡市域ポリゴン側の変数CSS_NAME（政令区名）が付与されていることが確認できます．レコード数は933から293に減りました．
#共通部分を取る前のデータのレコード数
dim(lp_fukuoka_xy)## [1] 933 141#地価ポイントと福岡市域ポリゴンの共通部分を取る
lp_fukuoka_xy_c <- sf::st_intersection(x=lp_fukuoka_xy,
y=ct_fukuoka_c_xy_ward)
head(lp_fukuoka_xy_c)## Simple feature collection with 6 features and 141 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58150.03 ymin: 63836.13 xmax: -56564 ymax: 66100.34
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## 336 000 001 000 001 2022 575000 8.9 1 false
## 337 000 002 000 002 2022 870000 5.5 1 false
## 338 000 003 000 003 2022 324000 12.1 1 false
## 339 000 004 000 004 2022 488000 6.6 1 false
## 340 000 005 000 005 2022 461000 7.7 1 false
## 341 000 006 000 006 2022 300000 12.8 1 false
## L01_010 L01_011 L01_012 L01_013 L01_014 L01_015 L01_016 L01_017 L01_018
## 336 false false false false false false false false false
## 337 false false false false false false false false false
## 338 false false false false false false false false false
## 339 false false false false false false false false false
## 340 false false false false false false false false false
## 341 false false false false false false false false false
## L01_019 L01_020 L01_021 L01_022 L01_023
## 336 false false false 40133 福岡中央
## 337 false false false 40133 福岡中央
## 338 false false false 40133 福岡中央
## 339 false false false 40133 福岡中央
## 340 false false false 40133 福岡中央
## 341 false false false 40133 福岡中央
## L01_024 L01_025 L01_026
## 336 福岡県 福岡市中央区赤坂２丁目２区５４番 赤坂２−２−１１ 737
## 337 福岡県 福岡市中央区大濠１丁目１３区１３３番 大濠１−１３−２６ 718
## 338 福岡県 福岡市中央区警固３丁目１０３番 _ 225
## 339 福岡県 福岡市中央区荒戸３丁目２９７番１ 荒戸３−５−５１ 1085
## 340 福岡県 福岡市中央区六本松４丁目６１番 六本松４−５−１８ 330
## 341 福岡県 福岡市中央区谷１丁目２９３番 谷１−９−２１ 161
## L01_027 L01_028 L01_029 L01_030 L01_031 L01_032 L01_033 L01_034
## 336 住宅,その他 共同住宅 001 RC true true true _
## 337 住宅 _ 001 W true true true _
## 338 空地 _ 004 その他 true true true 台形
## 339 住宅,その他 共同住宅 001 RC true true true _
## 340 住宅,その他 共同住宅 001 RC true true true _
## 341 住宅 _ 001 W true true true _
## L01_035 L01_036 L01_037 L01_038 L01_039 L01_040 L01_041 L01_042 L01_043
## 336 1 2.0 6 0 市道 南 5.9 _ _
## 337 1 1.0 2 0 市道 東 5.4 _ _
## 338 1 1.5 0 0 市道 北 4.0 _ _
## 339 1 2.0 7 0 市道 北 6.8 _ _
## 340 1 1.0 3 0 市道 南 7.0 _ _
## 341 1 2.0 2 0 市道 北東 5.6 _ _
## L01_044 L01_045 L01_046 L01_047
## 336 _ _ _ マンション等が建ち並ぶ都心に近い住宅地域
## 337 側道 南 _ 大規模一般住宅、マンションが建ち並ぶ住宅地域
## 338 _ _ _ 一般住宅の中に共同住宅が見られる住宅地域
## 339 _ _ _ 中規模住宅、マンション等が混在する住宅地域
## 340 _ _ _ 共同住宅の中に一般住宅も見られる住宅地域
## 341 _ _ _ 中規模住宅、アパート等が多い既成住宅地域
## L01_048 L01_049 L01_050 L01_051 L01_052 L01_053 L01_054 L01_055 L01_056
## 336 赤坂 560 2住居 準防 市街化 0 _ _ 60
## 337 唐人町 750 1住居 準防 市街化 0 _ _ 60
## 338 桜坂 420 1住居 準防 市街化 0 _ _ 60
## 339 大濠公園 800 1住居 準防 市街化 0 _ _ 60
## 340 六本松 500 1住居 _ 市街化 0 _ _ 60
## 341 六本松 430 1住居 _ 市街化 0 _ _ 60
## L01_057 L01_058 L01_059 L01_060 L01_061
## 336 200 false false 0000000000000000001111111111111111111111 0
## 337 200 false false 1111111111111111111111111111111111111111 252000
## 338 200 false false 0000000000000000000000000000000111111111 0
## 339 200 false false 0000000111111111111111111111111111111111 0
## 340 200 false true 0000001111111111111111111111111111111111 0
## 341 200 false false 0000000000000000000000000000001111111111 0
## L01_062 L01_063 L01_064 L01_065 L01_066 L01_067 L01_068 L01_069 L01_070
## 336 0 0 0 0 0 0 0 0 0
## 337 267000 273000 277000 305000 460000 550000 650000 750000 750000
## 338 0 0 0 0 0 0 0 0 0
## 339 0 0 0 0 0 0 450000 519000 519000
## 340 0 0 0 0 0 270000 300000 365000 365000
## 341 0 0 0 0 0 0 0 0 0
## L01_071 L01_072 L01_073 L01_074 L01_075 L01_076 L01_077 L01_078 L01_079
## 336 0 0 0 0 0 0 0 0 260000
## 337 650000 600000 555000 520000 475000 443000 423000 410000 400000
## 338 0 0 0 0 0 0 0 0 0
## 339 455000 385000 340000 308000 285000 270000 260000 248000 240000
## 340 328000 315000 275000 245000 227000 216000 212000 207000 203000
## 341 0 0 0 0 0 0 0 0 0
## L01_080 L01_081 L01_082 L01_083 L01_084 L01_085 L01_086 L01_087 L01_088
## 336 257000 253000 251000 252000 259000 300000 331000 308000 293000
## 337 393000 393000 393000 400000 420000 525000 603000 535000 494000
## 338 0 0 0 0 0 0 0 0 0
## 339 228000 220000 214000 214000 221000 284000 328000 310000 290000
## 340 195000 187000 181000 181000 182000 214000 235000 223000 215000
## 341 0 0 0 0 0 0 0 0 0
## L01_089 L01_090 L01_091 L01_092 L01_093 L01_094 L01_095 L01_096 L01_097
## 336 291000 292000 305000 321000 340000 359000 391000 425000 463000
## 337 494000 501000 518000 534000 566000 605000 650000 698000 756000
## 338 0 0 0 178000 186000 197000 208000 219000 232000
## 339 281000 282000 289000 299000 311000 329000 348000 373000 389000
## 340 208000 213000 217000 229000 242000 262000 288000 316000 346000
## 341 0 0 159000 164000 173000 187000 201000 212000 227000
## L01_098 L01_099 L01_100 L01_101 L01_102 L01_103
## 336 504000 528000 575000 00000000000000 00000000000000 00000000000000
## 337 800000 825000 870000 10000000000000 10000000000000 10000000000000
## 338 266000 289000 324000 00000000000000 00000000000000 00000000000000
## 339 438000 458000 488000 00000000000000 00000000000000 00000000000000
## 340 394000 428000 461000 00000000000000 00000000000000 00000000000000
## 341 249000 266000 300000 00000000000000 00000000000000 00000000000000
## L01_104 L01_105 L01_106 L01_107 L01_108
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 00000000000000 00000000000000 00000000000000 40000000000000 10010000000000
## 340 00000000000000 00000000000000 40000000000000 10000000000000 10010000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_109 L01_110 L01_111 L01_112 L01_113
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_114 L01_115 L01_116 L01_117 L01_118
## 336 00000000000000 00000000000000 00000000000000 00000000000000 40000000000000
## 337 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_119 L01_120 L01_121 L01_122 L01_123
## 336 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000010 10000000000000 10000000000000 10000010000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_124 L01_125 L01_126 L01_127 L01_128
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_129 L01_130 L01_131 L01_132 L01_133
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 338 00000000000000 00000000000000 40000000000000 10000000000000 10000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 40000000000000 10000000000000 10000000000000 10000000000000
## L01_134 L01_135 L01_136 L01_137 L01_138
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 10000000000000 10000000000000 10000010000000 10000000000000 11011000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 10000000000000 10000000000000 10000010000000 10000000000000 10000000000000
## L01_139 LID CSS_NAME geometry
## 336 10000000000000 336 中央区 POINT (-56754.69 65075.79)
## 337 10000000000000 337 中央区 POINT (-58150.03 65083.26)
## 338 10000000000000 338 中央区 POINT (-56564 64307.29)
## 339 10000000000000 339 中央区 POINT (-57971.52 66100.34)
## 340 10000000000000 340 中央区 POINT (-57697.73 63836.13)
## 341 10000000000000 341 中央区 POINT (-57460.85 64153.45)#共通部分を取った後のデータのレコード数
dim(lp_fukuoka_xy_c)## [1] 293 142次に，福岡市域ポリゴンの属性を地価ポイントに付与せずに共通部分を取る場合は，以下のようになります．属性付与を行う場合との違いは，引数yに市域ポリゴンそのものではなく，市域ポリゴンのgeometry属性を指定している点です．
#地価ポイントと福岡市域ポリゴンの共通部分を取る
lp_fukuoka_xy_c_geom <- sf::st_intersection(x=lp_fukuoka_xy,
y=sf::st_geometry(ct_fukuoka_c_xy_ward))
head(lp_fukuoka_xy_c_geom)## Simple feature collection with 6 features and 140 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58150.03 ymin: 63836.13 xmax: -56564 ymax: 66100.34
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## 336 000 001 000 001 2022 575000 8.9 1 false
## 337 000 002 000 002 2022 870000 5.5 1 false
## 338 000 003 000 003 2022 324000 12.1 1 false
## 339 000 004 000 004 2022 488000 6.6 1 false
## 340 000 005 000 005 2022 461000 7.7 1 false
## 341 000 006 000 006 2022 300000 12.8 1 false
## L01_010 L01_011 L01_012 L01_013 L01_014 L01_015 L01_016 L01_017 L01_018
## 336 false false false false false false false false false
## 337 false false false false false false false false false
## 338 false false false false false false false false false
## 339 false false false false false false false false false
## 340 false false false false false false false false false
## 341 false false false false false false false false false
## L01_019 L01_020 L01_021 L01_022 L01_023
## 336 false false false 40133 福岡中央
## 337 false false false 40133 福岡中央
## 338 false false false 40133 福岡中央
## 339 false false false 40133 福岡中央
## 340 false false false 40133 福岡中央
## 341 false false false 40133 福岡中央
## L01_024 L01_025 L01_026
## 336 福岡県 福岡市中央区赤坂２丁目２区５４番 赤坂２−２−１１ 737
## 337 福岡県 福岡市中央区大濠１丁目１３区１３３番 大濠１−１３−２６ 718
## 338 福岡県 福岡市中央区警固３丁目１０３番 _ 225
## 339 福岡県 福岡市中央区荒戸３丁目２９７番１ 荒戸３−５−５１ 1085
## 340 福岡県 福岡市中央区六本松４丁目６１番 六本松４−５−１８ 330
## 341 福岡県 福岡市中央区谷１丁目２９３番 谷１−９−２１ 161
## L01_027 L01_028 L01_029 L01_030 L01_031 L01_032 L01_033 L01_034
## 336 住宅,その他 共同住宅 001 RC true true true _
## 337 住宅 _ 001 W true true true _
## 338 空地 _ 004 その他 true true true 台形
## 339 住宅,その他 共同住宅 001 RC true true true _
## 340 住宅,その他 共同住宅 001 RC true true true _
## 341 住宅 _ 001 W true true true _
## L01_035 L01_036 L01_037 L01_038 L01_039 L01_040 L01_041 L01_042 L01_043
## 336 1 2.0 6 0 市道 南 5.9 _ _
## 337 1 1.0 2 0 市道 東 5.4 _ _
## 338 1 1.5 0 0 市道 北 4.0 _ _
## 339 1 2.0 7 0 市道 北 6.8 _ _
## 340 1 1.0 3 0 市道 南 7.0 _ _
## 341 1 2.0 2 0 市道 北東 5.6 _ _
## L01_044 L01_045 L01_046 L01_047
## 336 _ _ _ マンション等が建ち並ぶ都心に近い住宅地域
## 337 側道 南 _ 大規模一般住宅、マンションが建ち並ぶ住宅地域
## 338 _ _ _ 一般住宅の中に共同住宅が見られる住宅地域
## 339 _ _ _ 中規模住宅、マンション等が混在する住宅地域
## 340 _ _ _ 共同住宅の中に一般住宅も見られる住宅地域
## 341 _ _ _ 中規模住宅、アパート等が多い既成住宅地域
## L01_048 L01_049 L01_050 L01_051 L01_052 L01_053 L01_054 L01_055 L01_056
## 336 赤坂 560 2住居 準防 市街化 0 _ _ 60
## 337 唐人町 750 1住居 準防 市街化 0 _ _ 60
## 338 桜坂 420 1住居 準防 市街化 0 _ _ 60
## 339 大濠公園 800 1住居 準防 市街化 0 _ _ 60
## 340 六本松 500 1住居 _ 市街化 0 _ _ 60
## 341 六本松 430 1住居 _ 市街化 0 _ _ 60
## L01_057 L01_058 L01_059 L01_060 L01_061
## 336 200 false false 0000000000000000001111111111111111111111 0
## 337 200 false false 1111111111111111111111111111111111111111 252000
## 338 200 false false 0000000000000000000000000000000111111111 0
## 339 200 false false 0000000111111111111111111111111111111111 0
## 340 200 false true 0000001111111111111111111111111111111111 0
## 341 200 false false 0000000000000000000000000000001111111111 0
## L01_062 L01_063 L01_064 L01_065 L01_066 L01_067 L01_068 L01_069 L01_070
## 336 0 0 0 0 0 0 0 0 0
## 337 267000 273000 277000 305000 460000 550000 650000 750000 750000
## 338 0 0 0 0 0 0 0 0 0
## 339 0 0 0 0 0 0 450000 519000 519000
## 340 0 0 0 0 0 270000 300000 365000 365000
## 341 0 0 0 0 0 0 0 0 0
## L01_071 L01_072 L01_073 L01_074 L01_075 L01_076 L01_077 L01_078 L01_079
## 336 0 0 0 0 0 0 0 0 260000
## 337 650000 600000 555000 520000 475000 443000 423000 410000 400000
## 338 0 0 0 0 0 0 0 0 0
## 339 455000 385000 340000 308000 285000 270000 260000 248000 240000
## 340 328000 315000 275000 245000 227000 216000 212000 207000 203000
## 341 0 0 0 0 0 0 0 0 0
## L01_080 L01_081 L01_082 L01_083 L01_084 L01_085 L01_086 L01_087 L01_088
## 336 257000 253000 251000 252000 259000 300000 331000 308000 293000
## 337 393000 393000 393000 400000 420000 525000 603000 535000 494000
## 338 0 0 0 0 0 0 0 0 0
## 339 228000 220000 214000 214000 221000 284000 328000 310000 290000
## 340 195000 187000 181000 181000 182000 214000 235000 223000 215000
## 341 0 0 0 0 0 0 0 0 0
## L01_089 L01_090 L01_091 L01_092 L01_093 L01_094 L01_095 L01_096 L01_097
## 336 291000 292000 305000 321000 340000 359000 391000 425000 463000
## 337 494000 501000 518000 534000 566000 605000 650000 698000 756000
## 338 0 0 0 178000 186000 197000 208000 219000 232000
## 339 281000 282000 289000 299000 311000 329000 348000 373000 389000
## 340 208000 213000 217000 229000 242000 262000 288000 316000 346000
## 341 0 0 159000 164000 173000 187000 201000 212000 227000
## L01_098 L01_099 L01_100 L01_101 L01_102 L01_103
## 336 504000 528000 575000 00000000000000 00000000000000 00000000000000
## 337 800000 825000 870000 10000000000000 10000000000000 10000000000000
## 338 266000 289000 324000 00000000000000 00000000000000 00000000000000
## 339 438000 458000 488000 00000000000000 00000000000000 00000000000000
## 340 394000 428000 461000 00000000000000 00000000000000 00000000000000
## 341 249000 266000 300000 00000000000000 00000000000000 00000000000000
## L01_104 L01_105 L01_106 L01_107 L01_108
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 00000000000000 00000000000000 00000000000000 40000000000000 10010000000000
## 340 00000000000000 00000000000000 40000000000000 10000000000000 10010000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_109 L01_110 L01_111 L01_112 L01_113
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_114 L01_115 L01_116 L01_117 L01_118
## 336 00000000000000 00000000000000 00000000000000 00000000000000 40000000000000
## 337 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_119 L01_120 L01_121 L01_122 L01_123
## 336 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000010 10000000000000 10000000000000 10000010000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_124 L01_125 L01_126 L01_127 L01_128
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_129 L01_130 L01_131 L01_132 L01_133
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 338 00000000000000 00000000000000 40000000000000 10000000000000 10000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 40000000000000 10000000000000 10000000000000 10000000000000
## L01_134 L01_135 L01_136 L01_137 L01_138
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 10000000000000 10000000000000 10000010000000 10000000000000 11011000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 10000000000000 10000000000000 10000010000000 10000000000000 10000000000000
## L01_139 LID geometry
## 336 10000000000000 336 POINT (-56754.69 65075.79)
## 337 10000000000000 337 POINT (-58150.03 65083.26)
## 338 10000000000000 338 POINT (-56564 64307.29)
## 339 10000000000000 339 POINT (-57971.52 66100.34)
## 340 10000000000000 340 POINT (-57697.73 63836.13)
## 341 10000000000000 341 POINT (-57460.85 64153.45)sfオブジェクトのgeometry属性は，冒頭で説明したsfcというクラスに該当します．すなわち，地物属性なし・座標系情報ありの空間データです．
空間データから一定の距離圏にある別の空間データを抽出（バッファ）
バッファと共通部分抽出を組み合わせた処理
この空間データ操作は，ある地物から任意の距離圏内を覆うバッファというポリゴンを生成するものです．それ単独で用いるというよりは，上で紹介した共通部分を取る操作と組み合わせることで，そのバッファ内に含まれる別の地物を抽出する，という使い方をすることが多いです．ここでは例として，福岡市営地下鉄の各駅から500m圏内にある地価ポイントを抽出するという作業を実装してみます．
まずは，鉄道時系列データから，福岡市営地下鉄の駅だけを抽出してみます．各駅の運営事業者に関する変数N05_003が「福岡市」に等しく，かつ，設置終了時点の変数N05_005eが9999年（即ち現在も存続していることを表す）となっている駅ポイントのみを抽出します．
rail_jp_xy_fukuoka_sub <- rail_jp_xy %>%
dplyr::filter(N05_003=="福岡市"&N05_005e=="9999")次に，st_buffer関数を用いて各駅ポイントから500mのバッファを発生させます．引数xではバッファの発生元になる空間データ（今回は地下鉄駅ポイント），引数distではバッファが覆う距離圏を指定します．distの単位は平面直角座標系の単位（m）です．
rail_jp_xy_fukuoka_sub_buffer <- sf::st_buffer(x=rail_jp_xy_fukuoka_sub,dist=500)続いて，この500mバッファを用いて，地下鉄駅から500m圏内にある地価ポイントを抽出します．先ほどと同樣，st_intersection関数を用います．
lp_fukuoka_xy_c_sub_buffer_500m <- sf::st_intersection(x=lp_fukuoka_xy_c,
y=rail_jp_xy_fukuoka_sub_buffer)
head(lp_fukuoka_xy_c_sub_buffer_500m)## Simple feature collection with 6 features and 152 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -62892.13 ymin: 64452.54 xmax: -53032.13 ymax: 65884.26
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## 429 000 002 000 002 2022 223000 8.3 1 false
## 443 000 016 000 016 2022 258000 7.1 1 false
## 452 005 001 005 001 2022 525000 10.5 1 false
## 502 000 027 000 027 2022 492000 3.8 1 false
## 521 005 002 005 002 2022 688000 11.9 1 false
## 327 005 019 005 019 2022 352000 13.9 1 false
## L01_010 L01_011 L01_012 L01_013 L01_014 L01_015 L01_016 L01_017 L01_018
## 429 false false false false false false false false false
## 443 false false false false false false false false false
## 452 false false false false false false false false false
## 502 false false false false false false false false false
## 521 false false false false false false false false false
## 327 false false false false false false false false false
## L01_019 L01_020 L01_021 L01_022 L01_023
## 429 false false false 40135 福岡西
## 443 false false false 40135 福岡西
## 452 false false false 40135 福岡西
## 502 false false false 40137 福岡早良
## 521 false false false 40137 福岡早良
## 327 false false false 40132 福岡博多
## L01_024 L01_025 L01_026
## 429 福岡県 福岡市西区姪の浜５丁目９０６番１外 姪の浜５−４−１５ 1020
## 443 福岡県 福岡市西区姪浜駅南２丁目１１７番 姪浜駅南２−１６−１７ 132
## 452 福岡県 福岡市西区姪浜駅南１丁目７１番 姪浜駅南１−６−１９ 372
## 502 福岡県 福岡市早良区高取２丁目２２１番２外 高取２−５−１０ 1907
## 521 福岡県 福岡市早良区高取１丁目２９番３外 高取１−１−４２ 511
## 327 福岡県 福岡市博多区東比恵２丁目１２番 東比恵２−２−１０ 414
## L01_027 L01_028 L01_029 L01_030 L01_031 L01_032 L01_033
## 429 住宅,その他 共同住宅 001 RC true true true
## 443 住宅 _ 001 W true true true
## 452 住宅,店舗,事務所,その他 共同住宅 001 RC true true true
## 502 住宅,その他 共同住宅 001 RC true true true
## 521 事務所 _ 001 SRC true true true
## 327 住宅,事務所,その他 共同住宅 001 RC true true true
## L01_034 L01_035 L01_036 L01_037 L01_038 L01_039 L01_040 L01_041 L01_042
## 429 _ 1.0 1.5 9 0 市道 南 5 _
## 443 _ 1.5 1.0 3 0 市道 西 9 _
## 452 _ 1.0 1.0 7 0 県道 西 25 _
## 502 _ 1.5 1.0 8 0 市道 東 6 _
## 521 _ 1.0 1.2 6 0 市道 北 25 _
## 327 _ 1.0 2.5 6 0 市道 北西 11 _
## L01_043 L01_044 L01_045 L01_046
## 429 _ _ _ _
## 443 _ _ _ _
## 452 _ _ _ _
## 502 _ _ _ _
## 521 _ _ _ _
## 327 _ _ _ _
## L01_047 L01_048 L01_049 L01_050
## 429 共同住宅が多く見られる利便性のよい住宅地域 姪浜 450 1住居
## 443 一般住宅等が建ち並ぶ区画整然とした住宅地域 姪浜 500 1住居
## 452 中高層の店舗兼共同住宅等が建ち並ぶ商業地域 姪浜 170 商業
## 502 一般住宅、マンション等が建ち並ぶ既成住宅地域 藤崎 410 1住居
## 521 店舗、事務所、共同住宅が混在する商業地域 藤崎 250 商業
## 327 事務所、工場、共同住宅等が混在する商業地域 東比恵 550 準工
## L01_051 L01_052 L01_053 L01_054 L01_055 L01_056 L01_057 L01_058 L01_059
## 429 _ 市街化 0 _ _ 60 200 false false
## 443 _ 市街化 0 _ _ 60 200 false false
## 452 準防 市街化 0 _ _ 80 400 false false
## 502 _ 市街化 0 _ _ 60 200 false false
## 521 準防 市街化 0 _ _ 80 400 false false
## 327 _ 市街化 0 _ _ 60 200 false false
## L01_060 L01_061 L01_062 L01_063 L01_064
## 429 0000000000000000000000000000001111111111 0 0 0 0
## 443 0000000000000000000000111111111111111111 0 0 0 0
## 452 0000000000000000000001111111111111111111 0 0 0 0
## 502 0000000000000011111111111111111111111111 0 0 0 0
## 521 0000001111111111111111111111111111111111 0 0 0 0
## 327 0000000111111111111111111111111111111111 0 0 0 0
## L01_065 L01_066 L01_067 L01_068 L01_069 L01_070 L01_071 L01_072 L01_073
## 429 0 0 0 0 0 0 0 0 0
## 443 0 0 0 0 0 0 0 0 0
## 452 0 0 0 0 0 0 0 0 0
## 502 0 0 0 0 0 0 0 0 0
## 521 0 0 1100000 1210000 1350000 1350000 1200000 1030000 860000
## 327 0 0 0 500000 570000 570000 520000 465000 425000
## L01_074 L01_075 L01_076 L01_077 L01_078 L01_079 L01_080 L01_081 L01_082
## 429 0 0 0 0 0 0 0 0 0
## 443 0 0 0 0 0 0 0 0 0
## 452 0 0 0 0 0 0 0 0 340000
## 502 0 280000 260000 249000 240000 233000 225000 219000 212000
## 521 700000 590000 530000 475000 446000 426000 410000 393000 376000
## 327 370000 315000 290000 273000 256000 242000 230000 217000 202000
## L01_083 L01_084 L01_085 L01_086 L01_087 L01_088 L01_089 L01_090 L01_091
## 429 0 0 0 0 0 0 0 0 136000
## 443 170000 168000 170000 174000 173000 168000 166000 167000 171000
## 452 338000 338000 360000 380000 357000 335000 318000 316000 316000
## 502 208000 211000 227000 250000 239000 229000 229000 234000 259000
## 521 358000 345000 348000 359000 341000 315000 315000 318000 335000
## 327 189000 181000 182000 195000 185000 173000 164000 158000 158000
## L01_092 L01_093 L01_094 L01_095 L01_096 L01_097 L01_098 L01_099 L01_100
## 429 142000 147000 152000 160000 171000 183000 200000 206000 223000
## 443 178000 185000 191000 198000 206000 215000 235000 241000 258000
## 452 317000 318000 328000 345000 372000 400000 450000 475000 525000
## 502 274000 305000 340000 372000 407000 449000 471000 474000 492000
## 521 354000 385000 413000 441000 480000 526000 579000 615000 688000
## 327 165000 174000 186000 200000 222000 246000 281000 309000 352000
## L01_101 L01_102 L01_103 L01_104 L01_105
## 429 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 443 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 452 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 502 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 521 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 327 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_106 L01_107 L01_108 L01_109 L01_110
## 429 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 443 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 452 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 502 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 521 40000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 327 00000000000000 40000000000000 10010000000000 10010000000000 10000000000000
## L01_111 L01_112 L01_113 L01_114 L01_115
## 429 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 443 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 452 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 502 00000000000000 00000000000000 00000000000000 40000000000000 10000000000000
## 521 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 327 10000010000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_116 L01_117 L01_118 L01_119 L01_120
## 429 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 443 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 452 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 502 10000000000000 10000000000000 10000000000000 10000000000000 10000000000010
## 521 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 327 10000000000000 10000000000000 10000000000000 10000000000000 10000000000010
## L01_121 L01_122 L01_123 L01_124 L01_125
## 429 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 443 00000000000000 40000000000000 10000000000000 10000000000000 10000000000000
## 452 40000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 502 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 521 10000000000010 10000000000000 10000000000000 10000000000000 10000000000000
## 327 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_126 L01_127 L01_128 L01_129 L01_130
## 429 00000000000000 00000000000000 00000000000000 00000000000000 40000000000000
## 443 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 452 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 502 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 521 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 327 10000000000000 10000000000000 10000000000000 10000000000000 20000000000000
## L01_131 L01_132 L01_133 L01_134 L01_135
## 429 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 443 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 452 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 502 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 521 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 327 10000000000000 10000010000000 10000000000000 10000000000000 10000000000000
## L01_136 L01_137 L01_138 L01_139 LID CSS_NAME
## 429 10000000000000 10000000000000 10000000000000 10000000000000 429 西区
## 443 10000000000000 10000000000000 10000000000000 10000000000000 443 西区
## 452 10000000000000 10000000000000 10000000000000 10000000000000 452 西区
## 502 10000000000000 10000000000000 10000000000000 10000000000000 502 早良区
## 521 10000000000000 10000000000000 10000000000000 10000000000000 521 早良区
## 327 10000000000000 10000000000000 10000000000000 10000000000000 327 博多区
## N05_001 N05_002 N05_003 N05_004 N05_005b N05_005e N05_006 N05_007
## 429 3 1号線 福岡市 1983 1983 9999 EB03_22318001 <NA>
## 443 3 1号線 福岡市 1983 1983 9999 EB03_22318001 <NA>
## 452 3 1号線 福岡市 1983 1983 9999 EB03_22318001 <NA>
## 502 3 1号線 福岡市 1981 1981 9999 EB03_22318003 <NA>
## 521 3 1号線 福岡市 1981 1981 9999 EB03_22318003 <NA>
## 327 3 1号線 福岡市 1993 1993 9999 EB03_22318013 <NA>
## N05_008 N05_009 N05_011 geometry
## 429 <NA> <NA> 姪浜 POINT (-62892.13 65186.41)
## 443 <NA> <NA> 姪浜 POINT (-62733.53 64497.85)
## 452 <NA> <NA> 姪浜 POINT (-62764.13 64831.12)
## 502 <NA> <NA> 藤崎 POINT (-60278.71 64452.54)
## 521 <NA> <NA> 藤崎 POINT (-60153.28 64720.37)
## 327 <NA> <NA> 東比恵 POINT (-53032.13 65884.26)処理結果の集計・プロット（geom_sf関数使用）
福岡市域ポリゴンのプロットに，500mバッファと抽出された地価ポイントを重ね描きします．まずは，geom_sf関数を用いた結果を示します．今回は複数の地物のプロットを行いますので，地図を見やすくするために以下の工夫を行います．まず，各ポリゴン・ポイントの色分けを行います．色分けを指定する引数は，ポリゴンの場合はfill，ポイント・ラインの場合はcolorとなります．
ggplot2::ggplot() +
#福岡市域ポリゴンを白でプロット
ggplot2::geom_sf(data=ct_fukuoka_c_xy_ward,fill="white") +
#500mバッファを緑でプロット
ggplot2::geom_sf(data=rail_jp_xy_fukuoka_sub_buffer,fill="green") +
#地価ポイントを赤でプロット
ggplot2::geom_sf(data=lp_fukuoka_xy_c_sub_buffer_500m,color="red")
次に，プロット描画領域を，500mバッファを含むものだけに絞り込みます．空間データの空間範囲の取得には，sfのst_bbox関数を用います．例えば，500mバッファの空間範囲は以下のようになります．x（東西）,y（南北）軸の各々について，座標の最大値・最小値が格納されたベクトルが取得できます．
sf::st_bbox(obj=rail_jp_xy_fukuoka_sub_buffer)## xmin ymin xmax ymax
## -63551.92 60197.50 -50733.25 70720.93st_bbox関数を用いてプロットの空間範囲を限定する際には，ggplot2のcoord_sf関数を用います．引数xlimには最大値・最小値ベクトルの要素のうちx軸方向に関するもの，ylimにはy軸方向に関するものを指定します．
ggplot2::ggplot() +
#福岡市域ポリゴンを白でプロット
ggplot2::geom_sf(data=ct_fukuoka_c_xy_ward,fill="white") +
#500mバッファを緑でプロット
ggplot2::geom_sf(data=rail_jp_xy_fukuoka_sub_buffer,fill="green") +
#地価ポイントを赤でプロット
ggplot2::geom_sf(data=lp_fukuoka_xy_c_sub_buffer_500m,color="red") +
#xlimにx座標方向の最大・最小値を指定
ggplot2::coord_sf(xlim=sf::st_bbox(obj=rail_jp_xy_fukuoka_sub_buffer)[c("xmin","xmax")],
#ylimにy座標方向の最大・最小値を指定
ylim=sf::st_bbox(obj=rail_jp_xy_fukuoka_sub_buffer)[c("ymin","ymax")])
更なる応用として，バッファ毎に平均対数地価を計算し，その値に応じたプロットを行います．まず，dplyrの各変数を組み合わせることで，各バッファについて平均対数地価を計算します．その際geometry属性は不要なので，集計作業に先立ち，sfのst_set_geometry関数を用いて属性を削除します．
#各バッファ内の平均地価を計算
lp_fukuoka_xy_c_sub_buffer_500m_avg_tbl <- lp_fukuoka_xy_c_sub_buffer_500m %>%
#geometry属性を削除する
sf::st_set_geometry(value=NULL) %>%
#地価が入った変数L01_006が文字列変数になってしまっているので，数値型に変換
dplyr::mutate(lp=as.numeric(L01_006)) %>%
#変数lpの対数値を変数ln_lpとして追加
dplyr::mutate(ln_lp=log(lp)) %>%
#鉄道時系列データから結合された変数N05_011（駅名）でグループ化
dplyr::group_by(N05_011) %>%
#グループ毎に平均対数地価を計算
dplyr::summarise(ln_lp_avg=mean(ln_lp))次に，計算された平均対数地価を，500mバッファに左結合します．その際，地価観測地点が含まれないために平均対数地価がNAとなったバッファは削除します．
rail_jp_xy_fukuoka_sub_buffer_lp_avg <- rail_jp_xy_fukuoka_sub_buffer %>%
#変数N05_011（駅名）をキー変数として平均対数地価を結合
dplyr::left_join(y=lp_fukuoka_xy_c_sub_buffer_500m_avg_tbl,by="N05_011") %>%
#平均対数地価が結合されたバッファのみ残す
dplyr::filter(!is.na(ln_lp_avg))最後に，平均対数地価で色分けしたバッファのプロットを作成します．ある変数に基づいて空間データを色分けしたい場合には，geom_sf関数の引数fill（ポイント・ラインの場合はcolor）をaes()という関数で囲った上で，後続のプログラムで色分けの詳細（パレットや色分けの分級方法等）を指定します．色分けの指定をしているのがscale_fill_distiller関数で，これはポリゴンを何かしらの変数・パレットに基づいて色分けする方法を指定する関数です．パレットにはSpectral（青系→黄色→赤系とグラデーションが変化する）を用いることと，凡例の名前を「Log Land Price」とすることを指定しています．加えて，ggspatialのannotation_scale関数を用いて縮尺を，annotation_north_arrow関数を用いて方位を地図に追加します．
ggplot2::ggplot() +
#福岡市域ポリゴンを白でプロット
ggplot2::geom_sf(data=ct_fukuoka_c_xy_ward,fill="white") +
#「500mバッファを変数ln_lp_avgに基づいて色分けしプロットする」という指定
ggplot2::geom_sf(data=rail_jp_xy_fukuoka_sub_buffer_lp_avg,aes(fill=ln_lp_avg)) +
#Spectralというパレットを用いて色分け（同時に，凡例の名前を「Log Land Price」と指定）
ggplot2::scale_fill_distiller(palette="Spectral",name="Log Land Price") +
#xlimにx座標方向の最大・最小値を指定
ggplot2::coord_sf(xlim=sf::st_bbox(obj=rail_jp_xy_fukuoka_sub_buffer)[c("xmin","xmax")],
#ylimにy座標方向の最大・最小値を指定
ylim=sf::st_bbox(obj=rail_jp_xy_fukuoka_sub_buffer)[c("ymin","ymax")]) +
#縮尺記号をプロット
ggspatial::annotation_scale() +
#方位記号をプロット
ggspatial::annotation_north_arrow(location="tr")
上の例では，階級の色分けを自動で設定しましたが，自身で分級値を設定してプロットすることも可能です．今回は，分級値を変数ln_lp_avgの基本統計量を用いて設定します．まず，summary関数を用いて基本統計量を計算します．これを用い，geom_sf関数の引数fillで分級値を指定します．引数fillに渡されているcut関数は，数値ベクトルを任意の境界値で分級するための関数です．引数xには数値ベクトル，breaksには境界値を格納したベクトルを指定します．今回は，xに変数ln_lp_avgを，breaksに各種統計値を指定しています．
自身で階級を設定した場合，階級の色分けの設定はscale_fill_brewer関数を用いています．引数paletteには先ほどと同様Spectralを用います．scale_fill_brewer関数でSpectralを用いると，色の並びが暖色→寒色になってしまうので，引数directionを-1に指定することで，寒色→暖色に並び替えています．
#変数ln_lp_avg（平均対数地価）の基本統計量を計算
ln_lp_avg_sum <- summary(rail_jp_xy_fukuoka_sub_buffer_lp_avg$ln_lp_avg)
ln_lp_avg_sum## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 11.57 12.18 13.08 13.16 14.04 15.27ggplot2::ggplot() +
#福岡市域ポリゴンを白でプロット
ggplot2::geom_sf(data=ct_fukuoka_c_xy_ward,fill="white") +
#「500mバッファを変数ln_lp_avgに基づいて色分けしプロットする」という指定
#色分けの分級値は基本統計量ln_lp_avg_sumに基づいて設定し，階級数は4つにする
ggplot2::geom_sf(data=rail_jp_xy_fukuoka_sub_buffer_lp_avg,
aes(fill=cut(x=ln_lp_avg,breaks=c(ln_lp_avg_sum["Min."],
ln_lp_avg_sum["1st Qu."],
ln_lp_avg_sum["Median"],
ln_lp_avg_sum["3rd Qu."],
ln_lp_avg_sum["Max."]),
#最小値を含める
include.lowest=TRUE))) +
#階級の色分けには「Spectral」というパレットを使う
ggplot2::scale_fill_brewer(palette="Spectral",direction=-1) +
#図の凡例のタイトルを手動で指定する
ggplot2::guides(fill=ggplot2::guide_legend(title="Log Land Price")) +
#xlimにx座標方向の最大・最小値を指定
ggplot2::coord_sf(xlim=sf::st_bbox(obj=rail_jp_xy_fukuoka_sub_buffer)[c("xmin","xmax")],
#ylimにy座標方向の最大・最小値を指定
ylim=sf::st_bbox(obj=rail_jp_xy_fukuoka_sub_buffer)[c("ymin","ymax")]) +
#縮尺記号をプロット
ggspatial::annotation_scale() +
#方位記号をプロット
ggspatial::annotation_north_arrow(location="tr")
処理結果の集計・プロット（mapview関数使用）
mapview関数では，geom_sf関数を用いた場合ほど手の込んだ地図を作ることはできませんが，インタラクティブな操作が可能な地図を作ることができます．まず，福岡市域ポリゴン・500mバッファ・抽出された地価ポイントを重ねたプロットを作成します．空間データの色分けは，引数col.regionsで行います．mapview関数の場合，ポリゴン・ポイント内部の色は引数col.regionsで，それらの枠線の色は引数colで指定します．また，各空間データの凡例は今回は使わないので，引数legendをFALSEにしています．
mapview::mapview(x=ct_fukuoka_c_xy_ward,col.regions="white",legend=FALSE) +
mapview::mapview(x=rail_jp_xy_fukuoka_sub_buffer,col.regions="green",legend=FALSE) +
mapview::mapview(x=lp_fukuoka_xy_c_sub_buffer_500m,col.regions="red",legend=FALSE) 次に，対数地価で色分けしたバッファをプロットします．いくつかの引数を指定しますので，順を追って説明します．引数layer.nameでは，地図に表示される凡例の名前を指定します．引数zcolでは，色分けに用いる変数を指定します．今回は，平均対数地価ln_lp_avgに従って色分けするよう指定しています．引数col.regionsでは，色分けに用いる変数を幾つに分級するかと，それぞれの階級に対応するパレットを何にするかを指定しています．各階級にパレットの色を割り当てる際には，RColorBrewerのbrewer.pal関数を用います．引数nには階級数，nameにはパレット名を指定します．geom_sf関数の場合と同様に，パレットはSpectralを用います．brewer.pal関数をrev関数で囲っていますが，これはSpectralパレットの色の並びを暖色→寒色から，寒色→暖色に並び変えるためのものです．これにより，geom_sf関数の場合と同様な値と色の対応づけができます．
プロットされたバッファの中で，色が赤に近いものが平均対数地価が高いバッファになります．それをクリックすると，変数N05_011（駅名）が「博多」や「天神」であることが分かります．これら2つの駅は，福岡市内の二大ターミナル駅ですので，その周辺で地価が高くなることは，直感に合う結果だと思います．また，それら2つの駅から離れるほど，地価は低くなっていきます．
mapview::mapview(x=rail_jp_xy_fukuoka_sub_buffer_lp_avg,
layer.name="Log Land Price",
zcol="ln_lp_avg",
col.regions=rev(brewer.pal(n=8,name="Spectral")))距離的近さに基づいて空間データの属性を結合（空間結合）
この空間データ操作の用途は，例えば以下の2つです．
- 1．ポリゴン属性の結合：あるポリゴンに含まれる各ポイントに関して，そのポリゴンの属性を付与する．
- 2．最寄り地物情報の付与：2つのポイントデータ（例えば地価ポイントと駅ポイント）がある場合，片方のデータに，もう片方のデータに含まれる最寄りのポイントの属性を付与する．
いずれの空間データ操作も，st_join関数を用いて行うことができます．
ポイントへのポリゴン属性の結合
例として，福岡県内の地価ポイントに，福岡市域ポリゴンの属性を空間的に結合します．引数xでは属性が付与される空間データ（今回は福岡県内の地価ポイント），yでは付与する属性を持つ空間データ（今回は福岡市域ポリゴン）を指定します．すなわち，空間データxにyの属性を付与することになります．引数joinは空間結合の方法を指定するもので，デフォルトではst_intersectsが指定さています．他にもいくつかの結合方法が指定できますが，その一例は，後の「最寄り地物の属性情報の結合・最寄り地物までの距離計算」セクションで解説します．
引数joinをst_intersectsに指定してst_join関数を実行することで，あるポイントが属しているポリゴンの属性が付与されます．どのポリゴンにも属さないポイントでは，属性は付与されません．実行結果の先頭行を見てみると，福岡市域ポリゴン側の変数CSS_NAME（政令区名）がNAになっています．これは，地価データの冒頭にあるレコードは，福岡市ではなく北九州市にあるものだからです．
#地価ポイントに福岡市域ポリゴンの属性を付与する
lp_fukuoka_xy_c_join <- sf::st_join(x=lp_fukuoka_xy,
y=ct_fukuoka_c_xy_ward,
join=st_intersects)
head(lp_fukuoka_xy_c_join)## Simple feature collection with 6 features and 141 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -6571.936 ymin: 96836.62 xmax: 694.2906 ymax: 105154.8
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## 1 000 001 000 001 2022 92300 1.7 1 false
## 2 000 002 000 002 2022 59500 -1.2 1 false
## 3 000 003 000 003 2022 50600 -0.8 1 false
## 4 000 004 000 004 2022 75000 0.4 1 false
## 5 000 005 000 005 2022 26400 -1.1 1 false
## 6 000 006 000 006 2022 33200 -0.6 1 false
## L01_010 L01_011 L01_012 L01_013 L01_014 L01_015 L01_016 L01_017 L01_018
## 1 false false false false false false false false false
## 2 false false false false false false false false false
## 3 false false false false false false false false false
## 4 false false false false false false false false false
## 5 false false false false false false false false false
## 6 false false false false false false false false false
## L01_019 L01_020 L01_021 L01_022 L01_023
## 1 false false false 40101 北九州門司
## 2 false false false 40101 北九州門司
## 3 false false false 40101 北九州門司
## 4 false false false 40101 北九州門司
## 5 false false false 40101 北九州門司
## 6 false false false 40101 北九州門司
## L01_024 L01_025 L01_026
## 1 福岡県 北九州市門司区上馬寄１丁目２番２５ 上馬寄１−２−１６ 223
## 2 福岡県 北九州市門司区東門司２丁目１９５４番３ 東門司２−３−１８ 109
## 3 福岡県 北九州市門司区清滝３丁目１２番１３ 清滝３−１−５ 118
## 4 福岡県 北九州市門司区黄金町３１番８ 黄金町７−２０ 293
## 5 福岡県 北九州市門司区白野江４丁目１９８１番２ 白野江４−２−１５ 219
## 6 福岡県 北九州市門司区大字畑字今在家２１４４番１ _ 231
## L01_027 L01_028 L01_029 L01_030 L01_031 L01_032 L01_033 L01_034 L01_035
## 1 住宅 _ 001 W true true true _ 1.0
## 2 住宅 _ 001 W true true true _ 1.5
## 3 住宅 _ 001 W true true true _ 1.0
## 4 住宅 _ 001 W true true true _ 1.2
## 5 住宅 _ 001 W true false true _ 1.0
## 6 住宅 _ 001 W true false true 台形 1.0
## L01_036 L01_037 L01_038 L01_039 L01_040 L01_041 L01_042 L01_043 L01_044
## 1 1.5 2 0 市道 西 5.0 _ _ _
## 2 1.0 2 0 市道 南 6.0 _ _ _
## 3 1.2 2 0 市道 北 4.5 _ _ _
## 4 1.0 2 0 市道 東 6.0 _ _ _
## 5 1.5 2 0 県道 東 11.0 _ _ _
## 6 1.5 2 0 市道 西 4.6 _ _ _
## L01_045 L01_046 L01_047 L01_048 L01_049
## 1 _ _ 中規模一般住宅が区画整然と建ち並ぶ住宅地域 門司 1900
## 2 _ _ 一般住宅の中に小売店舗も見られる住宅地域 門司港 1400
## 3 _ _ 小規模一般住宅が多い閑静な住宅地域 門司港 900
## 4 _ _ 一般住宅のほかにマンションも多い住宅地域 門司 1100
## 5 _ _ 中規模一般住宅、農家住宅が混在する住宅地域 門司港 7100
## 6 _ _ 中規模一般住宅が建ち並ぶ住宅地域 門司 7500
## L01_050 L01_051 L01_052 L01_053 L01_054 L01_055 L01_056 L01_057 L01_058
## 1 1中専 _ 市街化 0 _ _ 60 200 false
## 2 1住居 準防 市街化 0 _ _ 60 200 false
## 3 1住居 _ 市街化 0 _ _ 60 200 false
## 4 1住居 準防 市街化 0 _ _ 60 200 false
## 5 1低専 _ 市街化 0 _ _ 50 80 false
## 6 1中専 _ 市街化 0 _ _ 60 200 false
## L01_059 L01_060 L01_061 L01_062 L01_063
## 1 true 0000000000000000000000000000001111111111 0 0 0
## 2 false 0000000011111111111111111111111111111111 0 0 0
## 3 false 1111111111111111111111111111111111111111 70600 73500 75500
## 4 false 0000000000001111111111111111111111111111 0 0 0
## 5 false 1111111111111111111111111111111111111111 40000 41500 42500
## 6 false 1111111111111111111111111111111111111111 45000 46500 47600
## L01_064 L01_065 L01_066 L01_067 L01_068 L01_069 L01_070 L01_071 L01_072
## 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 83000 85500 85800 85800
## 3 76500 76500 76500 77300 78100 82500 85000 85500 85500
## 4 0 0 0 0 0 0 0 0 0
## 5 43000 43000 43000 43000 43500 44800 46500 46600 46600
## 6 48000 48000 48700 49700 50500 56000 58000 59000 59000
## L01_073 L01_074 L01_075 L01_076 L01_077 L01_078 L01_079 L01_080 L01_081
## 1 0 0 0 0 0 0 0 0 0
## 2 85800 85800 85800 85800 85800 85800 84700 82800 80500
## 3 85500 85500 85500 85500 85500 85500 84700 84300 82600
## 4 110000 110000 110000 110000 110000 108000 106000 104000 102000
## 5 46600 46600 46600 46600 46600 46600 46200 45800 45300
## 6 59500 60100 60600 61000 61500 61500 60700 60000 58900
## L01_082 L01_083 L01_084 L01_085 L01_086 L01_087 L01_088 L01_089 L01_090
## 1 0 0 0 0 0 0 0 0 0
## 2 78000 75300 72400 69600 67100 65100 62800 61400 60900
## 3 79200 75600 72000 68500 66500 64800 62800 61000 59100
## 4 98000 94300 91000 87200 84100 81200 78400 76900 76000
## 5 44500 42500 40600 38600 36900 35200 33400 31900 30800
## 6 53500 50200 47200 44500 42900 41500 40000 38700 37500
## L01_091 L01_092 L01_093 L01_094 L01_095 L01_096 L01_097 L01_098 L01_099
## 1 84500 84500 84500 84800 85200 85400 87100 89900 90800
## 2 60600 60300 60000 59900 59900 59900 59900 60200 60200
## 3 57300 55900 54700 53800 53100 52400 51800 51400 51000
## 4 75200 74800 74500 74300 74300 74300 74300 74700 74700
## 5 29800 29100 28600 28200 27900 27600 27300 27000 26700
## 6 36300 35500 34800 34300 34000 33800 33600 33500 33400
## L01_100 L01_101 L01_102 L01_103 L01_104
## 1 92300 00000000000000 00000000000000 00000000000000 00000000000000
## 2 59500 00000000000000 00000000000000 00000000000000 00000000000000
## 3 50600 10000000000000 10000000000000 20000000000000 10000000000000
## 4 75000 00000000000000 00000000000000 00000000000000 00000000000000
## 5 26400 10000000000000 10011000000000 10000000000000 10000000000000
## 6 33200 10000000000000 10000000000000 10000000000000 10000000000000
## L01_105 L01_106 L01_107 L01_108 L01_109
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 00000000000000 00000000000000 00000000000000 40000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 5 11000000000000 10000000000000 10000000000000 10000000000000 10000100000000
## 6 10000100000010 10000000000000 10000100000000 10000000000000 10000000000000
## L01_110 L01_111 L01_112 L01_113 L01_114
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000001000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000001000000
## 4 00000000000000 00000000000000 40000000000000 10000000000000 10000001000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000001000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000001000010
## L01_115 L01_116 L01_117 L01_118 L01_119
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000010000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000010000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000001000010 10000000000000 10000000000000
## L01_120 L01_121 L01_122 L01_123 L01_124
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000010 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000010 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000010 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000011 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_125 L01_126 L01_127 L01_128 L01_129
## 1 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 20000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_130 L01_131 L01_132 L01_133 L01_134
## 1 40000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## L01_135 L01_136 L01_137 L01_138 L01_139
## 1 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 2 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 3 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 4 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 5 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 6 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## LID CSS_NAME geometry
## 1 1 <NA> POINT (-6571.936 98823.12)
## 2 2 <NA> POINT (-2522.357 105154.8)
## 3 3 <NA> POINT (-3584.012 104263.2)
## 4 4 <NA> POINT (-5503.101 100850.9)
## 5 5 <NA> POINT (694.2906 104778)
## 6 6 <NA> POINT (-2307.264 96836.62)変数CSS_NAME（政令区名）がNAでないレコードのみ抽出し，レコード数を確認すると293となり，st_intersection関数のレコード数と一致します．
#地価ポイントに福岡市域ポリゴンの属性を付与する
lp_fukuoka_xy_c_join <- sf::st_join(x=lp_fukuoka_xy,
y=ct_fukuoka_c_xy_ward,
join=st_intersects) %>%
#変数CSS_NAMEがNAでないレコードのみ残す
dplyr::filter(!is.na(CSS_NAME))
#レコード数を調べる
dim(lp_fukuoka_xy_c_join)## [1] 293 142最寄り地物の属性情報の結合・最寄り地物までの距離計算
例として，福岡市域の各地価ポイントから最寄りの駅ポイントの属性を付与する作業を実装します．今回は，最寄り駅ポイントの属性を付与するだけでなく，最寄り駅ポイントまでの直線距離を計算する作業も合わせて行います．まず，設置終了時点の変数N05_005eが9999年（即ち現在も存続していることを表す）となっている駅ポイントのみを抽出します．あわせて，抽出された駅ポイントにID変数（SID）を付与します．
rail_jp_xy_exist <- rail_jp_xy %>%
#現在も存続する駅のみ抽出
dplyr::filter(N05_005e=="9999") %>%
#駅IDを追加
dplyr::mutate(SID=dplyr::row_number())直線距離の計算を行うためには，各ポイントデータのgeometry属性を，四則演算などが可能な形式の変数として格納し直す作業が必要です．変数geometryからXY座標を取り出す際には，st_coordinates関数を用います．
#地価ポイントのXY座標を取り出し
lp_coords <- sf::st_coordinates(x=lp_fukuoka_xy_c)
head(lp_coords)## X Y
## 1 -56754.69 65075.79
## 2 -58150.03 65083.26
## 3 -56564.00 64307.29
## 4 -57971.52 66100.34
## 5 -57697.73 63836.13
## 6 -57460.85 64153.45#駅ポイントのXY座標を取り出し
rail_coords <- sf::st_coordinates(x=rail_jp_xy_exist)
head(rail_coords)## X Y
## 1 914730.2 1257807
## 2 884284.8 1245446
## 3 907786.9 1256819
## 4 889864.8 1249625
## 5 911706.7 1260321
## 6 882745.0 1238074続いて，取り出したXY座標を，各ポイントデータへ新たな変数として追加します．
#XY座標を地価ポイントに変数として追加
lp_fukuoka_xy_c_coords <- lp_fukuoka_xy_c %>%
dplyr::mutate(X_lp=lp_coords[,"X"],Y_lp=lp_coords[,"Y"])
#XY座標を駅ポイントに変数として追加
rail_jp_xy_exist_coords <- rail_jp_xy_exist %>%
dplyr::mutate(X_sta=rail_coords[,"X"],Y_sta=rail_coords[,"Y"])XY座標の作成から変数として追加までの手続きを，mutate関数の中で一括で行うことも可能です．少し複雑に見えますが，以下のプログラムによって実行できます．
#XY座標を地価ポイントに変数として追加
lp_fukuoka_xy_c_coords <- lp_fukuoka_xy_c %>%
dplyr::mutate(X_lp=sf::st_coordinates(.)[,"X"],Y_lp=sf::st_coordinates(.)[,"Y"])
head(lp_fukuoka_xy_c_coords)## Simple feature collection with 6 features and 143 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58150.03 ymin: 63836.13 xmax: -56564 ymax: 66100.34
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## 336 000 001 000 001 2022 575000 8.9 1 false
## 337 000 002 000 002 2022 870000 5.5 1 false
## 338 000 003 000 003 2022 324000 12.1 1 false
## 339 000 004 000 004 2022 488000 6.6 1 false
## 340 000 005 000 005 2022 461000 7.7 1 false
## 341 000 006 000 006 2022 300000 12.8 1 false
## L01_010 L01_011 L01_012 L01_013 L01_014 L01_015 L01_016 L01_017 L01_018
## 336 false false false false false false false false false
## 337 false false false false false false false false false
## 338 false false false false false false false false false
## 339 false false false false false false false false false
## 340 false false false false false false false false false
## 341 false false false false false false false false false
## L01_019 L01_020 L01_021 L01_022 L01_023
## 336 false false false 40133 福岡中央
## 337 false false false 40133 福岡中央
## 338 false false false 40133 福岡中央
## 339 false false false 40133 福岡中央
## 340 false false false 40133 福岡中央
## 341 false false false 40133 福岡中央
## L01_024 L01_025 L01_026
## 336 福岡県 福岡市中央区赤坂２丁目２区５４番 赤坂２−２−１１ 737
## 337 福岡県 福岡市中央区大濠１丁目１３区１３３番 大濠１−１３−２６ 718
## 338 福岡県 福岡市中央区警固３丁目１０３番 _ 225
## 339 福岡県 福岡市中央区荒戸３丁目２９７番１ 荒戸３−５−５１ 1085
## 340 福岡県 福岡市中央区六本松４丁目６１番 六本松４−５−１８ 330
## 341 福岡県 福岡市中央区谷１丁目２９３番 谷１−９−２１ 161
## L01_027 L01_028 L01_029 L01_030 L01_031 L01_032 L01_033 L01_034
## 336 住宅,その他 共同住宅 001 RC true true true _
## 337 住宅 _ 001 W true true true _
## 338 空地 _ 004 その他 true true true 台形
## 339 住宅,その他 共同住宅 001 RC true true true _
## 340 住宅,その他 共同住宅 001 RC true true true _
## 341 住宅 _ 001 W true true true _
## L01_035 L01_036 L01_037 L01_038 L01_039 L01_040 L01_041 L01_042 L01_043
## 336 1 2.0 6 0 市道 南 5.9 _ _
## 337 1 1.0 2 0 市道 東 5.4 _ _
## 338 1 1.5 0 0 市道 北 4.0 _ _
## 339 1 2.0 7 0 市道 北 6.8 _ _
## 340 1 1.0 3 0 市道 南 7.0 _ _
## 341 1 2.0 2 0 市道 北東 5.6 _ _
## L01_044 L01_045 L01_046 L01_047
## 336 _ _ _ マンション等が建ち並ぶ都心に近い住宅地域
## 337 側道 南 _ 大規模一般住宅、マンションが建ち並ぶ住宅地域
## 338 _ _ _ 一般住宅の中に共同住宅が見られる住宅地域
## 339 _ _ _ 中規模住宅、マンション等が混在する住宅地域
## 340 _ _ _ 共同住宅の中に一般住宅も見られる住宅地域
## 341 _ _ _ 中規模住宅、アパート等が多い既成住宅地域
## L01_048 L01_049 L01_050 L01_051 L01_052 L01_053 L01_054 L01_055 L01_056
## 336 赤坂 560 2住居 準防 市街化 0 _ _ 60
## 337 唐人町 750 1住居 準防 市街化 0 _ _ 60
## 338 桜坂 420 1住居 準防 市街化 0 _ _ 60
## 339 大濠公園 800 1住居 準防 市街化 0 _ _ 60
## 340 六本松 500 1住居 _ 市街化 0 _ _ 60
## 341 六本松 430 1住居 _ 市街化 0 _ _ 60
## L01_057 L01_058 L01_059 L01_060 L01_061
## 336 200 false false 0000000000000000001111111111111111111111 0
## 337 200 false false 1111111111111111111111111111111111111111 252000
## 338 200 false false 0000000000000000000000000000000111111111 0
## 339 200 false false 0000000111111111111111111111111111111111 0
## 340 200 false true 0000001111111111111111111111111111111111 0
## 341 200 false false 0000000000000000000000000000001111111111 0
## L01_062 L01_063 L01_064 L01_065 L01_066 L01_067 L01_068 L01_069 L01_070
## 336 0 0 0 0 0 0 0 0 0
## 337 267000 273000 277000 305000 460000 550000 650000 750000 750000
## 338 0 0 0 0 0 0 0 0 0
## 339 0 0 0 0 0 0 450000 519000 519000
## 340 0 0 0 0 0 270000 300000 365000 365000
## 341 0 0 0 0 0 0 0 0 0
## L01_071 L01_072 L01_073 L01_074 L01_075 L01_076 L01_077 L01_078 L01_079
## 336 0 0 0 0 0 0 0 0 260000
## 337 650000 600000 555000 520000 475000 443000 423000 410000 400000
## 338 0 0 0 0 0 0 0 0 0
## 339 455000 385000 340000 308000 285000 270000 260000 248000 240000
## 340 328000 315000 275000 245000 227000 216000 212000 207000 203000
## 341 0 0 0 0 0 0 0 0 0
## L01_080 L01_081 L01_082 L01_083 L01_084 L01_085 L01_086 L01_087 L01_088
## 336 257000 253000 251000 252000 259000 300000 331000 308000 293000
## 337 393000 393000 393000 400000 420000 525000 603000 535000 494000
## 338 0 0 0 0 0 0 0 0 0
## 339 228000 220000 214000 214000 221000 284000 328000 310000 290000
## 340 195000 187000 181000 181000 182000 214000 235000 223000 215000
## 341 0 0 0 0 0 0 0 0 0
## L01_089 L01_090 L01_091 L01_092 L01_093 L01_094 L01_095 L01_096 L01_097
## 336 291000 292000 305000 321000 340000 359000 391000 425000 463000
## 337 494000 501000 518000 534000 566000 605000 650000 698000 756000
## 338 0 0 0 178000 186000 197000 208000 219000 232000
## 339 281000 282000 289000 299000 311000 329000 348000 373000 389000
## 340 208000 213000 217000 229000 242000 262000 288000 316000 346000
## 341 0 0 159000 164000 173000 187000 201000 212000 227000
## L01_098 L01_099 L01_100 L01_101 L01_102 L01_103
## 336 504000 528000 575000 00000000000000 00000000000000 00000000000000
## 337 800000 825000 870000 10000000000000 10000000000000 10000000000000
## 338 266000 289000 324000 00000000000000 00000000000000 00000000000000
## 339 438000 458000 488000 00000000000000 00000000000000 00000000000000
## 340 394000 428000 461000 00000000000000 00000000000000 00000000000000
## 341 249000 266000 300000 00000000000000 00000000000000 00000000000000
## L01_104 L01_105 L01_106 L01_107 L01_108
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 00000000000000 00000000000000 00000000000000 40000000000000 10010000000000
## 340 00000000000000 00000000000000 40000000000000 10000000000000 10010000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_109 L01_110 L01_111 L01_112 L01_113
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_114 L01_115 L01_116 L01_117 L01_118
## 336 00000000000000 00000000000000 00000000000000 00000000000000 40000000000000
## 337 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_119 L01_120 L01_121 L01_122 L01_123
## 336 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000010 10000000000000 10000000000000 10000010000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_124 L01_125 L01_126 L01_127 L01_128
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_129 L01_130 L01_131 L01_132 L01_133
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 338 00000000000000 00000000000000 40000000000000 10000000000000 10000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 40000000000000 10000000000000 10000000000000 10000000000000
## L01_134 L01_135 L01_136 L01_137 L01_138
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 10000000000000 10000000000000 10000010000000 10000000000000 11011000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 10000000000000 10000000000000 10000010000000 10000000000000 10000000000000
## L01_139 LID CSS_NAME geometry X_lp Y_lp
## 336 10000000000000 336 中央区 POINT (-56754.69 65075.79) -56754.69 65075.79
## 337 10000000000000 337 中央区 POINT (-58150.03 65083.26) -58150.03 65083.26
## 338 10000000000000 338 中央区 POINT (-56564 64307.29) -56564.00 64307.29
## 339 10000000000000 339 中央区 POINT (-57971.52 66100.34) -57971.52 66100.34
## 340 10000000000000 340 中央区 POINT (-57697.73 63836.13) -57697.73 63836.13
## 341 10000000000000 341 中央区 POINT (-57460.85 64153.45) -57460.85 64153.45#XY座標を駅ポイントに変数として追加
rail_jp_xy_exist_coords <- rail_jp_xy_exist %>%
dplyr::mutate(X_sta=sf::st_coordinates(.)[,"X"],Y_sta=sf::st_coordinates(.)[,"Y"])
head(rail_jp_xy_exist_coords)## Simple feature collection with 6 features and 14 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 882745 ymin: 1238074 xmax: 914730.2 ymax: 1260321
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## N05_001 N05_002 N05_003 N05_004 N05_005b N05_005e
## 1 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 2 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 3 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 4 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## 5 2 函館線 北海道旅客鉄道（旧国鉄） 1911 1950 9999
## 6 2 函館線 北海道旅客鉄道（旧国鉄） 1898 1950 9999
## N05_006 N05_007 N05_008 N05_009 N05_011 geometry SID
## 1 EB03_11218106 <NA> <NA> <NA> 旭川 POINT (914730.2 1257807) 1
## 2 EB03_11218099 <NA> <NA> <NA> 妹背牛 POINT (884284.8 1245446) 2
## 3 EB03_11218103 <NA> <NA> <NA> 伊納 POINT (907786.9 1256819) 3
## 4 EB03_11218100 <NA> <NA> <NA> 深川 POINT (889864.8 1249625) 4
## 5 EB03_11218105 <NA> <NA> <NA> 近文 POINT (911706.7 1260321) 5
## 6 EB03_11218098 <NA> <NA> <NA> 江部乙 POINT (882745 1238074) 6
## X_sta Y_sta
## 1 914730.2 1257807
## 2 884284.8 1245446
## 3 907786.9 1256819
## 4 889864.8 1249625
## 5 911706.7 1260321
## 6 882745.0 1238074ここで，st_join関数を用いて，各地価ポイントに最寄りの駅ポイントの属性を結合します．ポリゴン属性結合の場合とは異なり，引数joinにはst_nearest_featureを指定することに注意してください．
#各地価ポイントから最寄りの駅ポイントの属性を付与
lp_fukuoka_xy_c_near_sta <- sf::st_join(x=lp_fukuoka_xy_c_coords,
y=rail_jp_xy_exist_coords,
join=st_nearest_feature)最後に，最寄り駅ポイントまでの直線距離を計算してみましょう．座標はm単位で定義されていますので，単に，地価ポイント・駅ポイントの座標からユークリッド距離（直線距離）を計算するだけです．後の分析のために，距離と地価の対数値も別の変数として追加します．
lp_fukuoka_xy_c_near_sta_dist <- lp_fukuoka_xy_c_near_sta %>%
#地価ポイントと駅ポイントの間の直線距離
dplyr::mutate(d_sta=sqrt((X_lp-X_sta)^2+(Y_lp-Y_sta)^2)) %>%
#直線距離の対数値も変数として追加
dplyr::mutate(ln_d_sta=log(d_sta)) %>%
#地価の対数値も変数として追加
dplyr::mutate(ln_lp=log(L01_006))
head(lp_fukuoka_xy_c_near_sta_dist)## Simple feature collection with 6 features and 160 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58150.03 ymin: 63836.13 xmax: -56564 ymax: 66100.34
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## 336 000 001 000 001 2022 575000 8.9 1 false
## 337 000 002 000 002 2022 870000 5.5 1 false
## 338 000 003 000 003 2022 324000 12.1 1 false
## 339 000 004 000 004 2022 488000 6.6 1 false
## 340 000 005 000 005 2022 461000 7.7 1 false
## 341 000 006 000 006 2022 300000 12.8 1 false
## L01_010 L01_011 L01_012 L01_013 L01_014 L01_015 L01_016 L01_017 L01_018
## 336 false false false false false false false false false
## 337 false false false false false false false false false
## 338 false false false false false false false false false
## 339 false false false false false false false false false
## 340 false false false false false false false false false
## 341 false false false false false false false false false
## L01_019 L01_020 L01_021 L01_022 L01_023
## 336 false false false 40133 福岡中央
## 337 false false false 40133 福岡中央
## 338 false false false 40133 福岡中央
## 339 false false false 40133 福岡中央
## 340 false false false 40133 福岡中央
## 341 false false false 40133 福岡中央
## L01_024 L01_025 L01_026
## 336 福岡県 福岡市中央区赤坂２丁目２区５４番 赤坂２−２−１１ 737
## 337 福岡県 福岡市中央区大濠１丁目１３区１３３番 大濠１−１３−２６ 718
## 338 福岡県 福岡市中央区警固３丁目１０３番 _ 225
## 339 福岡県 福岡市中央区荒戸３丁目２９７番１ 荒戸３−５−５１ 1085
## 340 福岡県 福岡市中央区六本松４丁目６１番 六本松４−５−１８ 330
## 341 福岡県 福岡市中央区谷１丁目２９３番 谷１−９−２１ 161
## L01_027 L01_028 L01_029 L01_030 L01_031 L01_032 L01_033 L01_034
## 336 住宅,その他 共同住宅 001 RC true true true _
## 337 住宅 _ 001 W true true true _
## 338 空地 _ 004 その他 true true true 台形
## 339 住宅,その他 共同住宅 001 RC true true true _
## 340 住宅,その他 共同住宅 001 RC true true true _
## 341 住宅 _ 001 W true true true _
## L01_035 L01_036 L01_037 L01_038 L01_039 L01_040 L01_041 L01_042 L01_043
## 336 1 2.0 6 0 市道 南 5.9 _ _
## 337 1 1.0 2 0 市道 東 5.4 _ _
## 338 1 1.5 0 0 市道 北 4.0 _ _
## 339 1 2.0 7 0 市道 北 6.8 _ _
## 340 1 1.0 3 0 市道 南 7.0 _ _
## 341 1 2.0 2 0 市道 北東 5.6 _ _
## L01_044 L01_045 L01_046 L01_047
## 336 _ _ _ マンション等が建ち並ぶ都心に近い住宅地域
## 337 側道 南 _ 大規模一般住宅、マンションが建ち並ぶ住宅地域
## 338 _ _ _ 一般住宅の中に共同住宅が見られる住宅地域
## 339 _ _ _ 中規模住宅、マンション等が混在する住宅地域
## 340 _ _ _ 共同住宅の中に一般住宅も見られる住宅地域
## 341 _ _ _ 中規模住宅、アパート等が多い既成住宅地域
## L01_048 L01_049 L01_050 L01_051 L01_052 L01_053 L01_054 L01_055 L01_056
## 336 赤坂 560 2住居 準防 市街化 0 _ _ 60
## 337 唐人町 750 1住居 準防 市街化 0 _ _ 60
## 338 桜坂 420 1住居 準防 市街化 0 _ _ 60
## 339 大濠公園 800 1住居 準防 市街化 0 _ _ 60
## 340 六本松 500 1住居 _ 市街化 0 _ _ 60
## 341 六本松 430 1住居 _ 市街化 0 _ _ 60
## L01_057 L01_058 L01_059 L01_060 L01_061
## 336 200 false false 0000000000000000001111111111111111111111 0
## 337 200 false false 1111111111111111111111111111111111111111 252000
## 338 200 false false 0000000000000000000000000000000111111111 0
## 339 200 false false 0000000111111111111111111111111111111111 0
## 340 200 false true 0000001111111111111111111111111111111111 0
## 341 200 false false 0000000000000000000000000000001111111111 0
## L01_062 L01_063 L01_064 L01_065 L01_066 L01_067 L01_068 L01_069 L01_070
## 336 0 0 0 0 0 0 0 0 0
## 337 267000 273000 277000 305000 460000 550000 650000 750000 750000
## 338 0 0 0 0 0 0 0 0 0
## 339 0 0 0 0 0 0 450000 519000 519000
## 340 0 0 0 0 0 270000 300000 365000 365000
## 341 0 0 0 0 0 0 0 0 0
## L01_071 L01_072 L01_073 L01_074 L01_075 L01_076 L01_077 L01_078 L01_079
## 336 0 0 0 0 0 0 0 0 260000
## 337 650000 600000 555000 520000 475000 443000 423000 410000 400000
## 338 0 0 0 0 0 0 0 0 0
## 339 455000 385000 340000 308000 285000 270000 260000 248000 240000
## 340 328000 315000 275000 245000 227000 216000 212000 207000 203000
## 341 0 0 0 0 0 0 0 0 0
## L01_080 L01_081 L01_082 L01_083 L01_084 L01_085 L01_086 L01_087 L01_088
## 336 257000 253000 251000 252000 259000 300000 331000 308000 293000
## 337 393000 393000 393000 400000 420000 525000 603000 535000 494000
## 338 0 0 0 0 0 0 0 0 0
## 339 228000 220000 214000 214000 221000 284000 328000 310000 290000
## 340 195000 187000 181000 181000 182000 214000 235000 223000 215000
## 341 0 0 0 0 0 0 0 0 0
## L01_089 L01_090 L01_091 L01_092 L01_093 L01_094 L01_095 L01_096 L01_097
## 336 291000 292000 305000 321000 340000 359000 391000 425000 463000
## 337 494000 501000 518000 534000 566000 605000 650000 698000 756000
## 338 0 0 0 178000 186000 197000 208000 219000 232000
## 339 281000 282000 289000 299000 311000 329000 348000 373000 389000
## 340 208000 213000 217000 229000 242000 262000 288000 316000 346000
## 341 0 0 159000 164000 173000 187000 201000 212000 227000
## L01_098 L01_099 L01_100 L01_101 L01_102 L01_103
## 336 504000 528000 575000 00000000000000 00000000000000 00000000000000
## 337 800000 825000 870000 10000000000000 10000000000000 10000000000000
## 338 266000 289000 324000 00000000000000 00000000000000 00000000000000
## 339 438000 458000 488000 00000000000000 00000000000000 00000000000000
## 340 394000 428000 461000 00000000000000 00000000000000 00000000000000
## 341 249000 266000 300000 00000000000000 00000000000000 00000000000000
## L01_104 L01_105 L01_106 L01_107 L01_108
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 00000000000000 00000000000000 00000000000000 40000000000000 10010000000000
## 340 00000000000000 00000000000000 40000000000000 10000000000000 10010000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_109 L01_110 L01_111 L01_112 L01_113
## 336 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_114 L01_115 L01_116 L01_117 L01_118
## 336 00000000000000 00000000000000 00000000000000 00000000000000 40000000000000
## 337 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000001000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_119 L01_120 L01_121 L01_122 L01_123
## 336 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000010 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000010 10000000000000 10000000000000 10000010000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_124 L01_125 L01_126 L01_127 L01_128
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 00000000000000 00000000000000 00000000000000 00000000000000
## L01_129 L01_130 L01_131 L01_132 L01_133
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000010000000 10000000000000
## 338 00000000000000 00000000000000 40000000000000 10000000000000 10000000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 00000000000000 40000000000000 10000000000000 10000000000000 10000000000000
## L01_134 L01_135 L01_136 L01_137 L01_138
## 336 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 337 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 338 10000000000000 10000000000000 10000010000000 10000000000000 11011000000000
## 339 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 340 10000000000000 10000000000000 10000000000000 10000000000000 10000000000000
## 341 10000000000000 10000000000000 10000010000000 10000000000000 10000000000000
## L01_139 LID CSS_NAME X_lp Y_lp N05_001 N05_002 N05_003
## 336 10000000000000 336 中央区 -56754.69 65075.79 3 1号線 福岡市
## 337 10000000000000 337 中央区 -58150.03 65083.26 3 1号線 福岡市
## 338 10000000000000 338 中央区 -56564.00 64307.29 3 3号線 福岡市
## 339 10000000000000 339 中央区 -57971.52 66100.34 3 1号線 福岡市
## 340 10000000000000 340 中央区 -57697.73 63836.13 3 3号線 福岡市
## 341 10000000000000 341 中央区 -57460.85 64153.45 3 3号線 福岡市
## N05_004 N05_005b N05_005e N05_006 N05_007 N05_008 N05_009 N05_011
## 336 1981 1981 9999 EB03_22318007 <NA> <NA> <NA> 赤坂
## 337 1981 1981 9999 EB03_22318005 <NA> <NA> <NA> 唐人町
## 338 2005 2005 9999 EB03_22320005 <NA> <NA> <NA> 桜坂
## 339 1981 1981 9999 EB03_22318006 <NA> <NA> <NA> 大濠公園
## 340 2005 2005 9999 EB03_22320016 <NA> <NA> <NA> 六本松
## 341 2005 2005 9999 EB03_22320016 <NA> <NA> <NA> 六本松
## SID X_sta Y_sta geometry d_sta ln_d_sta
## 336 9062 -56571.61 65488.53 POINT (-56754.69 65075.79) 451.5216 6.112623
## 337 9061 -58480.76 65619.18 POINT (-58150.03 65083.26) 629.7508 6.445324
## 338 9079 -56893.54 64184.46 POINT (-56564 64307.29) 351.6852 5.862736
## 339 9067 -57582.74 65597.68 POINT (-57971.52 66100.34) 635.4720 6.454368
## 340 9091 -57797.96 64228.67 POINT (-57697.73 63836.13) 405.1399 6.004233
## 341 9091 -57797.96 64228.67 POINT (-57460.85 64153.45) 345.4008 5.844705
## ln_lp
## 336 13.26213
## 337 13.67625
## 338 12.68850
## 339 13.09807
## 340 13.04115
## 341 12.61154参考までに，最寄り駅までの距離と地価の関係をプロットしてみましょう．plot関数（Rにデフォルトで用意されています）を使って散布図を描いてみます．引数xにx軸に持ってきたい変数（数値ベクトル），yにy軸に持ってきたい変数を指定します．最寄り駅までの距離が長くなるほど，地価は低くなる傾向がみられます．
plot(x=lp_fukuoka_xy_c_near_sta_dist$ln_lp,y=lp_fukuoka_xy_c_near_sta_dist$ln_d_sta)
デフォルトのplot関数による散布図のプロットは，手軽に作ることができる反面，細かい調整があまり利きません．ggplot2を用いると，よりフレキシブルに散布図を作ることができます．まず，ggplot2のgeom_point関数を用いて，シンプルな散布図をプロットします．引数dataにはデータセット，aesで囲んだ引数xにはx軸にしたい変数，yにはy軸にしたい変数を指定します．
ggplot2::ggplot() +
#x軸を変数ln_d_sta（最寄り駅までの距離の対数），y軸をln_lp（対数地価）とした点プロット
ggplot2::geom_point(data=lp_fukuoka_xy_c_near_sta_dist,aes(x=ln_d_sta,y=ln_lp))
次に，グループ毎に点の色を変えて散布図をプロットしてみます．ここでは変数CSS_NAME（政令区）をグループ変数としたものを作成します．先ほどのプログラムから2つ引数が増えており，引数groupにグループ変数として使いたい変数を，colorに色分けの基準に使いたい変数を指定します．MacOSを使う場合，ggplot2の関数群はそのままでは日本語を受け付けてくれないので，さらにtheme_gray関数を用いて日本語文字列（ヒラギノ角ゴ）を使うことを明示します．
ggplot2::ggplot() +
#変数CSS_NAMEを用いてグループ化＆色分けした点プロット
ggplot2::geom_point(data=lp_fukuoka_xy_c_near_sta_dist,aes(x=ln_d_sta,y=ln_lp,group=CSS_NAME,color=CSS_NAME)) +
#日本語文字を使うことを明示（Windowsで実行する際は不要）
ggplot2::theme_gray(base_family="HiraKakuPro-W3")
最寄り地物へのライン生成
各地価ポイントからの最寄駅を視覚的に把握する為の方法として，これら2つの空間データの間にラインを生成し可視化します．ラインの生成には，nngeoのst_connect関数を用います．引数xにはラインの発生元となる空間データ，yには最寄地物の候補となる空間データ，kには何番目に近い地物までのラインを生成するかを指定します．今回はk=1（1番近い地物までのラインを生成）としていますが，これを仮にk=3とすると，1番近い地物だけでなく，2・3番目に近い地物までのラインも生成できます．
#各地価ポイントから最寄りの駅ポイントへのラインを生成
lp_fukuoka_xy_c_near_sta_nn <- nngeo::st_connect(x=lp_fukuoka_xy_c,
y=rail_jp_xy_exist,
k=1)## Calculating nearest IDs
## Calculating lineshead(lp_fukuoka_xy_c_near_sta_nn)## Geometry set for 6 features
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -58480.76 ymin: 63836.13 xmax: -56564 ymax: 66100.34
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## First 5 geometries:生成されたデータはsfc形式で，そのままでも地図上にプロットすることができますが，始点・終点の属性は持たないので，もう一手間加えて属性の結合を行います．結合のキー変数となる最寄り駅のIDは，nngeoのst_nn関数を用いて取得できます．このIDは，最寄地物の候補となる空間データ内で一貫した（1から始まって途中で番号が飛んだりしない）ものである必要があります．今回駅ポイント内に振られているID（SID）は，この要件を満たします．st_nn関数の実行結果はlist形式ですが，unlist関数でベクトルに変換可能です．
#各地価ポイントから最寄りの駅のIDをlist形式で取得可能
nngeo::st_nn(x=lp_fukuoka_xy_c,y=rail_jp_xy_exist) %>%
#ベクトルに変換
unlist## [1] 9062 9061 9079 9067 9091 9091 9061 9061 9079 9091 9061 9079 9188 9061 9079
## [16] 9085 9079 9081 9091 9085 9067 9085 9064 9188 9062 9088 9087 9189 9188 9188
## [31] 9064 9064 9067 9062 9062 9067 9079 9085 9091 9067 9062 9064 9061 9064 9067
## [46] 9079 9188 9087 9088 9188 9190 8960 9197 9180 9197 9197 9180 9190 9180 9190
## [61] 9083 8960 9188 8960 8960 9197 9197 9180 9079 9197 9180 8395 9188 9190 9197
## [76] 9197 9180 9083 9190 9083 9083 9190 9197 9197 8960 9180 9190 9197 9180 9197
## [91] 9180 9180 9059 9216 8395 9063 9084 8376 9216 8385 8508 9216 8376 9073 9216
## [106] 8508 8376 8959 8397 8959 8397 9073 8959 9060 8959 8376 9066 9073 9059 9216
## [121] 9075 9073 9060 8959 9073 8395 9059 9216 9060 8395 9084 9059 9059 8508 8376
## [136] 9089 9089 9081 9083 9086 9083 9086 9086 9078 9090 9089 9081 9083 9083 9089
## [151] 9083 9083 9089 9083 9065 9068 9058 8644 9068 9058 9090 9068 9065 9090 9080
## [166] 9082 9080 9078 9082 9092 9092 9092 9080 9068 9077 9065 9092 9058 9078 9068
## [181] 9058 9065 9078 9082 9065 9065 9078 9058 9068 9081 9078 9090 9082 9082 9078
## [196] 9092 9068 9092 9065 9058 9081 9068 9082 9082 9065 9065 9065 9163 8488 9166
## [211] 8392 8485 8489 8494 8499 8450 8392 8394 8496 8489 9169 8399 8488 8399 8452
## [226] 9165 9166 8452 8484 8487 8496 9167 8450 9169 8488 9163 8450 9074 8488 9168
## [241] 8484 8499 9166 9164 9070 9076 8496 9169 9072 9072 9169 9167 9076 8487 9072
## [256] 8506 9161 9070 9070 9072 9074 8506 8506 8452 9068 8644 8644 8643 9077 8641
## [271] 8642 8642 8643 8643 8642 9077 8643 8644 8644 8644 8644 8644 9077 9077 9092
## [286] 8639 8643 8641 8644 9077 8644 8642 8639地価ポイントと最寄り駅ポイントの間にラインを生成し，各ポイントの属性を結合する一連の手続きを実装すると，以下のようになります．この方法で生成したラインの長さをsfのst_length関数で取得し，最寄り駅までの距離を求めることもできます（ただし，距離計算に際して新たな空間データを作成しない分，st_join関数を用いた前掲の方法の方が計算負荷は格段に小さいです）．
lp_fukuoka_xy_c_near_sta_line <- lp_fukuoka_xy_c %>%
#地価ポイントから最寄りの駅ポイントに直線を生成
dplyr::mutate(geometry=nngeo::st_connect(x=.,y=rail_jp_xy_exist)) %>%
#直線が引かれた先の駅のIDを新たな変数として追加
dplyr::mutate(SID=unlist(nngeo::st_nn(x=.,y=rail_jp_xy_exist))) %>%
#駅ポイントからgeometryを削除した上で，駅属性を駅IDをキーとして結合
dplyr::left_join(y=sf::st_set_geometry(x=rail_jp_xy_exist,NULL),by="SID") %>%
#生成した直線の長さを計算することで，最寄り駅までの距離が計算可能
dplyr::mutate(d_sta=sf::st_length(x=.))## Calculating nearest IDs
## Calculating lines
##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%可視化に先立ち，駅ポイント全体から，地価ポイントの最寄り駅と判定されたものだけを残します．
rail_jp_xy_exist_nn <- rail_jp_xy_exist %>%
#最寄り駅情報結合済みの地価ポイント内にある，通し番号変数SID（駅名）と値が一致するレコードのみ保持
dplyr::filter(SID%in%unique(lp_fukuoka_xy_c_near_sta_line$SID))mapview関数を用いて，地価ポイントデータ・駅ポイントデータ・地価ポイント-駅ラインデータを可視化します．
#地価ポイントデータ
mapview::mapview(x=lp_fukuoka_xy_c,
#ポイントの色は紫
col.region="purple",
#不透明度100%
alpha.regions=1,
#レイヤ名を「Obs Point」と設定
layer.name="Obs Point") +
#駅ポイントデータ
mapview::mapview(x=rail_jp_xy_exist_nn,
#ポイントの色は緑
col.region="green",
#不透明度100%
alpha.regions=1,
#レイヤ名を「Station」と設定
layer.name="Station") +
#地価ポイント-最寄駅間ラインデータ
mapview::mapview(x=lp_fukuoka_xy_c_near_sta_line,
#ラインの色は赤
color="red",
#レイヤ名を「Station-Obs Point」と設定
layer.name="Station-Obs Point")テキストデータ-空間データ間の変換
この空間データ操作では，csv等のテキスト形式で与えられた座標データがある場合に，それを空間データへと変換する方法を示します．ここでは例として，冒頭で読み込んだ，位置座標付きの福岡市周辺の店舗データfuk_p_lonlatを空間データに変換します．ポイントデータ作成には，sfのst_as_sf関数を用います．引数coordsには読み込んだ店舗データ内の座標を表す変数名を（X座標，Y座標の順），crsには座標系を指定します．店舗データの緯度経度座標はJGD2011（EPSG:6668）で定義されていますので，それを指定します．
fuk_p_lonlat <- fuk_p_lonlat %>%
#テキストデータをポイントデータに変換
sf::st_as_sf(coords=c("lon","lat"),crs=sf::st_crs(6668))最後に，他の空間データと座標系を揃えるため，JGD2000・平面直角座標第2系に投影変換します．
fuk_p_lonlat_xy <- fuk_p_lonlat %>%
#JGD2000に投影変換
sf::st_transform(sf::st_crs(4612)) %>%
#平面直角座標第2系に投影変換
sf::st_transform(sf::st_crs(2444))作成されたポイントデータを，mapview関数を用いてプロットすると以下のようになります．店舗は博多駅周辺や天神駅周辺といった福岡市の二大都心に多く集まっていることが分かります．
mapview::mapview(x=fuk_p_lonlat_xy)作成したポイントデータを元に，新たな変数を作成します．具体的には，各地価観測地点から1km圏内に含まれる店舗数をカウントします．そのためにまず，地価ポイントから1kmバッファを生成します．その際，通し番号変数LID以外の変数は使用しないので削除します．
lp_fukuoka_xy_c_buffer <- lp_fukuoka_xy_c_near_sta_dist %>%
dplyr::select(LID) %>%
sf::st_buffer(dist=1000)ここで，1kmバッファと店舗ポイントデータとの共通部分を取ります．
lp_fukuoka_xy_c_buffer_shop <- sf::st_intersection(x=lp_fukuoka_xy_c_buffer,
y=fuk_p_lonlat_xy)最後に，地価観測地点毎に1km圏内にある店舗の数を集計し，地価ポイントに集計結果を左結合します．結合に先立ち，1kmバッファを用いて抽出された店舗ポイントデータからgeometry属性を削除した上で，地価ポイントの通し番号変数LID毎にグループ化した上でレコード数をカウントし，結果を変数nshopとして追加します．
lp_fukuoka_xy_c_buffer_shop_tbl <- lp_fukuoka_xy_c_buffer_shop %>%
#geometry属性を削除
sf::st_set_geometry(NULL) %>%
#LIDでグループ化
dplyr::group_by(LID) %>%
#レコード数を集計
dplyr::summarise(nshop=n()) %>%
#グループ化情報を削除
dplyr::ungroup()通し番号変数LIDをキー変数として，集計結果を地価ポイントに左結合します．その際，店舗数が結合できなかった（変数nshopがNAの）レコードについては，0に置換します．また，後の分析のため，変数nshopに1足したものの対数値を別の変数として追加します（0の対数値は定義できないためです）．
lp_fukuoka_xy_c_near_sta_dist_shop <- lp_fukuoka_xy_c_near_sta_dist %>%
#集計結果を結合
dplyr::left_join(y=lp_fukuoka_xy_c_buffer_shop_tbl,by="LID") %>%
#変数nshopがNAでなければnshopのそのままの値を，さもなくば0とする
dplyr::mutate(nshop=ifelse(!is.na(nshop),nshop,0)) %>%
#変数nshop+1の対数値を新たな変数として追加
dplyr::mutate(ln_nshop=log(1+nshop))最寄り駅までの距離の場合と同様，1km圏内店舗数と地価の関係をプロットしてみましょう．plot関数を使って散布図を描いてみます．店舗数が多いほど，地価が高くなっている傾向が見て取れます．
plot(x=lp_fukuoka_xy_c_near_sta_dist_shop$ln_lp,y=lp_fukuoka_xy_c_near_sta_dist_shop$ln_nshop)
生成したデータを用いた重回帰分析
説明変数の作成
ここまでで作成した変数を用いて，地価と地価観測地点の属性との関係を検証します．具体的には，地価を被説明変数，地価観測地点の属性を説明変数とした重回帰分析を行います．地価ポイントに含まれる変数のうち，回帰分析で変数に使うのは以下です．標準値コードと供給施設有無（ガス）については，ダミー変数化した上で説明変数に用います．
- L01_001：標準地コード_見出し番号（地価観測地点周辺の土地利用状況を表す変数）
- L01_006：公示価格
- L01_026：地積
- L01_032：供給施設有無（ガス）
- L01_057：容積率
L01_001とL01_032以外の変数を，別の名前の変数として複製します，複製された変数lp（公示地価）とarea（地積）については，対数値の変数も合わせて用意します．
lp_fukuoka_xy_c_near_sta_dist_shop_reg <- lp_fukuoka_xy_c_near_sta_dist_shop %>%
dplyr::mutate(lp=L01_006,
area=L01_026,
far=L01_057) %>%
dplyr::mutate(ln_lp=log(lp),
ln_area=log(area))続いて，変数L01_001に基づくダミー変数を作成します．住宅地（L01_001が”000”）を参照カテゴリとして，それ以外の土地利用カテゴリについてダミー変数を作成します（ダミー変数の数は，最大でも[カテゴリの数-1]個です．詳細は計量経済学の教科書を参照）．地価ポイントの変数L01_001をユニーク化すると，変数L01_001は”000”,“005”,“009”のいずれかの値のみを取ることが分かります．そこで，後者2つのカテゴリについてダミー変数を作成します．
unique(lp_fukuoka_xy_c_near_sta_dist_shop_reg$L01_001)## [1] "000" "005" "009"lp_fukuoka_xy_c_near_sta_dist_shop_reg <- lp_fukuoka_xy_c_near_sta_dist_shop_reg %>%
dplyr::mutate(d_shogyo=ifelse(L01_001=="005",1,0),
d_kogyo=ifelse(L01_001=="009",1,0))最後に，変数L01_032に基づくダミー変数を作成します．もしガス供給施設があれば（L01_032が”true”ならば）1，そうでなければ0を取るダミー変数です．
lp_fukuoka_xy_c_near_sta_dist_shop_reg <- lp_fukuoka_xy_c_near_sta_dist_shop_reg %>%
dplyr::mutate(d_gus=ifelse(L01_032=="true",1,0))モデルの定式化・相関分析
今回推定する重回帰モデルを定式化すると以下のようになります（LaTex形式での数式入力の都合上，作成した変数名と若干違いがあります）．\(\ln lp_{i}\)は地価観測地点\(i\)における対数地価，\(\ln dsta _{i}\)は最寄り駅までの距離（対数），\(\ln nshop_{i}\)は地価観測地点から1km圏内に含まれる店舗の数（対数），\(\ln area_{i}\)は地積（対数値），\(far_{i}\)は容積率，\(dGus_{i}\)はガス供給施設ダミー，\(dShogyo_{i}\)は商業地ダミー，\(dKogyo_{i}\)は工業地ダミー，\(\epsilon_{i}\)は誤差項です．
\[ \ln lp_{i} = \beta_{0} + \ln dsta _{i} \beta_{1} + \ln nshop_{i} \beta_{2} + \ln area_{i} \beta_{3} + far_{i} \beta_{4} + dGus_{i} \beta_{5} + dShogyo_{i} \beta_{6} + dKogyo_{i} \beta_{7} + \epsilon_{i} \tag{1} \]
重回帰分析の事前分析として，変数間の相関係数を見ます．相関係数を行列形式で出力するには，Rのデフォルト関数であるcor関数を用います．一部の説明変数間（例えば，商業地ダミーd_shogyoと容積率far）に高い相関関係が見て取れるため，これら変数について，多重共線性が起きる可能性がある点に注意が必要です．
lp_fukuoka_xy_c_near_sta_dist_shop_reg %>%
#相関係数を計算したい変数のみ残す
dplyr::select(ln_lp,
ln_d_sta,
ln_nshop,
ln_area,
far,
d_gus,
d_shogyo,
d_kogyo) %>%
#geometry属性を削除
sf::st_set_geometry(NULL) %>%
#相関係数行列を作成
cor() %>%
#小数点第三位で四捨五入
round(digits=3)## ln_lp ln_d_sta ln_nshop ln_area far d_gus d_shogyo d_kogyo
## ln_lp 1.000 -0.546 0.620 0.370 0.881 0.340 0.667 -0.113
## ln_d_sta -0.546 1.000 -0.335 -0.134 -0.487 -0.230 -0.367 0.096
## ln_nshop 0.620 -0.335 1.000 0.161 0.550 0.275 0.393 0.016
## ln_area 0.370 -0.134 0.161 1.000 0.351 -0.169 0.258 0.428
## far 0.881 -0.487 0.550 0.351 1.000 0.179 0.769 -0.021
## d_gus 0.340 -0.230 0.275 -0.169 0.179 1.000 0.070 -0.343
## d_shogyo 0.667 -0.367 0.393 0.258 0.769 0.070 1.000 -0.121
## d_kogyo -0.113 0.096 0.016 0.428 -0.021 -0.343 -0.121 1.000重回帰モデルの推定
Rで重回帰分析を行う際には，デフォルト関数のlm関数を用います．引数formulaにはモデルの定式化を指定します．定式化は，
[被説明変数]~[説明変数1]+[説明変数2]+…
という風に行います．引数dataには，分析に用いるデータセットを指定します．推定された回帰係数や，それに対応する標準誤差・t値等は，推定結果のオブジェクトに対してsummary関数を適用すれば得られます．推定結果を解釈すると，
- 対数地価と統計的に有意な正の関係がある：店舗数（対数），地積（対数），容積率，ガス供給施設ダミー
- 対数地価と統計的に有意な負の関係がある：最寄り駅までの距離（対数），工業地ダミー
- 統計的に有意な関係が見られない：商業地ダミー
となり，おおむね直感と整合するものであると思います．商業地ダミーd_shogyoが有意でない理由として考えられるのは，容積率farとの高い相関関係です．modelsummary関数を用いると，より綺麗に推定結果を出力できます．
#回帰モデルを推定
reg <- lm(formula=ln_lp~ln_d_sta+ln_nshop+ln_area+far+d_gus+d_shogyo+d_kogyo,
data=lp_fukuoka_xy_c_near_sta_dist_shop_reg)
#回帰モデルの推定結果をまとめる
summary(reg)##
## Call:
## lm(formula = ln_lp ~ ln_d_sta + ln_nshop + ln_area + far + d_gus +
## d_shogyo + d_kogyo, data = lp_fukuoka_xy_c_near_sta_dist_shop_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24104 -0.25036 0.01053 0.23186 1.20156
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.893547 0.319150 31.000 < 2e-16 ***
## ln_d_sta -0.129703 0.031380 -4.133 4.71e-05 ***
## ln_nshop 0.159594 0.027514 5.800 1.76e-08 ***
## ln_area 0.221694 0.034114 6.499 3.61e-10 ***
## far 0.004767 0.000298 15.998 < 2e-16 ***
## d_gus 0.462953 0.082925 5.583 5.51e-08 ***
## d_shogyo -0.037589 0.082195 -0.457 0.648
## d_kogyo -0.690012 0.154532 -4.465 1.15e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3969 on 285 degrees of freedom
## Multiple R-squared: 0.8591, Adjusted R-squared: 0.8556
## F-statistic: 248.2 on 7 and 285 DF, p-value: < 2.2e-16#回帰モデルの推定結果をまとめる（modelsummary使用）
modelsummary::modelsummary(models=reg,stars=TRUE)| (1) | |
|---|---|
| (Intercept) | 9.894*** |
| (0.319) | |
| ln_d_sta | −0.130*** |
| (0.031) | |
| ln_nshop | 0.160*** |
| (0.028) | |
| ln_area | 0.222*** |
| (0.034) | |
| far | 0.005*** |
| (0.000) | |
| d_gus | 0.463*** |
| (0.083) | |
| d_shogyo | −0.038 |
| (0.082) | |
| d_kogyo | −0.690*** |
| (0.155) | |
| Num.Obs. | 293 |
| R2 | 0.859 |
| R2 Adj. | 0.856 |
| AIC | 299.8 |
| BIC | 332.9 |
| Log.Lik. | −140.911 |
| F | 248.180 |
| RMSE | 0.39 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
発展的な重回帰分析：パネルデータ分析
重回帰分析に関する懸念事項
先のセクションで重回帰分析を用いて明らかにした事実は，例えば「最寄り駅までの距離が小さい時に地価観測地点の地価は高い傾向にある」という相関関係でした．一方この相関関係を，「最寄り駅までの距離が短縮されたら地価が上昇する」という因果関係として解釈することは果たして妥当でしょうか．結論としては，必ずしも妥当とは言えません．そのことを理解する為に，「最寄り駅までの距離が小さい時に地価観測地点の地価は高い傾向にある」ということの背後に起こりうるいくつかのメカニズムを挙げてみましょう．まずは，先述した因果関係
- 最寄り駅までの距離が短縮されたら地価が上昇する
が挙げられると思いますが，そのほかにも以下のようなメカニズムが挙げられると思います．
- 1．駅までのアクセス性と地価に共通して影響する第三の要因により，あたかもアクセス性改善⇒地価上昇なる関係があるように見える
- 駅までのアクセス性が高い地域と地価が高い地域は，共通して地域経済の活性度や経済発展のポテンシャルが高い
- 駅までのアクセス性向上と地価上昇に共通してマクロ経済的要因（例えば景気・投資動向）が影響する
- 駅までのアクセス性が高い地域と地価が高い地域は，共通して地域経済の活性度や経済発展のポテンシャルが高い
- 2．地価上昇が鉄道整備という公共投資需要を喚起する（しない）
- 「地価が上昇傾向にある＝経済発展が進展している」地域と見なされて，優先的に公共投資が行われる
- 地価が高い場所は用地取得に費用がかかるので，地価が安い場所を通るルートが選ばれる
- 「地価が上昇傾向にある＝経済発展が進展している」地域と見なされて，優先的に公共投資が行われる
1のようなメカニズムは交絡，2のようなメカニズムは逆の因果関係と呼ばれるものであり，いずれも本来検証したい因果関係を計測する上での妨げとなります（Figure 3参照）．また，これら2つのメカニズムが引き起こす問題を，まとめて「内生性問題」とも呼びます．
Figure 3：最寄り駅までのアクセス改善と地価上昇の間のメカニズム
回帰分析の枠組みで，交絡因子や逆因果の問題を整理すると次のようになります．今，複数の観測地点・時点で地価が観測されているとしましょう．この時，先ほどの因果関係を検証するために，以下の回帰モデルを推定することを考えます．\(landPrice_{i,t}\)を地点\(i\)で時点\(t\)に観測される（対数）地価，\(access_{i,t}\)を駅までのアクセス性に関する変数，\(control_{i,t}\)をその他地価に影響する地域属性，\(\epsilon_{i,t}\)を誤差項として，
\[ landPrice_{i,t} = \beta_{0} + access_{i,t} \beta_{1} + control_{i,t} + \epsilon_{i,t} \tag{2} \]
必ずしも厳密な説明ではないのですが，この回帰モデルを用いて因果関係を識別するための条件は，アクセス性変数\(access_{i,t}\)と誤差項\(\epsilon_{i,t}\)が無相関，すなわち，
\[ Cov(access_{i,t},\epsilon_{i,t})=0 \tag{3} \]
という条件が成り立つ必要があります．一方，交絡因子や逆因果が存在している際には，この条件が成り立たないことが知られています．その結果，アクセス性変数\(access_{i,t}\)の因果効果を捉える回帰係数\(\beta_{1}\)の推定値はバイアスを持つため，効果計測が妨げられます．
以下では，パネルデータ（複数の個体について，複数の時点で観測値が得られているデータ）を用いて交絡因子の問題に対処するための代表的な計量経済学的手法である，固定効果法と差分の差分法について解説します．逆因果の問題に対処するための方法として操作変数法という手法がありますが，今回の文脈では適用が難しいので説明は省略します．
固定効果法（fixed-effects; FE）
交絡因子の問題が起きる原因のひとつは，本来考慮すべき変数が考慮されないという，除外変数の問題です．この問題を緩和するためのひとつの方法は，無視された変数全てをモデルに含めることですが，データの可用性からほとんどの場合難しいです．固定効果法のモチベーションは，これら除外変数のうち，以下の要因をコントロールすることにあります．
- 個体固定効果：観測個体毎に異なる時間不変の要因
- 例：自然地理的な優位性（海岸からの近さ，自然資源の可用性），都心からの距離，政治・制度的要因
- 時間固定効果：観測個体に共通する時間可変の要因
- 例：マクロ経済変動（景気），全国一律に適用される法制度の変化
個体固定効果を\(\kappa_{i}\)，時間固定効果を\(\rho_{t}\)として，これらを回帰モデルに導入した際の定式化は以下のようになります．
\[ landPrice_{i,t} = \kappa_{i} + \rho_{t} + access_{i,t} \beta_{1} + control_{i,t} + \epsilon_{i,t} \tag{4} \]
なお，個体固定効果を回帰モデルで実装する際の最も単純な方法は，個体\(i\)に該当するレコードで1，他では0を取るようなダミー変数を変数として加えることです．同じく，時間固定効果を実装する際には，時点\(t\)に該当するレコードで1，他では0を取るようなダミー変数を加えます．一方，ダミー変数を用いた推定法は，計算負荷の問題等から必ずしも用いられる訳ではなく，一階差分法や平均差分法と呼ばれる方法を用いて推定されることも多いです（これら推定法の詳細は計量経済学の教科書を参照）．
誤差項と説明変数に正の相関がある場合，すなわち\(Cov(access_{i,t},\epsilon_{i,t})>0\)の場合における，通常の重回帰分析で使われる推定法である最小二乗法（ordinary
least square;
OLS）と固定効果法の違いを図で表すと下の図のようになります．各個体に関して複数時点で観測されたデータを白抜きの青点で示しています．個体情報に関する制御（すなわち個体固定効果の導入）を行わずにOLSで重回帰分析を行った場合，各個体に関する回帰線よりも傾きが急な回帰線が推定されてしまいます．一方，個体固定効果を導入することで，個体間の時間不変な要因の差が制御され，各個体のデータは原点を中心とした場所に集約されます（白抜きでない青点）．集約されたデータに対して回帰線を推定した場合，各個体に関する回帰線と傾きがより近いものを推定することができます．
Figure 4：固定効果法のイメージ（『計量経済学 第2版』の図12-2をもとに筆者作成）
差分の差分法（difference-in-differences; DID）
DIDは，いわゆる自然実験（企業・労働者・場所などの経済主体が適応せざるを得ない環境の変化を実験として利用し，変化の影響を受けるグループを，影響を受けないグループから区別する方法）の手法のひとつです．影響を受けるグループを処置群（treatment
group），影響を受けないグループを対照群（control
group）とし，両群の間でアウトカムを比較することで，環境の変化・ショックの影響を検証できます．
DID推定の必須条件は，処置群・対照群両方のショック前後の状況が比較可能であることです．処置群か対照群かの2つの分類と，処置前か処置後かの2つの分類，合わせて2×2=4つの分類に関するデータがあれば，処置群・対照群間の属性差を制御することが可能（なことが多い）です．そのことをイメージ図で表したのが以下の図です．処置群と対照群について，処置前と処置後のアウトカムが観測されていたとします．DIDの考え方は，対照群のアウトカムの変化から，処置がなかった場合の処置群の仮想的なアウトカムの変化を作り出し，反実仮想的な分析を行うというものです．処置群における処置前後のアウトカムの差分と，対照群における処置前後のアウトカムの差分を計算し，さらにそれらの差分を取るという手続きで処置効果を推定することが，差分の差分法という名前のゆえんです．
Figure 5：DIDのイメージ（『An Introduction to Geographical and Urban Economics』のFigure 3.13を元に筆者作成）
ここで，アクセス性と地価の関係という文脈で，DIDを回帰モデルとして定式化します．\(treat_{i}\)を観測地点\(i\)が処置群（アクセス性が改善されたグループ）に含まれれば1を取るダミー変数，\(after_{t}\)を観測値が処置後（鉄道開業後）に観測されたものであれば1を取るダミー変数とすると，
\[ landPrice_{i,t} = \beta_{0} + treat_{i} \beta_{1} + after_{t} \beta_{2} + treat_{i} \times after_{t} \delta + control_{i,t} + \epsilon_{i,t} \tag{5} \]
この回帰モデルと，DIDのイメージ図（Figure 5）との対応関係を表すと，以下のようになります．
\[ landPrice_{i,t}= \left\{ \begin{array}{l} (A)treat_{i}=0, after_{t}=0:\beta_{0} \\ (B)treat_{i}=1, after_{t}=0:\beta_{0} + \beta_{1} \\ (E)treat_{i}=0, after_{t}=1:\beta_{0} + \beta_{2} \\ (C)treat_{i}=1, after_{t}=1:\beta_{0} + \beta_{1} + \beta_{2} + \delta \end{array} \right.\\ \delta=(C-E)-(B-A) \tag{6} \]
等式(5)での\(treat_{i}\)は，アクセス性が改善した全ての地価観測地点に共通の時間不変要因の平均的傾向を捉える変数です．その代わりに，観測地点レベルの個体固定効果をモデルに導入することによって，地価観測地点毎の時間不変要因の平均的傾向を捉えることができます．\(after_{t}\)も同様に，開業後の全ての時点・全ての地価観測地点に共通する時間可変要因の平均的傾向を捉える変数ですが，その代わりに時点毎に変動する効果を導入することで，開業前後の各時点で異なる時間可変要因の平均的傾向を捉えることができます．
これら個体固定効果・時間固定効果をモデルに導入した回帰モデルは，二方向の「固定効果」を制御するという意味で，two-ways
fixed effect
（TW-FE）の回帰分析と呼ばれます．TW-FEとDID分析を組み合わせることで，より精緻な前後比較が可能となります．個体固定効果を\(\kappa_{i}\)，時間固定効果を\(\rho_{t}\)としてモデルに導入した定式化は，以下のようになります．
\[ landPrice_{i,t} = \kappa_{i} + \rho_{t} + treat_{i} \times after_{t} \beta + control_{i,t} + \epsilon_{i,t} \tag{7} \]
DID推定に必要な別の条件として，平行トレンドの仮定があります．これは処置以前において，処置群と対照群のアウトカムとの「トレンド」が同様でなければならないという仮定です．対照群のデータを用いて処置群のアウトカムの反実仮想的変化を捉えるためには，処置の有無以外の要因が両群間で似通っている必要があります．そうした両群間の同質性があるか否かを確かめる方法（のひとつ）として，トレンドの比較を行うことは，推定された処置効果のバイアスの有無・程度を想定する上で非常に重要です．
平行トレンドの考え方を示したのが下の図です．左は平行トレンドが満たされていそうな場合，右は満たされていなさそうな場合のアウトカムのトレンドです．右のように，平行トレンドが満たされていない場合は，例えば時間可変なコントロール変数の追加や，マッチング法を用いた処置群・対照群間での属性のバランス化（詳しくは計量経済学の教科書を参照）を行い，対処する必要があります．
Figure 6：平行トレンドが満たされていそうな場合（左）・そうでない場合（右）
パネルデータ分析の実例
固定効果法・差分の差分法のようなパネルデータを用いた分析の実例として以下では，単純な定式化の下で，地下鉄開業の効果計測を行う方法を示します．具体的には，2005年に開業した，福岡市営地下鉄七隈線開業前後での地価の変化を検証します．七隈線は路線距離約12kmの地下鉄線で，福岡市の中心地区である天神地区と，福岡市南西地域（城南区・早良区）とを結びます．七隈線開業以前，福岡市南西地域は鉄道空白地帯でしたが，開業によって都心地域へのアクセス性が改善されました．更に近年では，天神地区から博多駅までの延伸工事が進められており，さらなる利便性向上が見込まれている路線です．
七隈線の駅の位置情報を把握するため，駅ポイントを地図上にプロットしてみます．そのためにまず，先ほど鉄道時系列データから作成した，現存する鉄道駅のポイントデータから，七隈線の駅ポイントを取り出します．変数N05_002（路線名）が「3号線」で，変数N05_003（運営会社）が「福岡市」であるレコードが，七隈線の駅ポイントに該当します．次に，取り出したポイントデータを，mapview関数でプロットします．
rail_jp_xy_exist_nanakuma <- rail_jp_xy_exist %>%
#七隈線に該当するポイントの抽出
dplyr::filter(N05_002=="3号線"&N05_003=="福岡市")
#駅ポイントのプロット
mapview::mapview(x=rail_jp_xy_exist_nanakuma)今回の分析では，処置群を七隈線の駅から1km圏内にある福岡市内の地価観測地点群，そうでない市内観測地点群を対照群と定義します．そのため，地価パネルデータの構築に先立ち，駅から1km圏内のバッファを作成しておきます．バッファを作成した後，後の分析のため，全てのレコードで1を取る変数treatを追加します．こちらも，駅ポイントと合わせて位置情報を確認しておきましょう．
rail_jp_xy_exist_nanakuma_buffer <- rail_jp_xy_exist_nanakuma %>%
#半径1000mのバッファを作成
sf::st_buffer(dist=1000) %>%
#全てのレコードで1を取る変数treatを作成
dplyr::mutate(treat=1) %>%
#変数treatだけを残す
dplyr::select(treat)
#バッファのプロット
mapview::mapview(x=rail_jp_xy_exist_nanakuma_buffer) +
#駅ポイントのプロット
mapview::mapview(x=rail_jp_xy_exist_nanakuma,col.regions="red")今回の分析の対象期間は，2002〜2008年と設定します．地価ポイントデータには，過去の地価も変数として含まれていますので，2002〜2008年の観測値を抽出した上で，パネルデータを作成します．1983〜2022年までの地価は，変数L01_061〜L01_100に格納されていますので，地価ポイントからgeometry属性を落とした上で，変数L01_061〜L01_100と各種必要な変数だけを残した新たなデータフレームを作成します．最後に，後の分析で扱いやすいように，変数名を変更します．
lp_fukuoka_xy_c_hist <- lp_fukuoka_xy_c %>%
#geometry属性を削除
sf::st_set_geometry(NULL) %>%
#必要な変数のみ残す
dplyr::select(LID,L01_024,CSS_NAME,L01_061:L01_100)
#変数名を変更
colnames(lp_fukuoka_xy_c_hist) <- c("LID","L01_024","CSS_NAME",1983:2022)
#先頭行を表示
head(lp_fukuoka_xy_c_hist)## LID L01_024 CSS_NAME 1983 1984
## 336 336 福岡県 福岡市中央区赤坂２丁目２区５４番 中央区 0 0
## 337 337 福岡県 福岡市中央区大濠１丁目１３区１３３番 中央区 252000 267000
## 338 338 福岡県 福岡市中央区警固３丁目１０３番 中央区 0 0
## 339 339 福岡県 福岡市中央区荒戸３丁目２９７番１ 中央区 0 0
## 340 340 福岡県 福岡市中央区六本松４丁目６１番 中央区 0 0
## 341 341 福岡県 福岡市中央区谷１丁目２９３番 中央区 0 0
## 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994
## 336 0 0 0 0 0 0 0 0 0 0
## 337 273000 277000 305000 460000 550000 650000 750000 750000 650000 600000
## 338 0 0 0 0 0 0 0 0 0 0
## 339 0 0 0 0 0 450000 519000 519000 455000 385000
## 340 0 0 0 0 270000 300000 365000 365000 328000 315000
## 341 0 0 0 0 0 0 0 0 0 0
## 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004
## 336 0 0 0 0 0 0 260000 257000 253000 251000
## 337 555000 520000 475000 443000 423000 410000 400000 393000 393000 393000
## 338 0 0 0 0 0 0 0 0 0 0
## 339 340000 308000 285000 270000 260000 248000 240000 228000 220000 214000
## 340 275000 245000 227000 216000 212000 207000 203000 195000 187000 181000
## 341 0 0 0 0 0 0 0 0 0 0
## 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014
## 336 252000 259000 300000 331000 308000 293000 291000 292000 305000 321000
## 337 400000 420000 525000 603000 535000 494000 494000 501000 518000 534000
## 338 0 0 0 0 0 0 0 0 0 178000
## 339 214000 221000 284000 328000 310000 290000 281000 282000 289000 299000
## 340 181000 182000 214000 235000 223000 215000 208000 213000 217000 229000
## 341 0 0 0 0 0 0 0 0 159000 164000
## 2015 2016 2017 2018 2019 2020 2021 2022
## 336 340000 359000 391000 425000 463000 504000 528000 575000
## 337 566000 605000 650000 698000 756000 800000 825000 870000
## 338 186000 197000 208000 219000 232000 266000 289000 324000
## 339 311000 329000 348000 373000 389000 438000 458000 488000
## 340 242000 262000 288000 316000 346000 394000 428000 461000
## 341 173000 187000 201000 212000 227000 249000 266000 300000現状のデータは，観測地点について，行方向に過去の地価が変数として並んでいる状態です（wide型という構造です）．このままでは，後のパネルデータ分析で扱えないので，long型（1行につき1つの観測値が格納されたデータ型）に変換します．まず，2002〜2008年の地価観測値の列を取り出し，全ての時点で地価が観測されている観測地点のみを抽出します．次に，tidyrパッケージ（tidyverseパッケージを構成するパッケージ群のひとつ）のpivot_longer関数を用いて，wide型からlong型へ変換を行います．引数dataはlong型に変換したいデータ，colsはlong型に変換したい列の名前（変数を-c()で囲うと，囲った変数「以外の」変数がlong型に変換される），names_toはlong型に変換される変数の元の名前を新たな列として加える際の列名，values_toはlong型に変換される変数の列名を指定します．long型に変換された後の地価観測値は変数lpに格納されますが，この対数値も新たな変数として追加します．また，後の分析のために，時点変数yearを数値型に変換します．
lp_fukuoka_xy_c_hist <- lp_fukuoka_xy_c_hist %>%
#変数LID,L01_024,CSS_NAMEと，過去の地価観測値の変数（2002から2008）のみ選択
dplyr::select(LID,L01_024,CSS_NAME,c("2002":"2008")) %>%
#変数2002から2008が全ての年で正のレコードのみ抽出
dplyr::filter(across(c("2002":"2008"),~.>0))
#wide型からlong型への変換
lp_fukuoka_xy_c_hist_long <- tidyr::pivot_longer(data=lp_fukuoka_xy_c_hist,
#LID,L01_024,CSS_NAME以外の変数をlong型に変換
#この場合，LID,L01_024,CSS_NAMEは時点ごとに繰り返される
cols=-c(LID,L01_024,CSS_NAME),
#元の変数名を格納する変数の名前
names_to="year",
#long型に変換された変数を格納する列の名前
values_to="lp") %>%
#地価lpの対数値を変数として追加
dplyr::mutate(ln_lp=log(lp)) %>%
#変数yearを数値型に変換
dplyr::mutate(year=as.numeric(year))
head(lp_fukuoka_xy_c_hist_long)## # A tibble: 6 × 6
## LID L01_024 CSS_NAME year lp ln_lp
## <int> <chr> <chr> <dbl> <int> <dbl>
## 1 336 福岡県 福岡市中央区赤坂２丁目２区５４番 中央区 2002 257000 12.5
## 2 336 福岡県 福岡市中央区赤坂２丁目２区５４番 中央区 2003 253000 12.4
## 3 336 福岡県 福岡市中央区赤坂２丁目２区５４番 中央区 2004 251000 12.4
## 4 336 福岡県 福岡市中央区赤坂２丁目２区５４番 中央区 2005 252000 12.4
## 5 336 福岡県 福岡市中央区赤坂２丁目２区５４番 中央区 2006 259000 12.5
## 6 336 福岡県 福岡市中央区赤坂２丁目２区５４番 中央区 2007 300000 12.6通し番号変数LIDのみを残した観測地点ポイントを作成し，上で作成した七隈線駅1kmバッファを空間結合します．その後，バッファが結合されなかった観測地点ポイントの変数treatを0に置換し，通し番号変数LIDに関してユニークにします（ユニークにする理由は，ひとつの観測地点が複数の駅バッファに属す際に変数LIDに重複が起きるためです）．geometry属性は以後使いませんので，削除します．
lp_fukuoka_xy_c_point <- lp_fukuoka_xy_c %>%
#通し番号変数LIDのみ選択
dplyr::select(LID)
#観測地点ポイントにバッファを空間結合
lp_fukuoka_xy_c_point_nanakuma_buffer <- sf::st_join(x=lp_fukuoka_xy_c_point,
y=rail_jp_xy_exist_nanakuma_buffer,
join=st_intersects) %>%
#バッファが結合されなかった観測地点ポイントの変数treatを0に置換
dplyr::mutate(treat=ifelse(!is.na(treat),treat,0)) %>%
#LIDについてユニークにする
dplyr::distinct(LID,.keep_all=TRUE) %>%
#geometry属性を削除
sf::st_set_geometry(value=NULL)long型の地価パネルデータに対して，空間結合の結果を左結合します．時点変数yearが2005以上の場合1を取るダミー変数afterも合わせて作成します．
lp_fukuoka_xy_c_hist_long_nanakuma <- lp_fukuoka_xy_c_hist_long %>%
#long型の地価パネルデータに空間結合結果を結合
dplyr::left_join(lp_fukuoka_xy_c_point_nanakuma_buffer,by="LID") %>%
#時点変数yearが2005以上の場合1を取るダミー変数を作成
dplyr::mutate(after=ifelse(year>=2005,1,0))DIDを推定するに先立ち，いくつかの記述統計を確認しておきます．まず，処置有無treat・時点yearでデータをグループ化した上で，各グループのレコード数，対数地価の平均値，標準偏差，標準誤差を計算します．また，後ほど行うプロットのため，集計結果内の変数treatを文字列型に変換します．
lp_fukuoka_xy_c_hist_long_nanakuma_sum <- lp_fukuoka_xy_c_hist_long_nanakuma %>%
#データを処置有無treat・時点yearでグループ化
dplyr::group_by(treat,year) %>%
#レコード数を計算
dplyr::summarise(obs=n(),
#平均値を計算
ln_lp_mean=mean(ln_lp),
#標準偏差を計算
ln_lp_sd=sd(ln_lp),
#標準誤差を計算
ln_lp_se=sd(ln_lp)/sqrt(obs)) %>%
#グループ化情報を削除
dplyr::ungroup() %>%
#処置変数treatを文字列型に変換
dplyr::mutate(treat=as.character(treat))
lp_fukuoka_xy_c_hist_long_nanakuma_sum## # A tibble: 14 × 6
## treat year obs ln_lp_mean ln_lp_sd ln_lp_se
## <chr> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0 2002 158 12.0 0.638 0.0507
## 2 0 2003 158 11.9 0.634 0.0505
## 3 0 2004 158 11.9 0.635 0.0505
## 4 0 2005 158 11.8 0.638 0.0508
## 5 0 2006 158 11.8 0.656 0.0522
## 6 0 2007 158 11.8 0.724 0.0576
## 7 0 2008 158 11.8 0.787 0.0626
## 8 1 2002 32 12.5 0.884 0.156
## 9 1 2003 32 12.4 0.891 0.158
## 10 1 2004 32 12.4 0.911 0.161
## 11 1 2005 32 12.3 0.923 0.163
## 12 1 2006 32 12.3 0.952 0.168
## 13 1 2007 32 12.4 1.04 0.184
## 14 1 2008 32 12.5 1.14 0.201集計結果をグラフとしてプロットします．処置群・対照群別に，年ごとの対数地価の平均と標準誤差を線グラフにして可視化します．開通年である2005年以前での処置群・対照群間でのアウトカムの傾向を見てみると，同程度に下降傾向になっているので，平行トレンドの仮定は満たされていそうな感じがします．2005年以降，両群で地価が上昇トレンドに転じていますが，上昇度は処置群の方でより大きいことが分かります．
#集計データを元にプロットするという宣言（x軸に変数year，y軸に平均値ln_lp_mean，処置変数treatでグラフをグループ化，treatの値毎に色分け）
ggplot2::ggplot(data=lp_fukuoka_xy_c_hist_long_nanakuma_sum,aes(x=year,y=ln_lp_mean,group=treat,color=treat)) +
#線グラフをプロット
ggplot2::geom_line() +
#線グラフ上にyear毎の点をプロット（点のサイズは3）
ggplot2::geom_point(size=3) +
#標準誤差はエラーバーでプロット
ggplot2::geom_errorbar(aes(ymin=ln_lp_mean-ln_lp_se,ymax=ln_lp_mean+ln_lp_se),width=0.1,color="black") +
#x=2005で垂直線を点線で引く
ggplot2::geom_vline(xintercept=2005,linetype="dashed")
等式(7)で定式化した固定効果DIDを推定します．固定効果モデルの推定には，fixestのfeols関数を用います．引数fmlでは，「|」よりも左側で定式化を，右側で含める固定効果（今回は変数LIDに基づく個体固定効果と変数yearに基づく時点固定効果の2つ）を指定します．引数dataにはデータセットを指定します．推定結果を表示するには，fixestのetable関数を用います．最初に推定したモデルのオブジェクトを，引数signifCodeでは有意水準とそれに対応した記号を，digitsでは回帰係数を小数点何桁までを表示するかを，digits.statsではt値等の各種統計量を小数点何桁までを表示するかを指定します．また，modelsummary関数による推定結果の表示も合わせて行います．
推定の結果，処置効果に関する変数treatxafterの回帰係数の推定値は，正かつ1%水準で統計的に有意です．故に，七隈線の開業が，駅周辺地域の地価を上昇させたことを示唆する結果が得られました．
#固定効果DIDの推定
base_fe <- fixest::feols(fml=ln_lp~treat:after|LID+year,data=lp_fukuoka_xy_c_hist_long_nanakuma)
#推定結果の表示（etable使用）
fixest::etable(base_fe,signifCode=c("***"=0.01,"**"=0.05,"*"=0.10),digits=3,digits.stats=3)## base_fe
## Dependent Var.: ln_lp
##
## treat x after 0.110*** (0.027)
## Fixed-Effects: ----------------
## LID Yes
## year Yes
## _______________ ________________
## S.E.: Clustered by: LID
## Observations 1,330
## R2 0.989
## Within R2 0.059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#推定結果の表示（modelsummary使用）
modelsummary::modelsummary(models=base_fe,stars=TRUE)| (1) | |
|---|---|
| treat × after | 0.110*** |
| (0.027) | |
| Num.Obs. | 1330 |
| R2 | 0.989 |
| R2 Adj. | 0.987 |
| R2 Within | 0.059 |
| R2 Within Adj. | 0.059 |
| AIC | −2517.8 |
| BIC | −1494.8 |
| RMSE | 0.08 |
| Std.Errors | by: LID |
| FE: LID | X |
| FE: year | X |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
ここで，今回の分析で得られた結果が頑健（robust）であることを確かめる際には，例えば以下のような点に対して更なる検証を行うことが重要です．
- 今回の分析では，1kmバッファという直線距離に依拠したアクセス性指標に基づき，処置群を定義しました．一方，別のバッファ距離（例えば1.5km）を使った場合でも，結果が大きく変わらないかが確認できると，結果はより説得的になります．また，駅周辺の道路ネットワーク等のデータが利用可能な場合，ネットワーク距離や所要時間など，より精緻なアクセス性の変数に基づき分析することが望ましいです．
- DID分析の前提条件のひとつである，「介入（新線開業）の影響は，介入を受けなかった集団（開業駅から1km以上離れた場所）には及ばない」というStable
Unit Treatment Value
Assumption（SUTVA）が満たされているかは，今回の分析デザインでは疑わしいです．新線開業効果が1kmという範囲を超えてスピルオーバーしたりすると，SUTVAは破れてしまう為，それに対処した分析を追加的に行う必要が出てきます．
- SUTVAの破れを検証する最もシンプルな方法として例えば，開業駅1km圏外だが，それと地理的に近い対照群のサンプルについて「偽の」処置変数を作成・モデルに変数として加えた上で，それが統計的に有意になるかどうか見るというものが挙げられます．仮に有意な結果が得られた場合，それら対照群のサンプルを外した上でも，結果が変わらないかを検証することが望ましいです．
- SUTVAの破れを検証する最もシンプルな方法として例えば，開業駅1km圏外だが，それと地理的に近い対照群のサンプルについて「偽の」処置変数を作成・モデルに変数として加えた上で，それが統計的に有意になるかどうか見るというものが挙げられます．仮に有意な結果が得られた場合，それら対照群のサンプルを外した上でも，結果が変わらないかを検証することが望ましいです．
- 今回のTW-FEのDID分析では，地価観測地点毎に異なる時間可変の地域要因については制御を行っていません．故に，それら地域要因に関する変数をコントロール変数として回帰モデルに追加しても結果が変わらないか，更に分析を行う必要があります．
- 駅の場所や地下鉄のルートが，各場所が持つ経済的属性（例えば，潜在的な開発需要）に応じて決定されている可能性があります．その場合，処置変数\(treat_{i} \times after_{t}\)と誤差項\(\epsilon_{i,t}\)の間に相関が生じてしまい，推定される介入効果が過大もしくは過小になります．これについても，例えば操作変数法と組み合わせたDID分析を行う等して，追加的な対処が必要です．
補足
任意の2つの地物間を結ぶラインの生成
本編では，st_join関数を用いて，各地価ポイントから最寄りの駅ポイントの属性を結合できることを示しました．また，st_connect関数を使うことで，k近傍の空間データとの間にラインを生成できることも，同じく示しました．ここでは最（k）近傍に限らず，任意の2つの地物間にラインを生成する方法を示します．
自作関数を定義
sfパッケージには，任意の2つの地物間のラインをdata.frameから一括して生成する関数は実装されていないので，いくつかの関数を組み合わせることで，この処理を行う関数を自作します．ラインデータの生成に先立ち，最寄りの駅ポイントの属性が付与された地価ポイントから，geometry属性を削除したデータフレームを新たに作成します．
#最寄りの駅ポイントの属性が付与された地価ポイントからgeometry属性を削除する
lp_fukuoka_xy_c_near_sta_dist_df <- lp_fukuoka_xy_c_near_sta_dist %>%
sf::st_set_geometry(value=NULL)関数のコンセプト
簡単な例として，上のデータフレームの1番目のレコードについて，ラインデータを生成する手順を示します．まず，1つ目のレコード内の地価ポイント・最寄り駅ポイントの座標値を取り出します．unlist関数を用いて取り出した座標値を数値ベクトルに変換したのち，行列形式に変形します．変形する際には，1行目に1つ目のポイントのXY座標，2行目に2つ目のポイントのXY座標が来るようにします．
coords_pair <- lp_fukuoka_xy_c_near_sta_dist_df[1, c("X_lp","Y_lp","X_sta","Y_sta")]
coords_pair## X_lp Y_lp X_sta Y_sta
## 336 -56754.69 65075.79 -56571.61 65488.53coords_pair_mat <- matrix(data=unlist(coords_pair),nrow=2,ncol=2,byrow=TRUE)
coords_pair_mat## [,1] [,2]
## [1,] -56754.69 65075.79
## [2,] -56571.61 65488.53この行列からラインデータを作成する際には，sfのst_linestring関数を用います．この関数では，与えられた2つの座標を結ぶsfg形式のラインデータを生成することができます．引数xには，座標値が入った行列を指定します．ラインデータの中身を見てみると，WKT形式で記述されたgeometry属性が確認できます．
coords_pair_line <- sf::st_linestring(x=coords_pair_mat) %>%
sf::st_sfc()
coords_pair_line## Geometry set for 1 feature
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -56754.69 ymin: 65075.79 xmax: -56571.61 ymax: 65488.53
## CRS: NA関数の実装
上で示した一連の流れを，地価ポイントの各レコードに並行して適用できるように関数st_xy_to_lineとして実装すると以下のようになります．その際，後の処理で扱いやすいよう，sfg形式からsfc形式にラインデータを変換する処理を組み込んでいます．
#引数を[始点のX座標],[始点のY座標],[終点のX座標],[終点のY座標]とする関数を定義
st_xy_to_line <- function(x_orig,y_orig,x_dest,y_dest){
#xがdata.frameの行である場合はベクトル化した上で（要素数は4），2×2の行列に変換する
line <- matrix(data=unlist(c(x_orig,y_orig,x_dest,y_dest)),nrow=2,ncol=2,byrow=TRUE) %>%
#行列からラインデータに変換
sf::st_linestring() %>%
#sfc形式に変換
sf::st_sfc()
#結果を返す
return(line)
}関数の適用
自作関数をmutate関数と組み合わせ，地価ポイントの各レコードに対して一括で適用します．そのためにまず，data.frame形式の地価ポイントを，tibbleという形式（data.frameの進化版で，list形式のオブジェクトを変数として保持できる等の特徴がある）に変換します．sfc形式の空間データの変数型は，R上ではlistとして扱われるため，変数として追加するに先立ちtibble形式に変換する必要があります．続いて，dplyrのrowwise関数を適用し，行ごとの関数適用を可能な形式にデータを変換します．
これら前準備を経て，ようやくmutate関数の中で自作関数st_xy_to_lineを適用し，地価観測地点の座標を始点座標，最寄駅の座標を終点座標としたラインデータが作成されます．最後に，sfのst_sf関数を用いて，現状全体としてはtibble形式として認識されているデータを，sf形式として認識させます．その際，引数sf_column_nameにはgeometry情報が入った変数，crsには座標系を指定します．
lp_fukuoka_xy_c_near_sta_dist_line <- lp_fukuoka_xy_c_near_sta_dist_df %>%
#tibbleに変換
dplyr::as_tibble() %>%
#行処理が可能な形に変換
dplyr::rowwise() %>%
#データ内の変数X_lp,Y_lp,X_sta,Y_staを引数とした関数st_xy_to_line（上で定義されたもの）の実行結果を，変数geom_lineとして追加
dplyr::mutate(geom_line=st_xy_to_line(x_orig=X_lp,
y_orig=Y_lp,
x_dest=X_sta,
y_dest=Y_sta)) %>%
#変数geom_lineはsfc形式になっているので，sf形式に変換する
sf::st_sf(sf_column_name="geom_line",crs=sf::st_crs(2444)) %>%
#グループ化を解除
dplyr::ungroup()
head(lp_fukuoka_xy_c_near_sta_dist_line)## Simple feature collection with 6 features and 160 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -58480.76 ymin: 63836.13 xmax: -56564 ymax: 66100.34
## Projected CRS: JGD2000 / Japan Plane Rectangular CS II
## # A tibble: 6 × 161
## L01_001 L01_002 L01_003 L01_004 L01_005 L01_006 L01_007 L01_008 L01_009
## <chr> <chr> <chr> <chr> <chr> <int> <dbl> <int> <chr>
## 1 000 001 000 001 2022 575000 8.9 1 false
## 2 000 002 000 002 2022 870000 5.5 1 false
## 3 000 003 000 003 2022 324000 12.1 1 false
## 4 000 004 000 004 2022 488000 6.6 1 false
## 5 000 005 000 005 2022 461000 7.7 1 false
## 6 000 006 000 006 2022 300000 12.8 1 false
## # … with 152 more variables: L01_010 <chr>, L01_011 <chr>, L01_012 <chr>,
## # L01_013 <chr>, L01_014 <chr>, L01_015 <chr>, L01_016 <chr>, L01_017 <chr>,
## # L01_018 <chr>, L01_019 <chr>, L01_020 <chr>, L01_021 <chr>, L01_022 <chr>,
## # L01_023 <chr>, L01_024 <chr>, L01_025 <chr>, L01_026 <int>, L01_027 <chr>,
## # L01_028 <chr>, L01_029 <chr>, L01_030 <chr>, L01_031 <chr>, L01_032 <chr>,
## # L01_033 <chr>, L01_034 <chr>, L01_035 <dbl>, L01_036 <dbl>, L01_037 <int>,
## # L01_038 <int>, L01_039 <chr>, L01_040 <chr>, L01_041 <dbl>, …最後に，地価ポイントデータ・駅ポイントデータとともに，作成した地価ポイント-駅ラインデータを可視化します．それに先立ち，全国をカバーした駅ポイントから，最寄り駅として地価ポイントに空間結合されたレコードのみを取り出します．
rail_jp_xy_exist_coords_nearest <- rail_jp_xy_exist_coords %>%
#最寄り駅情報結合済みの地価ポイント内にある，通し番号変数SID（駅名）と値が一致するレコードのみ保持
dplyr::filter(SID%in%unique(lp_fukuoka_xy_c_near_sta_dist$SID))mapview関数を用いて，地価ポイントデータ・駅ポイントデータ・地価ポイント-駅ラインデータを可視化します．
#地価ポイントデータ
mapview::mapview(x=lp_fukuoka_xy_c_near_sta_dist,
#ポイントの色は紫
col.region="purple",
#不透明度100%
alpha.regions=1,
#レイヤ名を「Obs Point」と設定
layer.name="Obs Point") +
#駅ポイントデータ
mapview::mapview(x=rail_jp_xy_exist_coords_nearest,
#ポイントの色は緑
col.region="green",
#不透明度100%
alpha.regions=1,
#レイヤ名を「Station」と設定
layer.name="Station") +
#地価ポイント-最寄駅間ラインデータ
mapview::mapview(x=lp_fukuoka_xy_c_near_sta_dist_line,
#ラインの色は赤
color="red",
#レイヤ名を「Station-Obs Point」と設定
layer.name="Station-Obs Point")