Simple feature | Rを使った地理空間データの可視化と分析

sfパッケージを用いてRでの空間データの取り扱いを学ぶ

Simple features あるいは simple feature access は、

現実世界の地物をコンピュータ上で表現する、また、そのようなオブジェクトがどのようにデータベースに格納・参照するための国際規格(ISO-19125-1:2004)です。

sfパッケージはsimple featureをRで扱うためのパッケージで、simple featureを作成するだけでなく、空間操作や多くの空間情報データを扱うファイルとの相互変換を可能にします。またPostGISからのデータ読み込みも可能なパッケージです。

library(sf)

sfパッケージは、GISでよく利用される、 GDALやGEOSをはじめとしたオープンソースのライブラリに依存します。そのためsfパッケージの利用にはこれらがあらかじめインストールされている必要があります。詳しくはsfのインストールについてをご覧ください。

sfパッケージの多くの関数はst_というプレフィックスがついています。これはspatial and temporalを意味しています。

simple feature geometry

主要なgeometryを紹介します。これらはGeoJSONでもサポートされているものです。ここからは実際にsfの関数を使ってsimple feature geometryの作成例を見ていきましょう。

種類	概略
`POINT`	単一の位置を示す
`LINESTRING`	2つ以上の位置の配列からなる線
`POLYGON`	2次元のジオメトリ。複数のポイントを結ぶ線で構成され、これらは互いに交差せずに閉じている。
`MULTIPOINT`	`POINT`の集合
`MULTILINESTRING`	`LINESTRING`の集合
`MULTIPOLYGON`	`POLYGON`の集合
`GEOMETRYCOLLECTION`	`GEOMETRYCOLLECTION`を除いたgeometryの組み合わせ

単一のポイントは数値ベクトル
ポイントを行列として表現することでLINESTRING, POLYGONを表現できる
他の組み合わせはリスト

ポイントを定義する関数はst_point()です。ポイントの位置を引数に与えます。引数には最低でも2つの数値ベクトルを与えなければいけません。st_point()に与えられる数値は、座標に変換されます。すなわち、2次元のジオメトリとして機能します。

# ポイントの作成
p1 <- 
  st_point(c(140.112, 36.083))

# 標準出力されます
p1

1
2
3
4
5
6

## POINT (140.112 36.083)

# ジオメトリの確認
p1 %>% 
  st_geometry_type()

1
2
3

## [1] POINT
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON ... TRIANGLE

1
2

さて、先ほどst_point()の実行に最低2つの数値ベクトルの入力が必要と述べました。ではそれ以上の数値を与えた場合はどうなるのでしょう。答えは3次元、4次元と次元をジオメトリの次元を増やす表現として機能します。

# st_point()にはポイントの座標を与えます
st_point(c(1, 2)) # 2つの値はXY座標として利用されます

1
2

## POINT (1 2)

st_point(c(1, 2, 3)) # XYに加え、Z座標を定義できます

## POINT Z (1 2 3)

st_point(c(1, 2, 3), dim = "XYM")

## POINT M (1 2 3)

# dim引数を利用し、明示的に次元とその値を指定できます
x <- 
  st_point(c(1, 2, 3, 4)) # 4次元の座標には4つの値を与えます

# Z/M座標を除外した座標値を返却します
st_zm(x)

1
2
3
4
5
6

## POINT (1 2)

set.seed(123)
mesh <- 
  jpmesh::rmesh(1) %>% 
  jpmesh::export_mesh()

st_bbox(mesh)

1
2
3
4
5
6

##      xmin      ymin      xmax      ymax 
## 129.96250  33.01667 129.97500  33.02500

1
2

st_area(mesh)

## 1078941 [m^2]

座標参照系の指定

sfgオブジェクトは座標を持ちますがあらかじめ座標参照系を含めることができません。

library(mapview) # 位置把握のために使います
st_crs(p1)

1
2

## Coordinate Reference System: NA

# mapview(p1)

st_sfc(p1, crs = 4326)

1
2
3

## Geometry set for 1 feature 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 140.112 ymin: 36.083 xmax: 140.112 ymax: 36.083
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

1
2
3
4
5
6

## POINT (140.112 36.083)

# st_sfc(p1, crs = 4326) %>%
#   mapview()

1
2

library(rnaturalearth)
library(ggplot2)

ne_jpn <- 
  ne_states(country = "Japan", returnclass = "sf") %>% 
  tibble::new_tibble(subclass = "sf")

ne_wld <- 
  ne_countries(returnclass = "sf") %>% 
  tibble::new_tibble(subclass = "sf")

plot(st_geometry(ne_wld))

