Polygon 면적 구하기: sf 와 raster 패키지

shapefile에 담겨져 있는 polygon의 면적을 구해보자 raster 패키지의 area 혹은 sf 패키지의 st_area 함수를 이용할 수 있다.

library(sf)
kor <- st_read("gadm41_KOR_1.shp")

Reading layer `gadm41_KOR_1' from data source 
  `C:\Users\jonghoon.kim\Documents\myblog\blog\posts\area-polygon\gadm41_KOR_1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 17 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 124.6097 ymin: 33.11236 xmax: 131.8715 ymax: 38.61177
Geodetic CRS:  WGS 84

# another way
# kor <- read_sf(dsn="C:/Users/jonghoon.kim/Documents/myblog/posts/area-polygon", layer="gadm41_KOR_1")
set.seed(42)
# 1e6 to get population density per 1 km^2
kor$area_sqkm1 <- as.numeric(st_area(kor))/1e6

# raster package way
# library(raster)
# kor <- shapefile("C:/Users/jonghoon.kim/Documents/myblog/posts/area-polygon/gadm41_KOR_1.shp")
# raster pkg works for the Spatial* object
# kor$area_sqkm2 <- raster::area(as(kor, 'Spatial'))/1e6

# plot population density
library(ggplot2)
theme_set(hrbrthemes::theme_ipsum_rc(base_size=14, subtitle_size=16, axis_title_size=12))
extrafont::loadfonts()

labels = expression(atop("Population density", per~1~km^2))
plt <- ggplot(kor)+
  geom_sf(aes(fill=area_sqkm1))+
  scale_fill_viridis_c(name=labels) +
  theme(legend.position="right")

# ggsave("southkorea_map.png", gg, units="in", width=3.4*2, height=2.7*2)  
plt

# use RColorBrewer 
library(RColorBrewer)
pal <- brewer.pal(9,"YlOrBr")
plt + scale_fill_gradientn(name=labels, colors=pal)