Mimicking ArcGIS Tools on R
Keisuke Takano
4/20/2020
Introduction
このWEBサイトに設置したリンクは，新しいタブ・ウィンドウで開くことを推奨します．バグの影響で，新しいタブ・ウィンドウでないと開けないリンクが多いためです．
Please open external links on this page in new tab or new window. Due to unknown bugs, many of links do not open in present window.
What We Do
このページでは，ArcGISを用いてしばしば行われるジオプロセシングの一部を，Rのsfパッケージ及びdplyrパッケージを用いて実装する方法を紹介します．
なお，dplyrの使い方について，WEBサイトであれば，例えばmatsuou1さんのdplyrを使いこなす！基礎編及びdplyrを使いこなす！JOIN編の内容を，書籍であれば，『RユーザのためのRStudio［実践］入門』の3・4章の内容を理解していることを前提として話を進めます．
また，このWEB教材の内容は，RStudioの公式サイトにあるsfとdplyrのCheatsheetsや，Sho Kuroda先生のRPubsRで行う公示地価パネルデータの整備と分析を下敷きとしています（cheetsheetのcheetsheetということになります，はい）．
今回は，下のアメリカ及びベトナムの空間情報データを用いて，Rでのジオプロセシングの実践方法を紹介していきます．なお，このWEBページの元となっているRmdファイルについては，筆者のWEBサイトに設置したこちらのRmdをダウンロードしてください．
＜アメリカ（ニューヨーク市）のデータ＞
- 地下鉄駅のポイントデータ
- 国勢調査区境界のポリゴンデータ
- 土地評価価格データのcsvデータ
＜ベトナムのデータ＞
- COVID-19感染状況が紐づいた行政区域ポリゴンデータ
- 夜間光量データ（人口衛星から撮影される夜間光のラスタデータ）
- 標高データ（標高のラスタデータ）
This page provides how to replicate the geoprocessing procedures in ArcGIS with R packages like sf and dplyr. Main readership of this page is those who already understand and handle dplyr package and basic concepts of GIS. For readers who have not learned dplyr, please read and try the contents of, for example, STAT 545 by Professor Jenny Bryan.
In making this page, I got inspirations from Cheetsheets of sf and dplyr on official WEB site of RStudio and RPubs by Professor Sho Kuroda.
I introduce the procedures with R utilizing several geographical datasets in the US and Vietnam below. Original Rmd file (click here) is distributed on author’s WEB page.
Basis of GIS
GISに関する各種概念の詳細な解説はこのページの中では行いませんが，ESRIジャパン株式会社さんのGIS 基礎解説より，このサイトで紹介されている方法を実践するにあたって事前に知っておくべき基礎概念を引用します．
以下に挙げる順に，各ページの内容に目をお通しください．
- GIS データモデルの基礎
- 座標系/空間参照
- EsriのGISデータフォーマット
- ArcGISと共通するGISの基礎知識
- 空間解析
- その他 データ/データモデル関連
その他，無償で利用可能なGISのチュートリアルとしては，QGISの利用を前提として書かれたものですが，科学研究費補助金基盤研究(A)「GISの標準コアカリキュラムと知識体系を踏まえた実習用オープン教材の開発」で開発された，GIS実習オープン教材が非常に有用です．
On this page, I will not describe the concepts of GIS in detail. Instead, the fundamental knowledge of GIS essential to learn the contents on this page is explained well on Intro to GIS and Spatial Analysis by Professor Manuel Gimond, especially from Chapter 1 to 9.
Concordance to ArcGIS
これから紹介するRでのジオプロセシングは，ArcGISのツールと粗方次のように対応づけることができます．
Geoprocessing procedures with R explained in this page are roughly corresponding to following ArcGIS tools.
- Define Projection：投影法の定義/Define Projection
- Transform Projection：投影変換/Project
- Dissolve：ディゾルブ/Dissolve
- Clip：クリップ/Clip
- Buffer：バッファ/Buffer,インターセクト/Intersect
- Spatial Join：空間結合/Spatial Join or テーブルの結合とリレート/Joining and Relating Tables
- XY to Point：X、Y 座標データをレイヤーとして追加する/Add x,y Coordinate Data as a Layer
- Geometry Calculation：ジオメトリ演算/Calculate Geometry
- Table Join：テーブルの結合とリレート/Joining and Relating Tables
- Summarizing Rasters with Polygons：Extraction
- Elevation：傾斜角/Slope,コスト距離解析/Cost Distance
Further Reading
sfパッケージとdplyrパッケージを組み合わせたGISに関して，さらに詳細・高度なジオプロセシングを解説しているWEBサイトをいくつか紹介します．日本語のサイトとしては，以下がおすすめです（作成者名の敬称略）．
- Rを使った地理空間データの可視化と分析 (S. Uryu)
These websites provide more advanced and detailed geoprocessing procedures with R.
- 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)
Updates
- 説明部分のタイポを修正しました（2020年4月26日）．
- ラスタデータ分析に関する内容を追加しました（2020年6月24日）．
- 筆者WEBページの移動に伴う各種リンクのURL変更に対応しました．また，外部チュートリアルサイトのまとめを更新しました（2020年8月7日）．
Setting
使うパッケージを読み込みます．なお，メインではないですが，以下のパッケージも読み込んでおきます．
インストールが済んでいない場合は，予めインストールを済ませておいてください．
- curl：各種ファイルをダウンロードするときに用いる
- ggplot2：読み込み・作成したshapefileをプロットするときに用いる
- units：距離単位（m,kmなど）を扱うときに用いる
- ggspatial：地図のプロットに方位縮尺記号を描画する際に用いる
- data.table：大型のcsvファイルを読み書きするときに用いる
- RColorBrewer：地図をプロットする際のカラーパレットを作る時に用いる
- geojsonsf：GeoJSONファイルをsfパッケージに読み込む時に用いる
- R.utils：特殊な圧縮形式のファイルを解凍するときに用いる
- lwgeom：不正な地物を修正するときに用いる
- movecost：コスト距離解析に用いる
- raster：ラスタデータのハンドリング・解析に用いる
- magrittr：パイプ処理内の各種操作に用いる
なお，作業ディレクトリは，このRmdファイルを置いた場所になるみたいです．
We will load packages below. If you have not finished installing them, please install in advance. Other than dplyr and sf, we also load these packages:
- curl: package for downloading files from website
- ggplot2: package for plotting imported or generalted shapefiles and other graphs
- units: package for handling variables with units [e.g., m^2]
- ggspatial: package for adding additional objects like scale bar and spatial-aware north arrow
- data.table: package for importing and exporting large csv file
- RColorBrewer: package for generating color palette to plot maps
- geojsonsf: package for importing GeoJSON file
- R.utils: package for handling special compressed file
- lwgeom: package for making an invalid geometry valid
- movecost: package for searching the least-cost path
- raster: package for handling and analyzing rastar data
- magrittr: package for various pipe operations
Working directory seems automatically set to the directory where we put Rmd file.
library(dplyr)
library(sf)
library(curl)
library(ggplot2)
library(units)
library(ggspatial)
library(data.table)
library(RColorBrewer)
library(geojsonsf)
library(R.utils)
library(lwgeom)
library(movecost)
library(raster)
library(magrittr)Download
分析用いる地下鉄駅・国勢調査区境界のshapefileを読み込みます
We use following shapefiles. Since land price dataset is available in csv format, we download and import it in XY to Point section.
##Subway Stations
#ダウンロードしたzipファイルを格納する一時ディレクトリを設定（set temporary directory to put downloaded zip file）
sta <- tempfile()
#zipファイルをダウンロード（download zip file）
curl::curl_download(url="https://data.cityofnewyork.us/api/geospatial/arq3-7z49?method=export&format=Shapefile",destfile=sta)
#zipを解凍するためのサブフォルダを作成（create subfolder to extract shapefile）
dir.create("sta_nyc")
#サブフォルダにzipファイルを解凍（extract shapefile from zip file to subfolder）
unzip(zipfile=sta,exdir=paste0(getwd(),"/sta_nyc"))
#一時ディレクトリを削除（remove temporary directory）
unlink(sta);rm(sta)
##Boundaries
boundary <- tempfile()
curl::curl_download(url="https://data.cityofnewyork.us/api/geospatial/fxpq-c8ku?method=export&format=Shapefile",destfile=boundary)
dir.create("boundary_nyc")
unzip(zipfile=boundary,exdir=paste0(getwd(),"/boundary_nyc"))
unlink(boundary);rm(boundary)Import Shapefiles
ダウンロード＆解凍したshapefileをst_read関数で読み込みます．
We will load extracted shapefiles with st_read function.
##Subway Stations
#shapefileを展開したフォルダに移動（move to subfolder where shapefile was extracted）
setwd("sta_nyc")
#展開先フォルダにあるファイルのリスト（create a list of files in subfolder）
list_files_sta <- list.files()
#リストから拡張子がshpのファイルの名前を取得（obtain name of file with shp extention）
shp_nam_sta <- list_files_sta[grepl(".shp",list_files_sta)]
#shapefileを読み込む．dsnにはshapefileの名前を指定する（import shapefile, file name is designated by dsn）
#読み込むshapefileの属性の中で文字列型のものをfactor型として扱わないので，stringsAsFactors=Fとする（import each character variable as is, not as factor variable）
sta_nyc <- sf::st_read(dsn=shp_nam_sta,stringsAsFactors=F,quiet=TRUE)
#shapefileの属性の先頭行を表示（print first 6 rows of attribute table）
head(sta_nyc)## Simple feature collection with 6 features and 5 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -74.00019 ymin: 40.66471 xmax: -73.89489 ymax: 40.88467
## epsg (SRID): 4326
## proj4string: +proj=longlat +ellps=WGS84 +no_defs
## line name
## 1 4-6-6 Express Astor Pl
## 2 4-6-6 Express Canal St
## 3 1-2 50th St
## 4 2-3-4 Bergen St
## 5 3-4 Pennsylvania Ave
## 6 1 238th St
## notes
## 1 4 nights, 6-all times, 6 Express-weekdays AM southbound, PM northbound
## 2 4 nights, 6-all times, 6 Express-weekdays AM southbound, PM northbound
## 3 1-all times, 2-nights
## 4 4-nights, 3-all other times, 2-all times
## 5 4-nights, 3-all other times
## 6 1-all times, exit only northbound
## objectid url geometry
## 1 1 http://web.mta.info/nyct/service/ POINT (-73.99107 40.73005)
## 2 2 http://web.mta.info/nyct/service/ POINT (-74.00019 40.7188)
## 3 3 http://web.mta.info/nyct/service/ POINT (-73.98385 40.76173)
## 4 4 http://web.mta.info/nyct/service/ POINT (-73.975 40.68086)
## 5 5 http://web.mta.info/nyct/service/ POINT (-73.89489 40.66471)
## 6 6 http://web.mta.info/nyct/service/ POINT (-73.90087 40.88467)#元の作業ディレクトリに戻る（move to original working directory）
setwd("..")
##Boundaries
setwd("boundary_nyc")
list_files_bound <- list.files()
shp_nam_bound <- list_files_bound[grepl(".shp",list_files_bound)]
bound_nyc <- sf::st_read(dsn=shp_nam_bound,stringsAsFactors=F,quiet=TRUE)
head(bound_nyc)## Simple feature collection with 6 features and 11 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -74.08721 ymin: 40.63981 xmax: -73.96433 ymax: 40.76364
## epsg (SRID): 4326
## proj4string: +proj=longlat +ellps=WGS84 +no_defs
## boro_code boro_ct201 boro_name cdeligibil ct2010 ctlabel ntacode
## 1 5 5000900 Staten Island E 000900 9 SI22
## 2 1 1009800 Manhattan I 009800 98 MN19
## 3 1 1010000 Manhattan I 010000 100 MN19
## 4 1 1010200 Manhattan I 010200 102 MN17
## 5 1 1010400 Manhattan I 010400 104 MN17
## 6 1 1011300 Manhattan I 011300 113 MN17
## ntaname puma shape_area shape_leng
## 1 West New Brighton-New Brighton-St. George 3903 2497010 7729.017
## 2 Turtle Bay-East Midtown 3808 1906016 5534.200
## 3 Turtle Bay-East Midtown 3808 1860938 5692.169
## 4 Midtown-Midtown South 3807 1860993 5687.802
## 5 Midtown-Midtown South 3807 1864600 5693.036
## 6 Midtown-Midtown South 3807 1890907 5699.861
## geometry
## 1 MULTIPOLYGON (((-74.07921 4...
## 2 MULTIPOLYGON (((-73.96433 4...
## 3 MULTIPOLYGON (((-73.96802 4...
## 4 MULTIPOLYGON (((-73.97124 4...
## 5 MULTIPOLYGON (((-73.97446 4...
## 6 MULTIPOLYGON (((-73.98412 4...setwd("..")既に読み込んだshapefileは以後使わないので削除します．
After importing shapefiles, we will remove downloaded shapefiles since we will not reimport them afterward.
##Subway Stations
#shapefileを展開したフォルダに移動（move to subfolder where shapefile was extracted）
setwd("sta_nyc")
#フォルダ内の全てのファイルを削除（remove all files in subfolder）
file.remove(list_files_sta)## [1] TRUE TRUE TRUE TRUE#元の作業ディレクトリに戻る（move to original working directory）
##Boundaries
setwd("..")
setwd("boundary_nyc")
file.remove(list_files_bound)## [1] TRUE TRUE TRUE TRUEsetwd("..")上に示した手順から分かるとおり，RでWEB上にあるshapefileを読み込む作業は若干煩雑です．幸いなことに，NYCのGISデータは，GeoJSON形式でも利用可能です．GeoJSON形式のデータの読み込みには，geojson_sf関数を用います．GeoJSON形式のデータは，後程ベトナムのデータを用いた解説の際に登場します．
As shown above, the procedures to import shapefiles on the website are a little complicated. Fortunately, GIS datasets of NYC are also available in GeoJSON format. To import GeoJSON file on a website or local directory, we can use geojsonsf function in geojsonsf package. We will use GeoJSON data in the later sections handling Vietnamese datasets.
# #[NOTRUN]Download GeoJSON directly
# sta_nyc <- geojsonsf::geojson_sf("https://data.cityofnewyork.us/api/geospatial/arq3-7z49?method=export&format=GeoJSON")
# bound_nyc <- geojsonsf::geojson_sf("https://data.cityofnewyork.us/api/geospatial/fxpq-c8ku?method=export&format=GeoJSON")Transform Projection
ここまでで，演習に使うshapefileを読み込み終えましたが，それらの座標は緯度経度で定義された状態です（st_is_longlat関数で，座標が緯度経度で定義されているか確認できます）．このままでは，距離計算等様々なジオプロセシングを適用できないので，m単位の座標（投影座標系）に変換を行う必要があります．
So far, we finished importing shapefiles. However, their coordinates are defined by latitude and longitude (you can check that with st_is_longlat function). Thus we have to transform their projection from geographic coordinate system (with units in degrees longitude and latitude) to projected coordinate system (with units of meters from a datum) for implementing each geoprocessing.
sf::st_is_longlat(sta_nyc)## [1] TRUEsf::st_is_longlat(bound_nyc)## [1] TRUE以下で，各shapefileについて，ニューヨーク市に対応する投影座標系であるUTM zone 18N（EPSG:6347）への投影法の変換を行います．投影変換は，st_transform関数を用いて行います．
For each shapefile, we transform its projection to UTM zone 18N, projected coordinate system corresponding to the area covering NYC. This transformation is implemented by st_transform function.
##Subway Stations
#座標系を変換（transform projection）
sta_nyc_xy <- sta_nyc %>% sf::st_transform("+init=epsg:6347")
#結果を確認（confirm a result）
sf::st_crs(sta_nyc_xy)## Coordinate Reference System:
## EPSG: 6347
## proj4string: "+proj=utm +zone=18 +ellps=GRS80 +units=m +no_defs"##Boundaries
bound_nyc_xy <- bound_nyc %>% sf::st_transform("+init=epsg:6347")
sf::st_crs(bound_nyc_xy)## Coordinate Reference System:
## EPSG: 6347
## proj4string: "+proj=utm +zone=18 +ellps=GRS80 +units=m +no_defs"XY to Point
このジオプロセシングでは，csv等のテキスト・表計算形式で与えられた座標データがある場合に，それをポイントデータへと変換する方法を示します．ここでは例として，NYCが公表している緯度経度座標付の土地評価価格データを使います．
この地価データのレコード数は80万件を超えるので，Rのデフォルト関数read.csvを使うと多大な時間がかかります．なので，読み込みにはdata.tableパッケージのfread関数を用います．
This geoprocessing procedure converts geotagged text datasets (datasets including coordinates as variables) into shapefiles. As an example, we convert geotagged dataset of assessed land value for the tax lot into point data tractable in R.
Because this land price dataset has more than 800,000 records, this procedure takes an enormous amount of time if we use default function like read.csv.
Thus we use fread function in data.table package to import.
#ダウンロードしたzipファイルを格納する一時ディレクトリを設定（set temporary directory to put downloaded zip file）
taxlot <- tempfile()
#zipファイルをダウンロード（download zip file）
curl::curl_download(url="https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nyc_pluto_20v2_csv.zip",destfile=taxlot)
#zipを解凍するためのサブフォルダを作成（create subfolder to extract csv）
dir.create("taxlot")
#サブフォルダにzipファイルを解凍（extract csv from zip file to subfolder）
unzip(zipfile=taxlot,exdir=paste0(getwd(),"/taxlot"))
#一時ディレクトリを削除（remove temporary directory）
unlink(taxlot);rm(taxlot)
#csvを展開したフォルダに移動（move to subfolder where csv was extracted）
setwd("taxlot")
#csvファイルを読み込み，緯度経度が欠損しておらず，かつ区名がBK（ブルックリン）のレコードのみ抽出
#（import csv file and extract records with no NA in longitude/latitude variables and those observed in Brooklyn "BK"）
data_taxlot_bk <- data.table::fread("pluto_20v2.csv") %>%
filter(!is.na(longitude)|!is.na(latitude)) %>% filter(borough=="BK")
#読み込んだcsvの先頭行を表示（print first 6 rows of imported csv）
head(data_taxlot_bk)## borough block lot cd ct2010 cb2010 schooldist council zipcode firecomp
## 1 BK 1549 4 303 301 3000 16 41 11233 L176
## 2 BK 850 69 307 104 3001 20 38 11220 L114
## 3 BK 574 32 306 53 1010 15 38 11231 L101
## 4 BK 5658 7 312 106 1000 20 38 11220 L114
## 5 BK 8202 42 318 966 2002 18 46 11236 L170
## 6 BK 464 7501 306 77 2001 15 39 11231 L131
## policeprct healtharea sanitboro sanitsub address zonedist1
## 1 81 3600 3 5B FULTON STREET R7D
## 2 72 6700 3 4E 725 58 STREET R6B
## 3 76 4100 3 1D 202 CONOVER STREET M2-1
## 4 66 6700 3 4E 5209 8 AVENUE R6
## 5 69 7510 3 2E 9410 FLATLANDS AVENUE R4-1
## 6 76 2500 3 1C 395 SMITH STREET R6B
## zonedist2 zonedist3 zonedist4 overlay1 overlay2 spdist1 spdist2 spdist3
## 1 C2-4 NA
## 2 NA
## 3 NA
## 4 C1-3 NA
## 5 NA
## 6 C2-4 NA
## ltdheight splitzone bldgclass landuse easements ownertype
## 1 N V1 11 0 C
## 2 N V0 11 0
## 3 N K2 5 0
## 4 N V1 11 0
## 5 N V0 11 0
## 6 N RD 3 0
## ownername lotarea bldgarea comarea resarea
## 1 NYC HOUSING PRESERVATION AND DEVELOPMENT 2000 0 NA NA
## 2 JWL REALTY LLC 2003 0 NA NA
## 3 202 CONOVER STREET LLC 2500 4314 4314 0
## 4 BROOKLYN CHINESE-AMERICAN ASSOCIATION, I NC. 1613 0 1600 0
## 5 REDDY BROOKLYN CENTER LLC 11984 0 NA NA
## 6 SMITH REGENCY CONDOMINIUM ASS 3682 13264 2560 10704
## officearea retailarea garagearea strgearea factryarea otherarea areasource
## 1 NA NA NA NA NA NA 7
## 2 NA NA NA NA NA NA 7
## 3 0 4314 0 0 0 0 2
## 4 0 0 0 0 0 1600 7
## 5 NA NA NA NA NA NA 7
## 6 0 0 0 0 0 2560 2
## numbldgs numfloors unitsres unitstotal lotfront lotdepth bldgfront bldgdepth
## 1 0 0 0 0 20.00 100.00 0 0
## 2 1 0 0 0 20.00 100.17 0 0
## 3 0 2 0 1 25.00 100.00 25 100
## 4 1 0 0 0 20.17 80.00 0 0
## 5 2 0 0 0 112.00 107.00 0 0
## 6 1 0 12 14 0.00 0.00 0 0
## ext proxcode irrlotcode lottype bsmtcode assessland assesstot exempttot
## 1 0 N 5 5 83250 83250 83250
## 2 0 N 5 5 16200 16200 0
## 3 1 N 3 0 168750 237600 0
## 4 0 N 5 5 99450 99450 0
## 5 0 N 3 5 43080 43080 0
## 6 0 N 5 5 50400 1242450 35294
## yearbuilt yearalter1 yearalter2 histdist landmark builtfar residfar commfar
## 1 0 0 0 0.00 4.20 0
## 2 0 2014 0 0.00 2.00 0
## 3 0 2018 0 1.73 0.00 2
## 4 0 2017 0 0.00 2.43 0
## 5 0 2014 0 0.00 0.75 0
## 6 0 1988 0 3.60 2.00 0
## facilfar borocode bbl condono tract2010 xcoord ycoord latitude
## 1 4.2 3 3015490004 NA 301 1006800 186469 40.67846
## 2 2.0 3 3008500069 NA 104 981560 171610 40.63771
## 3 0.0 3 3005740032 NA 53 980120 186237 40.67785
## 4 4.8 3 3056580007 NA 106 982896 172316 40.63964
## 5 2.0 3 3082020042 NA 966 1011337 173270 40.64222
## 6 2.0 3 3004647501 392 77 985434 186420 40.67836
## longitude zonemap zmcode sanborn taxmap edesignum appbbl appdate
## 1 -73.91870 17a 307 007 30603 NA <NA>
## 2 -74.00969 22a 306A033 30308 NA <NA>
## 3 -74.01489 16a 301 002 30207 NA <NA>
## 4 -74.00488 22a 306A034 31707 NA <NA>
## 5 -73.90240 23c 317 057 32405 NA <NA>
## 6 -73.99573 16c 301 049 30204 NA 3004640004 12/06/1991
## plutomapid version sanitdistrict healthcenterdistrict firm07_flag
## 1 1 20v2 3 32 NA
## 2 1 20v2 7 39 NA
## 3 1 20v2 6 38 NA
## 4 1 20v2 7 39 NA
## 5 1 20v2 18 33 NA
## 6 1 20v2 6 38 NA
## pfirm15_flag dcpedited notes
## 1 NA t NA
## 2 NA NA
## 3 1 NA
## 4 NA NA
## 5 NA NA
## 6 NA NA#元の作業ディレクトリに戻る（move to original working directory）
setwd("..")ポイントデータ作成には，st_as_sf関数を用います．ここで，coordsには読み込んだ地価データ内の座標を表す変数名を（X座標，Y座標の順），crsには座標系を指定します．座標系には，これまで読み込み・作成したshapefileに合わせてESPG:4326を指定します．
To generate point data, we use st_as_sf function, where coordinates are designated by coords (x-axis and y-axis, in order) and the coordinate system is designated by crs.
data_taxlot_bk_point <- data_taxlot_bk %>% sf::st_as_sf(coords=c("longitude","latitude"),crs="+init=epsg:4326")最後に，各種ジオプロセシングを行えるように，UTM zone 18N（EPSG:6347）への投影法の変換を行います．
Finally, we transform the projection to UTM zone 18N in order to apply each geoprocessing procedure.
data_taxlot_bk_point_xy <- data_taxlot_bk_point %>% sf::st_transform("+init=epsg:6347")
sf::st_crs(data_taxlot_bk_point_xy)## Coordinate Reference System:
## EPSG: 6347
## proj4string: "+proj=utm +zone=18 +ellps=GRS80 +units=m +no_defs"Dissolve
このセクションと次のセクション（Clip）を組み合わせて，ブルックリン区に含まれる地下鉄駅ポイントのみを抽出するジオプロセシングを行います．このセクションでは，ブルックリンに含まれる統計調査区（日本でいうと大字レベルに相当）ポリゴンを抽出した上で，1つのポリゴンに融合するジオプロセシングを実行します．
統計調査区shapefile内では，区名は変数boro_nameに格納されているので，boro_nameが「Brooklyn」であるポリゴンを，dplyrのfilter関数で抽出します．ここで示す通り，sfパッケージには，dplyrの関数群（今回はdplyr）をシームレスに使うことができる点です．
Through this section and the next section, we carry out a geoprocessing to extract subway station points included in Brooklyn. In this section, we firstly extract census tract polygons located in Brooklyn from boundary shapefile with filter function in dplyr package. As shown in this script, we can seamlessly use a series of dplyr functions combined with sf package.
bound_nyc_xy_bk_ct <- bound_nyc_xy %>% dplyr::filter(boro_name=="Brooklyn")抽出したポリゴンをプロットしてみましょう．その際，あるひとつの変数だけを選択した上でプロットを行った方が良いでしょう（変数を選択せずにplot関数を実行すると，全変数に関してのプロットがずらっと出力されてしまう）．
We will plot extracted polygons. When plotting one map only, we should select a variable to plot (otherwise all maps corresponding to all variables are plotted).
#例えばboro_nameを選択してプロットしてみる（select boro_name for example and plot）
plot(bound_nyc_xy_bk_ct %>% dplyr::select(boro_name))
次に，st_union関数を用いて，抽出した区単位のポリゴンを融合します．
Next, we firstly combine extracted census tract polygons in boundary shapefile into one polygon.
bound_nyc_xy_bk <- bound_nyc_xy_bk_ct %>% sf::st_union()融合したポリゴンをプロットしてみます．融合後のポリゴンには，属性情報は何もついていないので，区単位のポリゴンのプロットを行った場合と異なり，変数の選択は行いませんでした．
We will plot combined polygon. Since combined polygon has no attribute information, we do not select any variable, unlike uncombined polygons.
plot(bound_nyc_xy_bk)
Clip
上で融合したブルックリンポリゴンを用いて，NYC全体の地下鉄駅ポイントから，ブルックリンに含まれるポイントのみを抽出してみます．
この抽出作業を言い換えるならば，ブルックリンポリゴンと，NYC全体の地下鉄駅ポイントとの共通部分を取る作業になります．
このジオプロセシングは，st_intersection関数を用いて行います．
Using combined Brooklyn polygon, we extract subway station points located in Brooklyn from original point data.
In other words, we take the intersection between combined Brooklyn polygon and original point data. This geoprocessing can be implemented with st_intersection function.
sta_nyc_xy_bk <- sf::st_intersection(x=sta_nyc_xy,y=bound_nyc_xy_bk)ここまでで作成した各種shapefileを，一度にプロットしてみます．
この作業には，ggplot2の各種関数が必要ですが，それら関数に関する説明は最小限に留めておきます（詳細な説明は，例えば瓜生真也先生によるWEBサイトを参照）．
最初に融合済みのブルックリンポリゴンをプロットし（白地図），それに抽出された地下鉄駅ポイント（緑点）を重ね描きします．
また，ggspatialパッケージのannotation_scale関数を用いて縮尺を，annotation_north_arrow関数を用いて方位を示しておきましょう．
Finally, we plot generated shapefile on one map. This procedure requires a series of function in ggplot2, but I keep the description to a minimum (see, for example, this website).
After plotting a blank map of combined Brooklyn polygon, we overlay extracted subway station points (green).
Further, we add a scale bar and spatial-aware north arrow.
ggplot() +
geom_sf(data=bound_nyc_xy_bk,fill="white") +
geom_sf(data=sta_nyc_xy_bk,color="green",size=0.05) +
annotation_scale() +
annotation_north_arrow(location="tr")
Buffer
このジオプロセシングは，ある地物から任意の距離圏内を覆うバッファを発生するものです．それ単独で用いるというよりは，上で紹介したclip（やそれと似たintersect）と組み合わせることにより，そのバッファ内に含まれる別の地物を抽出する，という使い方をすることが多いです．
ここでは例として，各地下鉄駅から200m圏内にある地価ポイントを抽出するという作業を実装してみます．
This geoprocessing generates buffers covering an arbitrary range from each feature. We use this procedure by combining with clipping rather than using it alone. By doing so, we can extract features in another shapefile included in generated buffers.
As an example, we extract land price points within 200m from subway stations.
st_buffer関数を用いて各駅ポイントから200mのバッファを発生させます．
xでバッファ発生元になるshapefile（今回は地下鉄駅ポイント），distでバッファが覆う距離圏を指定します．
We can generate 200m buffers from stations with st_buffer function, where shapefile as an origin of buffers are designated by x, and a range of buffers are designated by dist.
sta_nyc_xy_bk_buffer_200m <- sf::st_buffer(x=sta_nyc_xy_bk,dist=units::set_units(200,m))続いて，この200mバッファを用いて，地下鉄駅から200m圏内にある地価ポイントを抽出します．先ほどと同樣，st_intersection関数を用います．ただし今回は，200mバッファが属性として持っていた駅情報が，地価ポイントの属性として新たに加わっている事に注意してください．融合されたブルックリンポリゴンのように，属性情報が何もないポリゴンでst_intersectionを実行した場合はArcGISのクリップと同様の操作になりますが，属性情報が付与されたポリゴンでst_intersectionを実行した場合はArcGISのインターセクトと同様の操作になります．
Next, we extract land price points within 200m from stations using generated buffers with st_intersection function. Please note that the attributes of subway stations included in 200m buffers are joined to extracted land price points. If we extract points with polygons that do not have an attribute, this operation is like Clip in ArcGIS while it is similar to intersect in ArcGIS if we extract with polygons having attributes.
data_taxlot_bk_point_xy_buffer_200m <- sf::st_intersection(x=data_taxlot_bk_point_xy,y=sta_nyc_xy_bk_buffer_200m)
head(data_taxlot_bk_point_xy_buffer_200m)## Simple feature collection with 6 features and 93 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 586441.1 ymin: 4503758 xmax: 586576.7 ymax: 4503932
## epsg (SRID): 6347
## proj4string: +proj=utm +zone=18 +ellps=GRS80 +units=m +no_defs
## borough block lot cd ct2010 cb2010 schooldist council zipcode firecomp
## 463 BK 931 10 306 129.02 1003 13 39 11217 L105
## 464 BK 931 22 306 129.02 1003 13 39 11217 L105
## 465 BK 931 27 306 129.02 1003 13 39 11217 L105
## 2310 BK 935 12 306 129.02 2002 13 39 11217 L105
## 2392 BK 935 27 306 129.02 2002 13 39 11217 L105
## 2393 BK 935 28 306 129.02 2002 13 39 11217 L105
## policeprct healtharea sanitboro sanitsub address zonedist1
## 463 78 2600 3 3D 450 DEAN STREET R7A
## 464 78 2600 3 3D 208 FLATBUSH AVENUE R7A
## 465 78 2600 3 3D 218 FLATBUSH AVENUE R7A
## 2310 78 2600 3 3D 12 ST MARKS AVENUE R6B
## 2392 78 2600 3 3D 40 ST MARKS AVENUE R6B
## 2393 78 2600 3 3D 42 ST MARKS AVENUE R6B
## zonedist2 zonedist3 zonedist4 overlay1 overlay2 spdist1 spdist2 spdist3
## 463 C2-4 NA
## 464 C2-4 NA
## 465 C2-4 NA
## 2310 R6A C1-4 NA
## 2392 NA
## 2393 NA
## ltdheight splitzone bldgclass landuse easements ownertype
## 463 N S5 4 0
## 464 N O9 5 0
## 465 N S2 4 0
## 2310 Y B2 1 0
## 2392 N B1 1 0
## 2393 N A4 1 0
## ownername lotarea bldgarea comarea resarea
## 463 DALEN LLC 2450 6372 1592 4780
## 464 208 FLATBUSH LLC 1914 4486 2653 1833
## 465 218-220 FLATBUSH AVENUE ASSOCIATES 1050 3030 1010 2020
## 2310 HOM THERESA 1506 3040 0 2280
## 2392 RUGG, THOMAS C 1035 2102 0 2102
## 2393 BODZIN,DANIEL,S 1040 2089 0 2089
## officearea retailarea garagearea strgearea factryarea otherarea areasource
## 463 0 1592 0 0 0 0 2
## 464 1327 1326 0 0 0 0 2
## 465 0 1010 0 0 0 0 2
## 2310 0 0 0 0 0 0 2
## 2392 0 0 0 0 0 0 2
## 2393 0 0 0 0 0 0 2
## numbldgs numfloors unitsres unitstotal lotfront lotdepth bldgfront
## 463 1 4 6 7 24.50 100.00 25.00
## 464 1 2 1 3 22.00 94.00 22.00
## 465 1 3 2 3 55.67 67.08 44.00
## 2310 1 3 2 2 18.75 80.33 18.75
## 2392 1 3 2 2 16.67 62.08 16.67
## 2393 1 3 1 1 16.67 62.42 16.67
## bldgdepth ext proxcode irrlotcode lottype bsmtcode assessland assesstot
## 463 66 N 3 N 5 2 20700 934650
## 464 92 N 3 Y 4 2 42750 499050
## 465 33 N 3 Y 3 2 14460 81300
## 2310 40 N 3 N 5 1 25620 185340
## 2392 40 N 3 N 5 2 19140 122040
## 2393 40 N 2 N 5 2 25320 129300
## exempttot yearbuilt yearalter1 yearalter2
## 463 0 1920 1999 0
## 464 0 1920 1985 0
## 465 0 1920 0 0
## 2310 1400 1880 0 0
## 2392 0 1880 0 0
## 2393 0 1880 0 0
## histdist landmark builtfar residfar
## 463 2.60 4
## 464 2.34 4
## 465 2.89 4
## 2310 Park Slope Historic District Extension II 2.02 2
## 2392 Park Slope Historic District Extension II 2.03 2
## 2393 Park Slope Historic District Extension II 2.01 2
## commfar facilfar borocode bbl condono tract2010 xcoord ycoord
## 463 0 4 3 3009310010 NA 12902 990835 187655
## 464 0 4 3 3009310022 NA 12902 990968 187467
## 465 0 4 3 3009310027 NA 12902 991045 187387
## 2310 0 2 3 3009350012 NA 12902 990598 187205
## 2392 0 2 3 3009350027 NA 12902 990845 187089
## 2393 0 2 3 3009350028 NA 12902 990859 187082
## zonemap zmcode sanborn taxmap edesignum appbbl appdate plutomapid version
## 463 16c 306 030 30401 NA <NA> 1 20v2
## 464 16c 306 030 30401 NA <NA> 1 20v2
## 465 16c 306 030 30401 NA <NA> 1 20v2
## 2310 16c 306 030 30401 NA <NA> 1 20v2
## 2392 16c 306 030 30401 NA <NA> 1 20v2
## 2393 16c 306 030 30401 NA <NA> 1 20v2
## sanitdistrict healthcenterdistrict firm07_flag pfirm15_flag dcpedited
## 463 6 38 NA NA
## 464 6 38 NA NA
## 465 6 38 NA NA
## 2310 6 38 NA NA t
## 2392 6 38 NA NA t
## 2393 6 38 NA NA t
## notes line name notes.1 objectid
## 463 NA 2-3-4 Bergen St 4-nights, 3-all other times, 2-all times 4
## 464 NA 2-3-4 Bergen St 4-nights, 3-all other times, 2-all times 4
## 465 NA 2-3-4 Bergen St 4-nights, 3-all other times, 2-all times 4
## 2310 NA 2-3-4 Bergen St 4-nights, 3-all other times, 2-all times 4
## 2392 NA 2-3-4 Bergen St 4-nights, 3-all other times, 2-all times 4
## 2393 NA 2-3-4 Bergen St 4-nights, 3-all other times, 2-all times 4
## url geometry
## 463 http://web.mta.info/nyct/service/ POINT (586511.8 4503932)
## 464 http://web.mta.info/nyct/service/ POINT (586552.9 4503876)
## 465 http://web.mta.info/nyct/service/ POINT (586576.7 4503852)
## 2310 http://web.mta.info/nyct/service/ POINT (586441.1 4503795)
## 2392 http://web.mta.info/nyct/service/ POINT (586516.8 4503760)
## 2393 http://web.mta.info/nyct/service/ POINT (586521.1 4503758)最後に，融合済みのブルックリンポリゴンのプロットに，200mバッファと抽出された地価ポイント（赤点）を重ね描きします．
また，ggspatialパッケージのannotation_scale関数を用いて縮尺を，annotation_north_arrow関数を用いて方位を示しておきましょう．
Finally, we plot generated shapefiles.
After plotting a blank map of combined Brooklyn polygon, we overlay 200m buffers and extracted land price point (red) in order.
ggplot() +
geom_sf(data=bound_nyc_xy_bk,fill="white") +
geom_sf(data=sta_nyc_xy_bk_buffer_200m,fill="green") +
geom_sf(data=data_taxlot_bk_point_xy_buffer_200m,color="red",size=0.0005) +
annotation_scale() +
annotation_north_arrow(location="tr")
Spatial Join
このジオプロセシングの用途は主に以下の2つです．
- 1．ポリゴン属性の付与：あるポリゴンに含まれる各ポイントに関して，そのポリゴンの属性を付与する
- 2．最寄り地物情報の付与：2つのshapefileの属性情報を地理的距離に基づいて結合する
いずれのジオプロセシングも，st_join関数を用いて行うことができます．
The main usage of this geoprocessing is as follows:
- Joining attributes of a polygon to each point included in it
- Joining the attributes of two shapefiles based on the location of the features in the shapefiles (attributes of the nearest feature are joined)
Both geoprocessing procedures can be done with st_join function.
1. Join Polygon Attributes
例として，先ほど作成した地下鉄駅ポイントに，それらが含まれる「結合前の」国勢調査区境界ポリゴンの属性を付与する作業を実装してみます．
st_joinでxに指定されるのは，属性を付与される側のshapefile（今回の例では地下鉄駅ポイント）です．
また，オプションjoinには，st_intersectsを指定してください．
For example, we implement the geoprocessing that joins the attributes of uncombined census tract boundary polygons to extracted subway stations.
Since st_join function implements left_join, we designate target shapefile (attributes will be joined to this shapefile) by x. In this case, we select “st_intersects” option in join argument.
#各駅ポイントが含まれる区単位ポリゴンの属性を付与（join polygon attributes to station points）
sta_nyc_xy_bk_bound <- sf::st_join(x=sta_nyc_xy_bk,y=bound_nyc_xy_bk_ct,join=st_intersects)
head(sta_nyc_xy_bk_bound)## Simple feature collection with 6 features and 16 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 586034.4 ymin: 4502125 xmax: 594678.2 ymax: 4504915
## epsg (SRID): 6347
## proj4string: +proj=utm +zone=18 +ellps=GRS80 +units=m +no_defs
## line name
## 4 2-3-4 Bergen St
## 5 3-4 Pennsylvania Ave
## 8 A-C Kingston - Throop Aves
## 12 J-Z Van Siclen Ave
## 13 J-Z Norwood Ave
## 15 B-D-N-Q-R DeKalb Ave
## notes
## 4 4-nights, 3-all other times, 2-all times
## 5 4-nights, 3-all other times
## 8 A-nights, C-all other times
## 12 Z-rush hours AM westbound, PM eastbound, J-all other times
## 13 Z-rush hours AM westbound, PM eastbound, J-all other times
## 15 B-weekdays and evenings, D-all other times, N-nights, R-all other times, Q-all times
## objectid url boro_code boro_ct201 boro_name
## 4 4 http://web.mta.info/nyct/service/ 3 3012902 Brooklyn
## 5 5 http://web.mta.info/nyct/service/ 3 3112600 Brooklyn
## 8 8 http://web.mta.info/nyct/service/ 3 3027100 Brooklyn
## 12 12 http://web.mta.info/nyct/service/ 3 3119800 Brooklyn
## 13 13 http://web.mta.info/nyct/service/ 3 3117400 Brooklyn
## 15 15 http://web.mta.info/nyct/service/ 3 3003100 Brooklyn
## cdeligibil ct2010 ctlabel ntacode ntaname puma
## 4 I 012902 129.02 BK37 Park Slope-Gowanus 4005
## 5 E 112600 1126 BK85 East New York (Pennsylvania Ave) 4007
## 8 E 027100 271 BK61 Crown Heights North 4006
## 12 I 119800 1198 BK82 East New York 4008
## 13 E 117400 1174 BK83 Cypress Hills-City Line 4008
## 15 I 003100 31 BK68 Fort Greene 4004
## shape_area shape_leng geometry
## 4 1372594 5163.445 POINT (586619.3 4503836)
## 5 1813472 5460.421 POINT (593412.1 4502125)
## 8 2767087 9484.846 POINT (589505.7 4503765)
## 12 3527652 8959.335 POINT (593666.3 4503607)
## 13 2295156 6842.142 POINT (594678.2 4504007)
## 15 3102639 7674.199 POINT (586034.4 4504915)2. Join Attrubtes of the Nearest Feature (Point)
例として，ブルックリンに位置する各地価ポイントから最寄りの駅ポイントの属性を付与する作業を実装してみます．今回は，ArcGISの空間結合ツールと同じアウトプットを作成するため，最寄り駅ポイントの属性を付与するだけでなく，最寄り駅ポイントまでの直線距離を計算する作業も合わせて行いたいと思います．
For example, we carry out the geoprocessing procedure to join the attributes of the nearest subway station to each land price point. To exactly mimic ArcGIS tool, we additionally calculate Euclidean distance to the nearest station as well as joining attributes.
これは完全にsfパッケージの仕様上の注意なのですが，直線距離の計算を行うためには，各shapefileのgeometry変数（座標を格納した変数）を，dplyrでハンドリング可能な形式の変数として格納し直す作業が必要です．
変数geometryからXY座標を取り出す際には，st_coordinates関数を用います．
This is just a notice due to the specification of sf package, but in order to calculate Euclidean distance, we have to generate coordinates variables tractable for data handling with dplyr from a geometry variable (a combination of x-axis and y-axis) in each shapefile.
To extract XY coordinates from geometry, we use st_coordinates function.
taxlot_coords <- sf::st_coordinates(x=data_taxlot_bk_point_xy)
head(taxlot_coords)## X Y
## 1 591380.2 4503626
## 2 583741.4 4499012
## 3 583252.0 4503463
## 4 584146.0 4499231
## 5 592808.3 4499621
## 6 584870.4 4503537sta_coords <- sf::st_coordinates(x=sta_nyc_xy)
head(sta_coords)## X Y
## 1 585198.5 4509281
## 2 584442.3 4508023
## 3 585767.5 4512804
## 4 586619.3 4503836
## 5 593412.1 4502125
## 6 592600.3 4526536続いて，取り出したXY座標を，各shapefileへ新たな変数として追加します．
Next, we add extracted coordinates to each shapefile as new variables.
data_taxlot_bk_point_xy_coords <- data_taxlot_bk_point_xy %>% dplyr::mutate(X_taxlot=taxlot_coords[,"X"],Y_taxlot=taxlot_coords[,"Y"])
sta_nyc_xy_coords <- sta_nyc_xy %>% dplyr::mutate(X_sta=sta_coords[,"X"],Y_sta=sta_coords[,"Y"])ここで，st_join関数を用いた，各地価ポイントから最寄りの駅ポイントの属性の付与を行います． ポリゴン属性付与とは異なり，オプションjoinには，st_nearest_featureを指定することに注意してください．
Here we join the attributes of the nearest subway station to each land price point with st_join function.
Note that we select “st_nearest_feature” option in join argument, unlike the case of joining polygon attributes.
#各地価ポイントから最寄りの駅ポイントの属性を付与（join the attributes of the nearest station to land price points）
data_taxlot_bk_point_xy_coords_near_sta <- sf::st_join(x=data_taxlot_bk_point_xy_coords,y=sta_nyc_xy_coords,st_nearest_feature)最後に，最寄り駅ポイントまでの直線距離を計算してみましょう．
座標はm単位で定義されていますので，単に，地価ポイント・駅ポイントの座標からユークリッド距離を計算するだけです．
Finally, we can calculate Euclidean distance to the nearest station in ease thanks to coordinates with units of meters.
data_taxlot_bk_point_xy_coords_near_sta_dist <- data_taxlot_bk_point_xy_coords_near_sta %>% dplyr::mutate(d_sta=sqrt((X_taxlot-X_sta)^2+(Y_taxlot-Y_sta)^2))
head(data_taxlot_bk_point_xy_coords_near_sta_dist)## Simple feature collection with 6 features and 98 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 583252 ymin: 4499012 xmax: 592808.3 ymax: 4503626
## epsg (SRID): 6347
## proj4string: +proj=utm +zone=18 +ellps=GRS80 +units=m +no_defs
## borough block lot cd ct2010 cb2010 schooldist council zipcode firecomp
## 1 BK 1549 4 303 301 3000 16 41 11233 L176
## 2 BK 850 69 307 104 3001 20 38 11220 L114
## 3 BK 574 32 306 53 1010 15 38 11231 L101
## 4 BK 5658 7 312 106 1000 20 38 11220 L114
## 5 BK 8202 42 318 966 2002 18 46 11236 L170
## 6 BK 464 7501 306 77 2001 15 39 11231 L131
## policeprct healtharea sanitboro sanitsub address zonedist1
## 1 81 3600 3 5B FULTON STREET R7D
## 2 72 6700 3 4E 725 58 STREET R6B
## 3 76 4100 3 1D 202 CONOVER STREET M2-1
## 4 66 6700 3 4E 5209 8 AVENUE R6
## 5 69 7510 3 2E 9410 FLATLANDS AVENUE R4-1
## 6 76 2500 3 1C 395 SMITH STREET R6B
## zonedist2 zonedist3 zonedist4 overlay1 overlay2 spdist1 spdist2 spdist3
## 1 C2-4 NA
## 2 NA
## 3 NA
## 4 C1-3 NA
## 5 NA
## 6 C2-4 NA
## ltdheight splitzone bldgclass landuse easements ownertype
## 1 N V1 11 0 C
## 2 N V0 11 0
## 3 N K2 5 0
## 4 N V1 11 0
## 5 N V0 11 0
## 6 N RD 3 0
## ownername lotarea bldgarea comarea resarea
## 1 NYC HOUSING PRESERVATION AND DEVELOPMENT 2000 0 NA NA
## 2 JWL REALTY LLC 2003 0 NA NA
## 3 202 CONOVER STREET LLC 2500 4314 4314 0
## 4 BROOKLYN CHINESE-AMERICAN ASSOCIATION, I NC. 1613 0 1600 0
## 5 REDDY BROOKLYN CENTER LLC 11984 0 NA NA
## 6 SMITH REGENCY CONDOMINIUM ASS 3682 13264 2560 10704
## officearea retailarea garagearea strgearea factryarea otherarea areasource
## 1 NA NA NA NA NA NA 7
## 2 NA NA NA NA NA NA 7
## 3 0 4314 0 0 0 0 2
## 4 0 0 0 0 0 1600 7
## 5 NA NA NA NA NA NA 7
## 6 0 0 0 0 0 2560 2
## numbldgs numfloors unitsres unitstotal lotfront lotdepth bldgfront bldgdepth
## 1 0 0 0 0 20.00 100.00 0 0
## 2 1 0 0 0 20.00 100.17 0 0
## 3 0 2 0 1 25.00 100.00 25 100
## 4 1 0 0 0 20.17 80.00 0 0
## 5 2 0 0 0 112.00 107.00 0 0
## 6 1 0 12 14 0.00 0.00 0 0
## ext proxcode irrlotcode lottype bsmtcode assessland assesstot exempttot
## 1 0 N 5 5 83250 83250 83250
## 2 0 N 5 5 16200 16200 0
## 3 1 N 3 0 168750 237600 0
## 4 0 N 5 5 99450 99450 0
## 5 0 N 3 5 43080 43080 0
## 6 0 N 5 5 50400 1242450 35294
## yearbuilt yearalter1 yearalter2 histdist landmark builtfar residfar commfar
## 1 0 0 0 0.00 4.20 0
## 2 0 2014 0 0.00 2.00 0
## 3 0 2018 0 1.73 0.00 2
## 4 0 2017 0 0.00 2.43 0
## 5 0 2014 0 0.00 0.75 0
## 6 0 1988 0 3.60 2.00 0
## facilfar borocode bbl condono tract2010 xcoord ycoord zonemap zmcode
## 1 4.2 3 3015490004 NA 301 1006800 186469 17a
## 2 2.0 3 3008500069 NA 104 981560 171610 22a
## 3 0.0 3 3005740032 NA 53 980120 186237 16a
## 4 4.8 3 3056580007 NA 106 982896 172316 22a
## 5 2.0 3 3082020042 NA 966 1011337 173270 23c
## 6 2.0 3 3004647501 392 77 985434 186420 16c
## sanborn taxmap edesignum appbbl appdate plutomapid version
## 1 307 007 30603 NA <NA> 1 20v2
## 2 306A033 30308 NA <NA> 1 20v2
## 3 301 002 30207 NA <NA> 1 20v2
## 4 306A034 31707 NA <NA> 1 20v2
## 5 317 057 32405 NA <NA> 1 20v2
## 6 301 049 30204 NA 3004640004 12/06/1991 1 20v2
## sanitdistrict healthcenterdistrict firm07_flag pfirm15_flag dcpedited notes.x
## 1 3 32 NA NA t NA
## 2 7 39 NA NA NA
## 3 6 38 NA 1 NA
## 4 7 39 NA NA NA
## 5 18 33 NA NA NA
## 6 6 38 NA NA NA
## X_taxlot Y_taxlot line name notes.y
## 1 591380.2 4503626 A-C Ralph Ave A-nights, C-all other times
## 2 583741.4 4499012 N 8th Ave N-all times
## 3 583252.0 4503463 F-G Smith - 9th Sts F,G-all times
## 4 584146.0 4499231 N 8th Ave N-all times
## 5 592808.3 4499621 L Canarsie - Rockaway Pkwy L-all times
## 6 584870.4 4503537 F-G Carroll St F,G-all times
## objectid url X_sta Y_sta
## 1 38 http://web.mta.info/nyct/service/ 591203.5 4503664
## 2 77 http://web.mta.info/nyct/service/ 583590.6 4498706
## 3 375 http://web.mta.info/nyct/service/ 584862.9 4503014
## 4 77 http://web.mta.info/nyct/service/ 583590.6 4498706
## 5 219 http://web.mta.info/nyct/service/ 592848.4 4500113
## 6 394 http://web.mta.info/nyct/service/ 584934.3 4503751
## geometry d_sta
## 1 POINT (591380.2 4503626) 180.7055
## 2 POINT (583741.4 4499012) 340.5397
## 3 POINT (583252 4503463) 1672.4063
## 4 POINT (584146 4499231) 764.3173
## 5 POINT (592808.3 4499621) 494.2243
## 6 POINT (584870.4 4503537) 222.8846Geometry Calculation, Summarize, Table Join
ここでは，ブルックリンに含まれる地価ポイントを用いて，各種集計変数を計算する作業を示します．
今回計算する集計変数は，以下の通りです．
- 統計調査区別地価ポイント数
- 統計調査区別100m2あたり地価ポイント数（地価ポイント密度）
In this section, we calculate several summary values using land price points located in Brooklyn. We calculate the following variables.
- Total number of land price points by census tracts
- Total number of land price points by census tracts per 100m2 (land price point density)
最初に計算する集計変数は，統計調査区別地価ポイント数です．
地価ポイントデータには既に，集計単位である統計調査区番号は付与された状態なので，geometry属性が付与された（shapefile）状態で集計を行うのではなく，geometry属性を削除したdata.frame形式に変換した後で集計変数を計算します．
geometry属性の削除は，st_set_geomerty関数でvalueをNULLと指定することで行えます．
Firstly we calculate the total number of land price points by census tracts.
Since land price point data already has a variable about census tract, we do not have to use geometry information afterwards.
We can remove geometry information with st_set_geomerty function, where we set argument value NULL.
data_taxlot_bk_point_xy_df <- data_taxlot_bk_point_xy %>% sf::st_set_geometry(value=NULL) 次に，dplyrパッケージの関数をいくつか用いて，統計調査区別地価ポイント数を計算します．具体的には，調査区コードct2010に基づき，レコードをgroup_by関数でグループ化した後，summarise関数で調査区別地価ポイント数の合計を計算します．その後，集計結果をdata.frameに変換します．
後程結合する国勢調査区境界ポリゴンデータでは，調査区コードct2010に対応する変数の名前はctlabelとなっており，かつ，変数の型は文字列なので，集計結果data.frame内でも，同じ名前・型の変数を作成します．
変数名の重複を避けるために，集計結果data.frameから変数ct2010を削除します（調査区境界ポリゴンデータの中に，名前は同じct2010だが，全く定義が異なる変数が含まれている）．
Next, we calculate the total number with several functions in dplyr package. Specifically, we firstly group records by variable ct2010 (census tract IDs) with group_by function, and we calculate the total number of points with summarise function.
Since the name of a variable corresponding to ct2010 is “ctlabel” in census tract boundary data we join later and the type of ctlabel is character, we also generate variable ctlabel in summarised data whose type is character.
For avoiding the duplication of variables, we remove variable ct2010 from summarised data (boundary data also has a variable named “ct2010” based on completely irrelevant definition).
data_taxlot_bk_point_xy_df_sum <- data_taxlot_bk_point_xy_df %>% dplyr::group_by(ct2010) %>%
dplyr::summarise(taxlot_sum=n()) %>% as.data.frame %>%
dplyr::mutate(ctlabel=as.character(ct2010)) %>%
dplyr::select(-ct2010)統計調査区別100m2あたり地価ポイント数（地価ポイント密度）を計算するため，dplyrのleft_join関数を用いて，国勢調査区境界ポリゴンに集計結果を左結合します．左結合のキー変数にはctlabelを指定します．
In order to calculate the total number of land price points by census tracts per 100m2 (land price point density), we join summarised data to census tract boundary polygons by variable ctlabel with left_join function in dplyr package.
bound_nyc_xy_bk_ct_taxlot_sum <- dplyr::left_join(x=bound_nyc_xy_bk_ct,y=data_taxlot_bk_point_xy_df_sum,by=c("ctlabel"))地価ポイント密度を計算するには，各区の面積が必要です．
shapefileに含まれる各ポリゴンの面積を計算する際には，st_area関数を用います．今回は，unitsパッケージのset_units関数を用いて，面積の単位を1km2に指定します．
For calculating land price point density, the area of each boundary polygon is necessary. To calculate area, we use st_area package. We set units km2 with set_units function in units function.
bound_nyc_xy_bk_ct_taxlot_sum_area <- sf::st_area(x=bound_nyc_xy_bk_ct_taxlot_sum) %>% units::set_units(km^2) dplyrのmutate関数を用いて，先ほど計算した面積（数値ベクトル）をarea_ctという変数に代入し，その変数で区別地価ポイント数taxlot_sumを割り，統計調査区別1km2あたり地価ポイント数を計算します．最後に，この1km2あたりポイント数を100で割ることで，統計調査区別100m2あたり地価ポイント数を計算します．
ちなみに，作成した変数taxlot_sumは，units（単位）が付与された数値ベクトル（区別面積）を元にした変数なので，それ自身にもunitsが付与されている状態です．unitsがついたままの状態では，後程行う集計変数のマッピングができないので，unitsパッケージのset_units関数を用いて，unitsを削除します．
Using mutate function in dplyr package, we define a variable area_ct with the calculated area. After that, we define a new variable taxlot_dens defined with taxlot_sum/area_ct. Finally, we divide taxlot_dens by 100 and remove units (this is because we cannot apply several functions in ggplot2 to variables with units) and assign the result to variable taxlot_dens_wo_units.
bound_nyc_xy_bk_ct_taxlot_sum_area_dens <- bound_nyc_xy_bk_ct_taxlot_sum %>%
dplyr::mutate(area_ct=bound_nyc_xy_bk_ct_taxlot_sum_area) %>%
dplyr::mutate(taxlot_dens=taxlot_sum/area_ct) %>%
dplyr::mutate(taxlot_dens_wo_units=units::set_units(taxlot_dens,value=NULL)/100 %>% round(digits=2))国勢調査区境界ポリゴンを地価ポイント密度で階級分けした上でプロットした後に，地下鉄駅ポイントを重ね描きします．
After plotting census tract boundary polygons setting fill color based on the value of land price point density, we overlay subway station points.
ggplot() +
geom_sf(data=bound_nyc_xy_bk_ct_taxlot_sum_area_dens,aes(fill=ggplot2::cut_interval(taxlot_dens_wo_units,n=5))) +
scale_fill_brewer(palette="Spectral",direction=-1,name="Lot Density[100m^2]") +
geom_sf(data=sta_nyc_xy_bk,color="black",size=1) +
annotation_scale() +
annotation_north_arrow(location="tr")
Raster Data
ここまでで紹介してきたジオプロセシングの方法は，全てベクタデータを対象にしたものでした．以下では，ラスタデータとベクタデータを組み合わせたジオプロセシングの方法を見て行きます．なお，今回扱うラスタデータを用いたジオプロセシングには，30分以上時間がかかるものや，PCの空き容量が十分ないと（少なくとも16GB以上）実行できないものもありますので，ご注意ください．
Geoprocessing techniques I explained above were for vector data. Hereafter, we will see the techniques convining raster data with vector data. Please note that several geoprocessing procedures below take more than 30 minutes to implement or require sufficient free disk space (at least more than 16GB).
1. Summarizing Rasters with Polygons
以下では，ラスタデータの情報を集計し，ポリゴンに結合する方法を見ていきます．具体的には，ベトナムの行政区域ポリゴンに，夜間光量ラスタデータを集計する作業を示します．属性がついていないポリゴンデータを使っても良いのですが，社会情勢を考慮し，COVID-19の地域別発生件数が属性として付与されたポリゴンデータに結合してみます．
Here we will learn how to summarize raster data with polygons. Specifically, we summarize the night-time light raster data with Vietnamese provincial administration polygons. Considering the social condition in 2020, we dare to summarize the raster data with the administration polygons having the number of COVID-19 cases as attributes, although it is also fine even if we use polygons without the attribute.
WEBページ上から，夜間光量のラスタデータが入った圧縮ファイルをダウンロードし，作業ディレクトリに解凍します．ここでは単に，shapefileの時と同様の操作を行っています．解凍が完了したら，rasterパッケージのraster関数を用いて，データを読み込みます．ただし，データのダウンロードには，早くても25分かかることに留意してください．
We download zipped night-time light data from the web page and unzip on our working directory. We simply repeat the same exercise with that for shapefiles. After that, we import the raster data with raster function in raster package. Please note that it takes at least 25 minutes to download the data.
#ダウンロードしたラスタが入ったtarファイルを格納する一時ファイルを設定（set temporary directory to put downloaded tar file）
nl <- tempfile()
#tarファイルをダウンロード（2013年版）（download 2013 version tar file）
curl::curl_download(url="https://ngdc.noaa.gov/eog/data/web_data/v4composites/F182013.v4.tar",destfile=nl)
#現在の作業ディレクトリにtarファイルを解凍（unter on current working directory）
untar(tarfile=nl,exdir=getwd())
R.utils::gunzip("F182013.v4c_web.stable_lights.avg_vis.tif.gz",overwrite=T)
#一時ファイルを削除（delete temporary file）
unlink(nl);rm(nl)
#夜間光量（全世界）を読み込み（import global night-time light data）
nightlight <- raster::raster("F182013.v4c_web.stable_lights.avg_vis.tif") 次に，COVID-19の地域別発生件数が集計されたベトナムの行政区域ポリゴンをWEBページからダウンロード・読み込みします．このポリゴンはGeoJSON形式でも利用可能なので，今回はgeojson_sf関数を用いて読み込みましょう．ウィキペディアによると，GeoJSONは，JavaScript Object Notation (JSON) を用いて空間データをエンコードし非空間属性を関連付けるファイルフォーマットで，sfオブジェクトとして読み込んだ後は，shapefileと同じように操作することができます．
読み込んだGeoJSONの地理座標系は既にWGS84（EPSG:4326）で定義されているので，後のジオプロセシングを行うために，座標系を投影座標系に変換します．今回は，UTM zone 49N（EPSG:3149）を用います．最後に，なぜか感染者数が重複して観測される地域のデータを除外します．
Next, we download and import Vietnamese provincial administration polygons from the web page. The polygons are also available in GeoJSON format, we can import them with geojson_sf function. According to Wikipedia GeoJSON is an open standard format designed for representing simple geographical features, along with their non-spatial attributes. After importing GeoJSON datasets as sf objects, we can handle them as with shapefile. Since the geographic coordinate system is already defined as WGS84 (EPSG:4326), we transform the system to projected one for the later geoprocessing procedures. We set the system to UTM zone 49N (EPSG:3149). Finally, we drop several regions whose number of infected people is duplicately observed due to some reason.
#直接GeoJSONをダウンロード・読み込み（download and import GeoJSON directly）
viet_covid <- geojsonsf::geojson_sf("https://data.opendevelopmentmekong.net/dataset/100d98c8-4ffc-4b5b-b2a5-c6b71d2eff66/resource/c8b14875-6ab9-4a3e-99e9-8352f695cd77/download/covid-19-provinces.geojson") %>%
sf::st_transform(crs="+init=epsg:3149") %>%
dplyr::filter(!(Code%in%c(48,56,77,91))) %>%
dplyr::mutate(infected=as.numeric(infected))
ggplot2::ggplot() +
ggplot2::geom_sf(data=viet_covid,aes(fill=ggplot2::cut_interval(infected,n=6))) +
ggplot2::scale_fill_brewer(palette="Spectral",direction=-1,name="No. of Infected") +
ggspatial::annotation_scale() +
ggspatial::annotation_north_arrow(location="tr")
読み込んだラスタデータから，必要な領域だけを切り出すときには，rasterパッケージのcrop関数を用います．引数xには切り抜きを行いたいラスタデータを，引数yには，切り抜き領域を指定します．最初に，ベトナムを含む大まかな範囲を緯度経度座標で指定します．最初に大まかな範囲で切り抜きを行うのは，世界版の夜間光量データのサイズが，切り抜きを行わずに座標系を定義・変換するには大きすぎるためです．
ラスタデータの座標系の定義は，projectRaster関数で行えます．今回は，WGS84（EPSG:4326）で定義します．地理座標系を定義した後，後のジオプロセシングを行うために，投影座標系に変換します．最後に，行政区域ポリゴンの正確な範囲で，夜間光量データを切り抜きます．その際，引数yには行政区域ポリゴンを指定します．
We use crop function when extracting the area necessary for geoprocessing procedures. We assign the cropped raster data to argument x and the cropping extent expressed by coordinates to argument y. We firstly crop the data with rough extent since the size is too large to define and transform the coordinate system directly.
Then we define the geographic coordinate system for the data. We set the system to WGS84 (EPSG:4326). After that, we transform the system to projected one for the later geoprocessing procedures. Finally, again we crop the raster data exactly covered with Vietnamese provincial administration by assigning the polygon data y.
#地理座標系を定義し，ベトナムの範囲だけ抽出する（define the geographic coordinate system, transform to the projected coordinate system, and extract the raster data covering Vietnam and ）
viet_nightlight <- raster::crop(x=nightlight,y=extent(102,110,8,24)) %>%
raster::projectRaster(crs=sp::CRS("+init=epsg:4326")) %>%
raster::projectRaster(crs=sp::CRS("+init=epsg:3149")) %>%
raster::crop(y=viet_covid)
plot(viet_nightlight)
ラスタデータのポリゴンへの集計には，rasterパッケージのextract関数を用います．引数xには集計したいラスタデータ，引数yには集計先のポリゴンデータを指定します．また，どのような集計を行うかはfunで指定します．今回は，メッシュ毎に夜間光量の合計を求めたいので，sumで指定します．ポリゴンのジオメトリを保持したい場合には，spにTRUEを指定してください．最後に，夜間光量の合計を表す変数名layerをntlに変更します（ggplot2にlayerという名前の関数があるので）．
We can summarize the raster data with polygons with extract function. We assign the raster data to argument x and polygons summarizing rasters to argument y. We specify the function to summarize the values with argument fun. Since we calculate the total amount of night-time light, we use “sum”. If we want to keep the geometry of polygons, we set sp to TRUE. Finally, we rename the variable name since “layer” is a function of ggplot2.
#メッシュ毎に夜間光量を集計する．足しあわせる過程でsfオブジェクトがspオブジェクトに変わるので元に戻す（summarize rasters with polygons. Convert output generated as sp object to sf object）
viet_covid_nightlight <- raster::extract(x=viet_nightlight,y=viet_covid,fun=sum,na.rm=T,sp=T) %>%
sf::st_as_sf(.) %>% sf::st_set_crs(sf::st_crs(viet_covid)$epsg) %>%
rename(ntl=layer)
#夜間光量のメッシュプロット
ggplot2::ggplot() +
ggplot2::geom_sf(data=viet_covid_nightlight,aes(fill=ntl),size=0.05) +
ggspatial::annotation_scale() +
ggspatial::annotation_north_arrow(location="tr")
ここでは簡単な分析として，夜間光量とCOVID-19の感染者数との関係を見てみます．縦軸を感染者数，横軸を夜間光量としてプロットを描くと，多数の感染者数が0である地域と少数の感染者が発生した地域があることが分かります（幸いにもベトナムは感染拡大抑制に最も成功した国の1つです）．このような分布形を踏まえれば，線形回帰モデルではなく，ポアソン回帰モデルを用いて関係を検証することが妥当そうです． glm関数を用いてポアソン回帰分析を行うと，夜間光量を表す変数ntlの回帰係数はとても強く正に有意です．擬似決定係数の値もまあまあなので，経済活動が活発な地域ではCOVID-19の感染者数が多いという関係があることが分かります．
As a simple example, we can test the correlation between the number of infected people and the strength of night-time light. From the plot with the number of infected people as the vertical axis and ntl, the strength of night-time light, as the horizontal axis, we can observe the majority of regions whose number of infected people is 0 and a few regions with the positive number of infected people (fortunately, Vietnam is one of the countries that succeeded most to prevent the spread of infection). Considering this distribution, it might be plausible to examine the association by using Poisson regression rather than linear regression.
The result of Poisson regression shows that the regression coefficient of ntl is strongly and positively significant, and the value of Pseudo-R-square is moderate. Thus we can say that a region with a higher level of economic activity has more infected people.
plot(viet_covid_nightlight$ntl,viet_covid_nightlight$infected)
pois_reg <- glm(infected~ntl,data=viet_covid_nightlight,family=poisson)
summary(pois_reg)##
## Call:
## glm(formula = infected ~ ntl, family = poisson, data = viet_covid_nightlight)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -8.8213 -2.4955 -1.8346 0.0782 9.1215
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.862e-01 1.012e-01 3.814 0.000136 ***
## ntl 3.511e-05 1.399e-06 25.091 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1217.92 on 58 degrees of freedom
## Residual deviance: 681.78 on 57 degrees of freedom
## AIC: 784.62
##
## Number of Fisher Scoring iterations: 6PR2 <- 1-(pois_reg$deviance/pois_reg$null.deviance)
PR2## [1] 0.44020642. Elevation
よりベーシックなラスタデータのひとつに，標高データがあります．ここでは，標高データの基本的なハンドリング手順に加え，標高データを用いた空間解析の一例として，最小コスト距離解析を紹介します．処理能力の制約から，今回はベトナム南部に焦点を絞ってジオプロセシングの手順を示します．
まずは，標高データをWEBページからダウンロード・読み込みします．その後，ベトナム南部の領域に対応する標高データを抽出し，座標系を変換します．
As a more basic raster data, we may use the elevation data. Here we will learn how to handle the elevation data. We also learn the least-cost path analysis as an example of spatial analysis using them. Due to the computational capacity, I will illustrate geoprocessing procedures focusing on the southern part of Vietnam. Firstly, we download and import the elevation data from the web page. Then we extract the raster data covering the southern part and assign the coordinate system.
#ラスタをダウンロード・読み込み，ベトナム南部のみに絞り込み，座標系を変換（download and import elevation data, crop raster data covering southern part, and assign coordinate system）
elev_s_viet <- raster::raster("https://github.com/globalmaps/gmvn20/raw/master/el_vnm.tif") %>%
raster::crop(y=extent(102,110,8,14))
#(I separately convert projection due to some bug in pipe operation)
raster::crs(elev_s_viet) <- sp::CRS("+init=epsg:4326")
#値が入っていない・水面をNAに置換（set value of rasters on water surface or rasters that do not have value to NA）
elev_s_viet[elev_s_viet>4000] <- NA
elev_s_viet[elev_s_viet<0] <- NA
#直角座標系に変換（convert coordinate system to projected one）
elev_s_viet_xy <- elev_s_viet %>%
raster::projectRaster(crs=CRS("+init=epsg:3149"))
plot(elev_s_viet_xy)
単純なジオプロセシングの例として，各地点の傾斜を計算してみましょう．傾斜の計算は，rasterパッケージのterrain関数で行うことができます．引数optにはslopeを，傾斜を度数法で表したい場合には，引数unitにdegreesを指定します．
As an example of simple geoprocessing, we can calculate the slope using terrain function. We assign “slope” to argument opt and “degrees” to the argument unit if we want to express the slope with degree measure.
#傾斜区分図を度数法で作成（plot inclination classification in degree measure）
elev_s_viet_xy_slope <- elev_s_viet_xy %>% raster::terrain(opt="slope",unit="degrees")
plot(elev_s_viet_xy_slope)
標高データをインプットとしたコスト距離解析の準備をします．ここでは具体的に，ホーチミンのサイゴン駅を出発地点，ニャチャンのニャチャン駅・カントーのカントー国際空港のそれぞれを到着地点として，コスト最小の経路を探索します．ここでのコストとは，傾斜や地表面の条件，地下構造などを指しますが，今回は傾斜をコスト要因として考慮した経路探索を行います．
まず，各駅・空港の緯度経度座標をsfオブジェクトに変換し，到着点・出発点オブジェクトに格納します．コスト距離解析に用いるパッケージmovecostはsfオブジェクトには対応していないので，格納時にas関数を用いてspオブジェクトに変換します．
We prepare for the least-cost path analysis using the elevation data as inputs. Specifically, we search the routes from Saigon station in Ho Chi Minh to Nha Trang station in Nha Trang or Can Tho International Airport in Can Tho, respectively. In this context, “cost” means the slope, ground-surface conditions, or underground structure. Here we will search for routes considering the slope as a cost factor. Firstly, we convert the coordinates of each point to sf object and assign to origin or destination object. Since movecost package that we will use cannot handle sf object, we convert the objects to sp objects with as function.
#出発・到着地点の緯度経度をsfオブジェクトに変換後，spオブジェクトに変換（convert coordinates of each origin or destination point to sf object and re-convert them to sp objects）
orig_point <- data.frame(matrix(c(106.677,10.7822),nrow=1,ncol=2)) %>%
magrittr::set_colnames(c("x","y")) %>%
sf::st_as_sf(coords=c("x","y"),crs="+init=epsg:4326") %>%
sf::st_transform(crs="+init=epsg:3149") %>%
dplyr::mutate(objname="Saigon station") %>%
methods::as("Spatial")
dest_point <- data.frame(matrix(c(105.711944,10.085278,109.184167,12.248333),nrow=2,ncol=2,byrow=T)) %>%
magrittr::set_colnames(c("x","y")) %>%
sf::st_as_sf(coords=c("x","y"),crs="+init=epsg:4326") %>%
sf::st_transform(crs="+init=epsg:3149") %>%
dplyr::mutate(objname=c("Can Tho International Airport","Nha Trang station")) %>%
methods::as("Spatial")解析に用いる関数は，movecost関数です．引数dtmには標高ラスタデータ，引数originには出発地点ポイントを表すspオブジェクト，3番目の引数destinには到着地点ポイントを表すspオブジェクトを指定します．引数moveでは，どの範囲のセルに移動できる設定にするか（今回は，縦横斜の8方向）を指定します．探索の結果得られたルートは，st_as_sf関数を用いて，sfオブジェクトに変換することができます．
We can search the routes with movecost function. We assign the elevation raster data to argument dtm, origin sp object to argument origin, and destination object to argument destin. We can select the number of directions in which cells are connected (8 directions in our case) with argument move. Obtained routes can be converted to sf objects with st_as_sf function.
#徒歩で移動する場合の最小距離（search least-cost hiking routes）
result_walk <- movecost::movecost(dtm=elev_s_viet_xy,origin=orig_point,destin=dest_point,move=8)
#ルートをsfオブジェクトに変換（convert routes to df objects）
path_walk <- sf::st_as_sf(result_walk$LCPs)