ggplot() + 
  geom_sf(data = st_transform(ne_wld, 
                              crs = "+proj=laea +y_0=0 +lon_0=155 +lat_0=-90 +ellps=WGS84 +no_defs"))

st_crs(ne_jpn)

ne_jpn %>% 
  st_geometry() %>% 
  plot(graticule = st_crs(4326))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

sf_poi <- 
  tibble::data_frame(
  id = seq(1, 2),
  lng = c(140.112, 140.12),
  lat = c(36.083, 36.05)) %>% 
  st_as_sf(coords = c("lng", "lat"), 
           crs = 4326)

sf_poi

1
2
3
4
5
6
7
8
9

## Simple feature collection with 2 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 140.112 ymin: 36.05 xmax: 140.12 ymax: 36.083
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 2 x 2
##      id         geometry
##   <int>      <POINT [°]>
## 1     1 (140.112 36.083)
## 2     2   (140.12 36.05)

1
2
3
4
5
6
7
8
9
10
11

# sf_poi %>% 
#   mapview::mapview()

sf_poi %>% 
  st_coordinates()

1
2
3
4
5

##         X      Y
## 1 140.112 36.083
## 2 140.120 36.050

1
2
3

sf_poi %>% 
  dplyr::mutate(lng = st_coordinates(geometry)[, 1],
                lat = st_coordinates(geometry)[, 2]) %>% 
  st_set_geometry(NULL)

1
2
3
4

## # A tibble: 2 x 3
##      id   lng   lat
## * <int> <dbl> <dbl>
## 1     1  140.  36.1
## 2     2  140.  36.0

1
2
3
4
5

空間操作

ジオメトリを操作することで、2点間の距離を求めたり、ポイントにバッファを待たせたりという処理が可能になります。ジオメトリの操作には次のものが含まれます。対応するsfパッケージの関数名を併記します。

面積の計算... st_area()
重心(セントロイド)の抽出... st_centroid()
距離計算... st_distance()
属性テーブルの結合... st_join()
ディゾルブ... st_union()
バッファリング... st_buffer()
ボロノイ分割... st_voronoi()

library(units)

# 距離の算出には横メルカトル投影法の平面直角座標にしておきましょう
p1 <- 
  st_point(c(140.112, 36.083)) %>% 
  st_sfc(crs = 2451)
p2 <- 
  st_point(c(139.677, 35.661)) %>% 
  st_sfc(crs = 2451)

st_distance(p1, p2)

1
2
3
4
5
6
7
8
9

## Units: [m]
##           [,1]
## [1,] 0.6060602

1
2
3

単位はメートルですが、これをキロメートルあるいはマイルに変換するにはunitsパッケージの関数を利用します。

units::set_units(st_distance(p1, p2), 
                 km)

1
2

## Units: [km]
##              [,1]
## [1,] 0.0006060602

1
2
3

units::set_units(st_distance(p1, p2), 
                 mile)

1
2

## Units: [mile]
##              [,1]
## [1,] 0.0003765884

1
2
3

空間隣接行列

空間オブジェクト間の隣接関係を表す指標として「空間隣接行列」があります。

これはオブジェクトの近接関係

空間自己相関を扱う上で基礎となる理解となります。

当該メッシュに対して

隣接している = 1
隣接指定ない = 0

とすると隣接行列は

p1 <-
  st_point(c(0, 0))
p2 <-
  st_point(c(2, 2))

st_relate(p1, p2)

1
2
3
4
5
6

##      [,1]       
## [1,] "FF0FFF0F2"

1
2

3 * 3の格子状区画(グリッド)を考えます。

sfc <-
  st_sfc(st_point(c(0, 0)), 
         st_point(c(3, 3)))

sfc_grid <- 
  sfc %>% 
  st_make_grid(n = c(3, 3), 
               what = "polygons")

1
2
3
4
5
6
7
8

st_make_grid()は任意の区画数のグリッドまたは区画の重心点を作成するのに便利な関数です。

# plot(sfc_grid)

sf_grid <- 
  st_sf(sfc_grid) %>%
  rename(geometry = sfc_grid) %>% 
  mutate(id = paste0("x", 
                     c(seq(7, 9),
                       seq(4, 6),
                       seq(1, 3))))

# plot(sf_grid)
# invisible(
#   text(st_coordinates(st_centroid(sf_grid))[, 1], 
#        st_coordinates(st_centroid(sf_grid))[, 2], 
#        labels = 
#          as.character(sf_grid$id), 
#        cex = 2))

