RでGIS
GISというとArcGISもしくはQGISを想定しがちだけど、個人的にはこれらGUIに依拠したSoftware（マウスでカチカチする系）は好みではないし、おすすめしない。というのも、何度もマウスでクリックしながら作業を進めるので、以下のリスクがある。
作業行程の再現性を担保しにくい
エラーがあった場合、すべての行程を繰り返す必要がある
うっかりミスしやすい
これらの理由から、私はGISもスクリプトベースで作業行程を記録すべきと思う。スクリプトベースのGISというとPostGISが標準だったと思うが、最近ではこれらの関数群はほぼすべてRで使えるようになっているし、工夫すれば計算速度もArcGISと遜色ない。しかし、スクリプトベースのGISに関する資料は英語によるものがほとんどのため、アクセスしにくいといえばそうかもしれない。そこで、RのGIS処理について、いくつかの記事にわけて紹介したいと思う。なお、このブログで書かれていることの大半はTaro Mienoさん(University of Nebraska Lincoln)のオンライン資料で勉強させていただいた。とても事細かな説明があるので、この記事を読んで興味がわいた方はぜひこちらで勉強するといいと思う。
GIS用のRパッケージ
どのようなタイプのレイヤー（ベクター、ラスター）を扱うかによって、必要となるパッケージは変わる。基本的な作業はsf, raster, terraあたりで済むが、多少込み入った作業には他の補助的なパッケージも必要になる。以下の表に簡単にまとめる。
| パッケージ名 | レイヤタイプ | 備考 |
|---|---|---|
sf |
ベクター | ベクター処理の場合はsf一択。大体これで足りる（Website）。 |
rmapshaper |
ベクター | フィーチャーが重いときに頂点を削って計算を高速化する。（Website） |
raster |
ラスター | ラスター処理のデファクトスタンダードだったが、C++ベースのterraに移行中（Website）。 |
terra |
ラスター | ラスター処理の新しいデファクトスタンダードになりつつある（Website）。 |
stars |
ベクター、ラスター | 時空間情報のついたラスター処理に特化。ラスター｜ベクター変換の際にも重宝する（Website）。 |
exactextractr |
ベクター、ラスター | フィーチャーごとにラスターデータを集計するときに使う（Website）。 |
tmap |
NA | GISレイヤーの描画（Website） |
sfによるベクター処理
データを読み込む
ベクター処理の鉄板であるsfを使ってみる。sfについてくるデータセットを使う。拡張子を見てもらうとわかるが、呼び出されるファイルはシェイプファイル（.shp）。
library(sf)
fname <- system.file("shape/nc.shp", package="sf")
print(fname)## [1] "C:/Users/Akira/OneDrive - UNCG/Documents/R/win-library/4.0/sf/shape/nc.shp"このファイルパスをst_read()で読み込む。ここでは例を示すためにsf内のデータセットを呼び出しているが、read_csv()のような感じでシェイプファイルを読める（詳しくはこちら）。プリントしてみると、何やらデータフレームのようなものが吐き出されるが、これはGISでいう属性テーブルに対応している。フィーチャーについて、座標系、Geometry type、CoordinateなどGIS作業になくてはならない情報もついている。
sf_nc <- st_read(fname)## Reading layer `nc' from data source
## `C:\Users\Akira\OneDrive - UNCG\Documents\R\win-library\4.0\sf\shape\nc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27print(sf_nc)## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
## First 10 features:
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## 3 0.143 1.630 1828 1828 Surry 37171 37171 86 3188 5
## 4 0.070 2.968 1831 1831 Currituck 37053 37053 27 508 1
## 5 0.153 2.206 1832 1832 Northampton 37131 37131 66 1421 9
## 6 0.097 1.670 1833 1833 Hertford 37091 37091 46 1452 7
## 7 0.062 1.547 1834 1834 Camden 37029 37029 15 286 0
## 8 0.091 1.284 1835 1835 Gates 37073 37073 37 420 0
## 9 0.118 1.421 1836 1836 Warren 37185 37185 93 968 4
## 10 0.124 1.428 1837 1837 Stokes 37169 37169 85 1612 1
## NWBIR74 BIR79 SID79 NWBIR79 geometry
## 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
## 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
## 3 208 3616 6 260 MULTIPOLYGON (((-80.45634 3...
## 4 123 830 2 145 MULTIPOLYGON (((-76.00897 3...
## 5 1066 1606 3 1197 MULTIPOLYGON (((-77.21767 3...
## 6 954 1838 5 1237 MULTIPOLYGON (((-76.74506 3...
## 7 115 350 2 139 MULTIPOLYGON (((-76.00897 3...
## 8 254 594 2 371 MULTIPOLYGON (((-76.56251 3...
## 9 748 1190 2 844 MULTIPOLYGON (((-78.30876 3...
## 10 160 2038 5 176 MULTIPOLYGON (((-80.02567 3...これだけみてもよくわらからないので、図示する。tmapパッケージの関数群を使えば、ggplotのように地図をRで書ける。どうやらこれはアメリカ、ノースカロライナ州のCounty（行政区のようなもの）のポリゴンのようだ。
library(tmap)
tmap_mode("view")
tm_shape(sf_nc) +
tm_polygons()フィールド編集
sfを使う最大のメリットの一つは、レイヤーをそのままR上で編集できることだろう。しかもtidyverseとの相性が最高で、だいたいのdplyrの関数はそのまま使える。ここでは練習のため、もともと入っていた列情報を一旦削除し、逐一情報を追加する。
# remove all columns except geometry
sf_nc0 <- sf_nc %>%
dplyr::select(NULL)
print(sf_nc0)## Simple feature collection with 100 features and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
## First 10 features:
## geometry
## 1 MULTIPOLYGON (((-81.47276 3...
## 2 MULTIPOLYGON (((-81.23989 3...
## 3 MULTIPOLYGON (((-80.45634 3...
## 4 MULTIPOLYGON (((-76.00897 3...
## 5 MULTIPOLYGON (((-77.21767 3...
## 6 MULTIPOLYGON (((-76.74506 3...
## 7 MULTIPOLYGON (((-76.00897 3...
## 8 MULTIPOLYGON (((-76.56251 3...
## 9 MULTIPOLYGON (((-78.30876 3...
## 10 MULTIPOLYGON (((-80.02567 3...まず、投影変換して投影座標系に直し、それぞれのポリゴンの面積を計算する。その後、面積順にID番号を与えてみよう。投影変換はst_transform() 、面積はst_area() で計算できる（そのほかにも長さst_length()、重心st_centroid()、バッファーst_buffer()なども計算できる。関数一覧はここ）。アウトプットを見ると、座標系が変換され、新しい列も追加されていることがわかる。
sf_nc1 <- sf_nc0 %>%
st_transform(crs = 32617) %>% # EPSG: 32617 = WGS84/UTM Zone 17
mutate(area = st_area(.)) %>% # calculate area
arrange(desc(area)) %>% # arrange by area; decreasing order
mutate(area_rank = row_number()) # assign row numbers as it appears
print(sf_nc1)## Simple feature collection with 100 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 196602 ymin: 3751723 xmax: 1002267 ymax: 4057841
## Projected CRS: WGS 84 / UTM zone 17N
## First 10 features:
## area geometry area_rank
## 1 2451191239 [m^2] MULTIPOLYGON (((716665.9 37... 1
## 2 2439491560 [m^2] MULTIPOLYGON (((696137.5 38... 2
## 3 2439475275 [m^2] MULTIPOLYGON (((764328.5 38... 3
## 4 2289719021 [m^2] MULTIPOLYGON (((751778.2 38... 4
## 5 2195301272 [m^2] MULTIPOLYGON (((688395.9 39... 5
## 6 2183764603 [m^2] MULTIPOLYGON (((773656.8 38... 6
## 7 2167579313 [m^2] MULTIPOLYGON (((716665.9 37... 7
## 8 2083418432 [m^2] MULTIPOLYGON (((723776.1 39... 8
## 9 2074226184 [m^2] MULTIPOLYGON (((803168.3 38... 9
## 10 2051040434 [m^2] MULTIPOLYGON (((853257.9 39... 10このままだとarea の列がわかりにくいので、単位をkm\(^2\)に変換する。
sf_nc1 <- sf_nc1 %>%
mutate(area = units::set_units(area, "km^2"))
print(sf_nc1)## Simple feature collection with 100 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 196602 ymin: 3751723 xmax: 1002267 ymax: 4057841
## Projected CRS: WGS 84 / UTM zone 17N
## First 10 features:
## area geometry area_rank
## 1 2451.191 [km^2] MULTIPOLYGON (((716665.9 37... 1
## 2 2439.492 [km^2] MULTIPOLYGON (((696137.5 38... 2
## 3 2439.475 [km^2] MULTIPOLYGON (((764328.5 38... 3
## 4 2289.719 [km^2] MULTIPOLYGON (((751778.2 38... 4
## 5 2195.301 [km^2] MULTIPOLYGON (((688395.9 39... 5
## 6 2183.765 [km^2] MULTIPOLYGON (((773656.8 38... 6
## 7 2167.579 [km^2] MULTIPOLYGON (((716665.9 37... 7
## 8 2083.418 [km^2] MULTIPOLYGON (((723776.1 39... 8
## 9 2074.226 [km^2] MULTIPOLYGON (((803168.3 38... 9
## 10 2051.040 [km^2] MULTIPOLYGON (((853257.9 39... 10この情報をggplotにそのまま引き渡し、図示もできる。
sf_nc1 %>%
mutate(area = as.numeric(area)) %>% # remove unit [km^2]
ggplot(aes(x = area)) +
geom_histogram()
あるいはfilter()に引き渡して、特定の特徴をもつポリゴンだけを抽出することも可能。
sf_nc1 %>%
dplyr::filter(area > units::set_units(1000, "km^2")) %>% # select polygons with > 1000 km^2
tm_shape() +
tm_polygons() # show mapレイヤー間の紐づけ
どのポイントがどのポリゴンに属しているのか、という情報を抽出したい場合も多い。この場合、st_join()が有効。先のポリゴンレイヤーの範囲内に、ランダムなポイントデータを生成する。
# generate random points ----
sf_point <- tibble() %>%
bind_rows(c(st_bbox(sf_nc1))) %>% # extract the boundary box (bbox) with "st_bbox()"
summarize(Y = runif(100, ymin, ymax),
X = runif(100, xmin, xmax)) %>% # random points within the bbox
mutate(point_id = row_number()) %>%
st_as_sf(coords = c("X", "Y")) %>% # define as sf object
st_set_crs(st_crs(sf_nc1)) # define CRS
# show map ----
tmap_mode("view")
tm_shape(shp = sf_nc1) +
tm_polygons() +
tm_shape(shp = sf_point) +
tm_dots(col = "black")ポイントをポリゴンと紐づける。x で指定されるレイヤーに、yのレイヤーの情報を紐づける。紐づけの条件はいろいろ変えれるが、デフォルトでは両者のフィーチャが重なっている場合に紐づけされる。この例では、ポリゴンの中に含まれたポイントには情報が付加されたが、ポリゴン外にあるものについてはNAが充てられている。
sf_point_join <- st_join(x = sf_point, y = sf_nc1)
print(sf_point_join)## Simple feature collection with 100 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 198793.6 ymin: 3753806 xmax: 1000434 ymax: 4057215
## Projected CRS: WGS 84 / UTM zone 17N
## # A tibble: 100 × 4
## point_id geometry area area_rank
## * <int> <POINT [m]> [km^2] <int>
## 1 1 (256849.5 3848752) NA NA
## 2 2 (929876 3964631) 992. 68
## 3 3 (739275.8 3974891) 2195. 5
## 4 4 (384740.5 3900150) 601. 95
## 5 5 (876817.5 3768223) NA NA
## 6 6 (880047.9 3895897) 836. 78
## 7 7 (635296.7 3972869) 1110. 59
## 8 8 (868596.4 3953171) 2051. 10
## 9 9 (415847.7 3960724) 1340. 40
## 10 10 (956227.4 3938149) NA NA
## # … with 90 more rowsareaの列がNAになっているポイントを取り除き、再度図示する。途中、drop_na()を使うとsfオブジェクトがtibbleに強制変換されるようなので、再度st_as_sf()を使ってsfに戻している。
# filter points ----
sf_point_filter <- sf_point_join %>%
drop_na(area) %>% # forced to be a "tibble"
st_as_sf() # back convert to sf
# view map ----
tmap_mode("view")
tm_shape(sf_nc1) +
tm_polygons() +
tm_shape(sf_point_filter) +
tm_dots(col = "red")まとめ
sfの使い方についてごく初歩的な部分に触れたが、本当にいろんなことができるので、ぜひsfパッケージや他のページを参照してみてほしい。