kh-blog

今話題のsfパッケージ、まさか知らない人はいないでしょう??。
今回はじめて使いましたので、簡単にメモ。

library(dplyr)

## Warning: パッケージ 'dplyr' はバージョン 3.4.2 の R の下で造られました

## 
##  次のパッケージを付け加えます: 'dplyr'

##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag

##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union

library(sf)

## Warning: パッケージ 'sf' はバージョン 3.4.4 の R の下で造られました

## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3

Shapeファイルからの読み込みの例が多かったので、
data.frameから作る場合をやってみます。 irisに適当に座標をつけます。
この場合はst_as_sf()で簡単に作ることができます。 SP*オブジェクトよりもスッキリしてる気がします。

iris_sf_wgs <- iris %>% 
  mutate(x = rnorm(dim(iris)[1],140,1),
         y = rnorm(dim(iris)[1],35,1)) %>% 
  st_as_sf(coords = c("x","y"),
           crs = 4326)
iris_sf_wgs %>% 
  head()

## Simple feature collection with 6 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 138.3046 ymin: 33.53104 xmax: 141.2912 ymax: 35.23003
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
##                    geometry
## 1 POINT (140.3085 33.53104)
## 2  POINT (140.6285 34.3082)
## 3 POINT (139.1832 34.60905)
## 4 POINT (138.3046 34.79331)
## 5 POINT (141.2912 34.66879)
## 6 POINT (139.8825 35.23003)

geometryに座標が入りました。postGISみたいなやつですね。 投影変換はst_transformです。crsに直接EPSGコードをぶち込めます。

iris_sf_wgs %>% 
  st_transform(crs = 32654) %>% 
  head()

## Simple feature collection with 6 features and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 253387.4 ymin: 3710374 xmax: 526682.7 ymax: 3899125
## epsg (SRID):    32654
## proj4string:    +proj=utm +zone=54 +datum=WGS84 +units=m +no_defs
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
##                   geometry
## 1 POINT (435796.6 3710374)
## 2 POINT (465816.2 3796392)
## 3 POINT (333421.1 3831190)
## 4 POINT (253387.4 3853434)
## 5 POINT (526682.7 3836353)
## 6 POINT (398313.7 3899125)

plotでの描画はこの通り。

plot(iris_sf_wgs %>% select(Species))

ggplotで扱うにはまだ開発中であるgeom_sfがあります。

Happy mapping!