library(ggplot2)
ggplot() +
  geom_sf(data = sf_grid) +
  geom_text(data = sf_grid,
            aes(x = st_coordinates(st_centroid(geometry))[, 1],
                y = st_coordinates(st_centroid(geometry))[, 2],
                label = id)) +
  xlab(NULL) +
  ylab(NULL) +
  coord_sf(datum = NULL) +
  theme_void()


st_relate(sf_grid[5, ], 
          sf_grid[1, ])

st_relate(sf_grid[8, ], 
          sf_grid[1, ])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

隣接行列の表現方法としてエッジ、ルーク、クイーンがあります。

重ね合わせ

ポリゴンのオーバーレイ

論理積: ポリゴンxに対して別のポリゴンyを重ね合わせ、各オブジェクトが交差している部分を抽出して新しいポリゴンを作る

x <- 
  st_polygon(list(
  rbind(
  st_point(c(1, 1)),
  st_point(c(2, 1)),
  st_point(c(2, 2)),
  st_point(c(1, 2)),
  st_point(c(1 ,1))
)))


plot(x)

y <- 
  geom_rot2(x / 1.2, 45) + 0.75

ggplot() + 
  geom_sf(data  = x) +
  geom_sf(data = y) +
  coord_sf(datum = NULL) +
  theme_void()

x <- 
  x %>% 
  st_sfc() %>% 
  st_sf()

y <- 
  st_sf(y) %>% 
  st_cast("POLYGON")

st_join(x, y) %>% 
  plot()
st_join(x, y, join = st_disjoint) %>% 
  plot()

plot(x, col = "red")
plot(y, col = "blue", add = TRUE)

plot(st_union(x, y))
plot(st_union(x, y, by_feature = TRUE)) # 同じ?

st_join(x, y, join = st_overlaps) %>% 
  plot()

# 論理積
st_intersection(x, y) %>% 
  plot()
# 削除
st_difference(x, y) %>% 
  plot()
st_sym_difference(x, y) %>% 
  plot()
st_snap(x, y, tolerance = 0.5) %>% 
  plot()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

sfオブジェクトの作成

データフレームの列に地物の形状 geometry と属性 attributes を格納する
simple featureのための3つのクラスを提供する
- sf: simple featureとobjectの属性情報を紐づける
- sfc: simple featureのgeometry情報をlist-columnとして持つ
- sfg: 個別のsimple featureを扱うfeature geometry
- 大きさとしては sfg < sfc < sf
sfオブジェクトのgeometry列およびsfcオブジェクトは属性データをもたない。地物の形状 (geometry) を記録するだけ
- sfではgeometryをテキスト形式のwell-known text: (WKT)で記述する。

sfg: simple feature geometry

ポイント

sfg_pt <- 
  st_point(c(0, 0))

sfg_pt

1
2
3
4

## POINT (0 0)

st_geometry_type(sfg_pt)

## [1] POINT
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON ... TRIANGLE

1
2

# st_point()にはポイントの座標を与えます
st_point(c(1, 2)) # 2つの値はXY座標として利用されます

1
2

## POINT (1 2)

st_point(c(1, 2, 3)) # XYに加え、Z座標を定義できます

## POINT Z (1 2 3)

st_point(c(1, 2, 3), dim = "XYM") # dim引数を利用し、明示的に次元とその値を指定できます

## POINT M (1 2 3)

x <- 
  st_point(c(1, 2, 3, 4)) # 4次元の座標には4つの値を与えます

st_zm(x)

1
2
3
4

## POINT (1 2)

sfg_mpt <- 
  c(st_point(c(0, 0)),
    st_point(c(1, 1)))

sfg_mpt

1
2
3
4
5

## MULTIPOINT (0 0, 1 1)

sfgオブジェクトでは演算処理が可能です。

c(sfg_pt, sfg_pt + 1)

## MULTIPOINT (0 0, 1 1)

またsfgオブジェクトにはplot()が適用できます。これにより、視覚的にジオメトリの形状を把握することが可能になります。plot()の詳細は後述します。

plot(sfg_mpt)

plot of chunk unnamed-chunk-19

plot(c(sfg_pt, sfg_pt + 1), 
     col = "red", 
     pch = 16)

1
2
3

plot of chunk unnamed-chunk-19

all.equal(sfg_mpt,
          c(sfg_pt, sfg_pt + 1))

1
2

## [1] TRUE

ライン

sfg_line <- 
  st_linestring(sfg_mpt)

sfg_line

1
2
3
4

## LINESTRING (0 0, 1 1)

sfg_mline <- 
  c(sfg_line,
  sfg_line + 1)

sfg_mline

1
2
3
4
5

## MULTILINESTRING ((0 0, 1 1), (1 1, 2 2))

plot(sfg_mline, add = TRUE)

ポリゴン

# 起点（終点にもなる）を用意
# 起点と終点は同じ座標
sfg_pt <- 
  st_point(c(1, 1))

sfg_poly <- 
  rbind(
  sfg_pt, # 起点
  st_point(c(2, 1)),
  st_point(c(2, 2)),
  st_point(c(1, 2)),
  sfg_pt # 終点
  ) %>% 
  print() %>%
  list() %>% 
  st_polygon()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

##        [,1] [,2]
## sfg_pt    1    1
##           2    1
##           2    2
##           1    2
## sfg_pt    1    1

1
2
3
4
5
6

plot(sfg_poly)

# 閉じていないと怒られます
rbind(
  sfg_pt, # 起点
  st_point(c(2, 1)),
  st_point(c(2, 2)),
  st_point(c(1, 2))
  ) %>% 
  list() %>% 
  st_polygon()
# Error in MtrxSet(x, dim, type = "POLYGON", needClosed = TRUE) : 
#   polygons not (all) closed

1
2
3
4
5
6
7
8
9
10
11

sfg_poly_fault <- rbind(
  sfg_pt,
  st_point(c(2, 1)),
  st_point(c(2, 2)),
  st_point(c(0, 2)),
  st_point(c(1, 2)),
  sfg_pt) %>% 
  print() %>%
  list() %>% 
  st_polygon()

1
2
3
4
5
6
7
8
9
10

##        [,1] [,2]
## sfg_pt    1    1
##           2    1
##           2    2
##           0    2
##           1    2
## sfg_pt    1    1

1
2
3
4
5
6
7

st_is_valid(sfg_poly)

## [1] TRUE

st_is_valid(sfg_poly_fault)

## [1] FALSE

plot(sfg_poly_fault)

sfc: simple feature list column

sfcはsfgオブジェクトから構成されるオブジェクトです。st_sfc()にジオメトリを与えて作成します。引数にはジオメトリの他に、参照座標系を指定することができます。

複数のジオメトリをもつこともでき、その実態はジオメトリの種類、参照座標系、矩形の範囲、欠損の数といった要素をもつリストです。ジオメトリの形状と位置に関する情報を持っています。

st_sfc(sfg_pt)

## Geometry set for 1 feature 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 1 ymin: 1 xmax: 1 ymax: 1
## epsg (SRID):    NA
## proj4string:    NA

1
2
3
4
5
6

## POINT (1 1)

st_sfc(st_point(c(140.112, 36.083)), 
       crs = 4326)

1
2

## Geometry set for 1 feature 
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 140.112 ymin: 36.083 xmax: 140.112 ymax: 36.083
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

1
2
3
4
5
6

## POINT (140.112 36.083)

sf: simple feature

sfcに属性情報を追加したものがsfオブジェクトとして扱われます。st_sf()がsfオブジェクトを作成する関数となります。

pt1 <- st_point(c(0, 1))
pt2 <- st_point(c(1, 1))
sfc <- st_sfc(pt1, pt2)

sfc %>% 
  st_sf()

1
2
3
4
5
6

## Simple feature collection with 2 features and 0 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 0 ymin: 1 xmax: 1 ymax: 1
## epsg (SRID):    NA
## proj4string:    NA
##      geometry
## 1 POINT (0 1)
## 2 POINT (1 1)

1
2
3
4
5
6
7
8
9

# a という属性情報を付与したsfオブジェクトを作成
data.frame(a = 1:2, 
           geometry = sfc) %>% 
  st_sf()

1
2
3
4

## Simple feature collection with 2 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 0 ymin: 1 xmax: 1 ymax: 1
## epsg (SRID):    NA
## proj4string:    NA
##   a    geometry
## 1 1 POINT (0 1)
## 2 2 POINT (1 1)

1
2
3
4
5
6
7
8
9

読み込み

sfでの地理空間データの読み込み関数としてst_read()があります。この関数を使うことで、地物情報を記録したファイルやPostGISデータベース上のデータをRで操作可能になります。

日本語の文字化け

一部のファイルでは日本語等のマルチバイト文字列が「文字化け」を起こしていることがあります。その場合、正しく読み込むために、次のようにoptions引数を利用してファイルのエンコード形式を変更する必要があります。

# ENCODING=CP932 の中にスペースを入れない
st_read(dsn, options = c("ENCODING=CP932"))

1
2

データベース

他のクラスとsfへの相互変換

csvファイル等に座標の値が列に記録されている場合、データフレームとして読み込んでおくことでsfオブジェクトとして扱うことができます。これには任意の列を座標として扱うst_as_sf()を利用します。

sf_poi <- 
  tibble::data_frame(
  id = seq(1, 2),
  lng = c(140.112, 140.12),
  lat = c(36.083, 36.05)) %>% 
  st_as_sf(coords = c("lng", "lat"),
           crs = 4326)

sf_poi

1
2
3
4
5
6
7
8
9

## Simple feature collection with 2 features and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 140.112 ymin: 36.05 xmax: 140.12 ymax: 36.083
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 2 x 2
##      id         geometry
##   <int>      <POINT [°]>
## 1     1 (140.112 36.083)
## 2     2   (140.12 36.05)

1
2
3
4
5
6
7
8
9
10
11

また、sfオブジェクトをデータフレームとして扱うにはst_set_geometry()関数にNULLを指定してジオメトリ情報を除外します。

sf_poi %>% 
  st_set_geometry(value = NULL)

1
2

## # A tibble: 2 x 1
##      id
## * <int>
## 1     1
## 2     2

1
2
3
4
5

一方でこの方法では座標情報が失われてしまうので、列として座標値を残したい場合には次のようにする必要があります。st_coordinate()はジオメトリの座標を求める関数です。ジオメトリがポイントの場合はXY座標を返します。対象のジオメトリがポリゴンの場合にはst_centroid()を使い、重心の座標に変換しておく等の処理が必要になります。

sf_poi %>% 
  dplyr::mutate(lng = st_coordinates(geometry)[, 1],
                lat = st_coordinates(geometry)[, 2]) %>% 
  st_set_geometry(NULL)

1
2
3
4

## # A tibble: 2 x 3
##      id   lng   lat
## * <int> <dbl> <dbl>
## 1     1  140.  36.1
## 2     2  140.  36.0

1
2
3
4
5

またsfでは、従来Rで地理空間情報データ（特にベクトルデータ）を扱うのに用いられているspパッケージのオブジェクトとの相互互換のための関数が用意されています。変換は次のように行います。

library(spatstat)
data(chicago)
class(chicago)

1
2
3

## [1] "lpp" "ppx"

st_as_sf(chicago)

## Simple feature collection with 620 features and 4 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: 0.3893523 ymin: 153.1034 xmax: 1281.986 ymax: 1276.56
## epsg (SRID):    NA
## proj4string:    NA
## First 10 features:
##      label seg tp marks                           geom
## 1   window  NA NA  <NA> POLYGON ((1281.986 153.1035...
## 2  segment  NA NA  <NA> LINESTRING (0.3894739 1253....
## 3  segment  NA NA  <NA> LINESTRING (109.683 1251.77...
## 4  segment  NA NA  <NA> LINESTRING (109.683 1251.77...
## 5  segment  NA NA  <NA> LINESTRING (198.1486 1276.5...
## 6  segment  NA NA  <NA> LINESTRING (197.9988 1251.1...
## 7  segment  NA NA  <NA> LINESTRING (290.4787 1276.5...
## 8  segment  NA NA  <NA> LINESTRING (288.9907 1250.5...
## 9  segment  NA NA  <NA> LINESTRING (380.1326 1276.5...
## 10 segment  NA NA  <NA> LINESTRING (379.9827 1249.8...

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

as_Spatial(sf_poi)

## class       : SpatialPointsDataFrame 
## features    : 2 
## extent      : 140.112, 140.12, 36.05, 36.083  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : id 
## min values  :  1 
## max values  :  2

1
2
3
4
5
6
7
8

空間解析パッケージはsfをベースにしたものが増えていますが、現状ではまだまだspベースのパッケージが多いです。用途に応じてオブジェクトの変換を行うのが良いでしょう。

# sfパッケージを用いてRでの空間データの取り扱いを学ぶ

# simple feature geometry

# 座標参照系の指定

# 空間操作

# 空間隣接行列

# 重ね合わせ

# sfオブジェクトの作成

# sfg: simple feature geometry

# ポイント

# ライン

# ポリゴン

# sfc: simple feature list column

# sf: simple feature

# 読み込み

# 日本語の文字化け

# データベース

# 他のクラスとsfへの相互変換

sfパッケージを用いてRでの空間データの取り扱いを学ぶ

simple feature geometry

座標参照系の指定

空間操作

空間隣接行列

重ね合わせ

sfオブジェクトの作成

sfg: simple feature geometry

ポイント

ライン

ポリゴン

sfc: simple feature list column

sf: simple feature

読み込み

日本語の文字化け

データベース

他のクラスとsfへの相互変換