Split the map into a grid using sf

sf 패키지를 이용해서 지도를 grid로 나누고 특정 좌표가 몇번째 grid에 있는지 알아보기

sf
Author

Don Don

Published

March 9, 2024

library(sp)
library(sf)
library(tidyverse)
library(data.table)
library(plotly)
theme_set(theme_bw())

따릉이 대여소 정보 데이터 불러오기

서울시 공공 자전거 대여소 정보는 https://data.seoul.go.kr/dataList/OA-13252/F/1/datasetView.do 에서 다운로드가 가능하다.

station <- read_csv('~/Desktop/quarto_blog/posts/data/공공자전거 대여소 정보.csv', locale = locale(encoding = 'CP949'))
colnames(station) <- c('대여소번호', '대여소명', '자치구', '상세주소', 'lat', 'long', '설치시기', 'LCD거치대수', 'QR거치대수', '운영방식')
station %>% head()
# A tibble: 6 × 10
  대여소번호 대여소명         자치구 상세주소   lat  long 설치시기   LCD거치대수
       <dbl> <chr>            <chr>  <chr>    <dbl> <dbl> <date>           <dbl>
1       2301 현대고등학교 건… 강남구 서울특…   37.5  127. 2017-06-13          10
2       2302 교보타워 버스정… 강남구 서울특…   37.5  127. 2017-06-13          10
3       2303 논현역 7번출구   강남구 서울특…   37.5  127. 2017-06-13          15
4       2304 신영 ROYAL PALA… 강남구 서울특…   37.5  127. 2017-06-13          10
5       2305 MCM 본사 직영점… 강남구 서울특…   37.5  127. 2017-06-13          10
6       2306 압구정역 2번 출… 강남구 서울특…   37.5  127. 2017-06-13          30
# ℹ 2 more variables: QR거치대수 <dbl>, 운영방식 <chr>
station_lcd <- station %>% filter(운영방식 == 'LCD') %>% select(대여소번호, 자치구, lat, long)
station_QR <- station %>% filter(운영방식 == 'QR') %>% select(대여소번호, 자치구, lat, long)

좌표계 변환

이전 포스팅 참고 https://statistics-dondon-7781aa.netlify.app/2021/01/16/map-visualization-using-sf-package/

temp <- station_lcd

colnames(temp)[c(3,4)] = c("Y","X")
convertCoordSystem <- function(data = temp, from.crs = from.crs, to.crs = to.crs){
  xy = data.frame(long = data$X, lat = data$Y)
  xy = as.data.frame(apply(xy, 2, as.numeric))
  
  coordinates(xy) = ~long+lat
  
  from.crs = CRS(from.crs)
  from.coordinates = SpatialPoints(xy, proj4string=from.crs)
  
  to.crs = CRS(to.crs)
  changed = as.data.frame(SpatialPoints(spTransform(from.coordinates, to.crs)))
  names(changed) = c("long", "lat")
  
  return(changed)
}

from.crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
to.crs = "+proj=tmerc +lat_0=38 +lon_0=127.5 +k=0.9996 +x_0=1000000 +y_0=2000000 +ellps=GRS80 +units=m +no_defs"
coord = convertCoordSystem(data = temp, from.crs = from.crs, to.crs = to.crs)
temp[,c(3,4)] = coord 
colnames(temp)[c(3,4)] = c("X","Y")
coordinates(temp) = ~X+Y
temp = st_as_sf(temp , coord = c("X", "Y") , crs = 5179)

서울 경계구역 데이터 불러오기

대한민국 경계구역 데이터는 이전 블로그 포스팅을 참고하면 불러올 수 있다. https://statistics-dondon-7781aa.netlify.app/2021/01/15/map-visualization/

대한민국 경계 데이터에서 서울 데이터만 추출하였다. shp 파일에 한글이 섞여있는데 icov 함수를 이용해서 인코딩을 UTF-8로 변경하면 에러 없이 불러올 수 있다.

nc <- st_read('~/Desktop/quarto_blog/posts/data/CTPRVN_202005/CTPRVN.shp') %>% st_transform()
Reading layer `CTPRVN' from data source 
  `/Users/sangdon/Desktop/quarto_blog/posts/data/CTPRVN_202005/CTPRVN.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 17 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 746110.3 ymin: 1458754 xmax: 1387950 ymax: 2068444
Projected CRS: PCS_ITRF2000_TM
nc$CTP_KOR_NM <- iconv(nc$CTP_KOR_NM, from = "CP949", to = "UTF-8", sub = NA, mark = TRUE, toRaw = FALSE)
nc <- nc %>% filter(str_detect(CTP_ENG_NM, 'Seoul'))

CRS 값 변경

서로 다른 데이터를 시각화할 때 좌표계가 일치하지 않으면 지도에 시각화를 할 수 없다. 따라서 반드시 좌표계를 일치시켜야만 같은 지도에 올바른 위치에 점이 찍히게 된다. st_crs 함수를 이용해서 따릉이 대여소 데이터의 좌표계와 서울 경계구역 데이터의 좌표계를 일치시켜주었다.

st_crs(temp) = st_crs(nc)

Grid 개수 지정

지도를 몇 개의 grid로 나눌건지를 설정해야한다. 나는 서울을 20×10=200 grid로 나누었다.

pts <- temp %>% st_sf

grid_50 <- st_make_grid(nc, n = c(20, 10)) %>% 
  st_sf(grid_id = 1:length(.))
pts
Simple feature collection with 1531 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 938070.3 ymin: 1937734 xmax: 971761.8 ymax: 1965672
Projected CRS: PCS_ITRF2000_TM
# A tibble: 1,531 × 3
   대여소번호 자치구           geometry
 *      <dbl> <chr>         <POINT [m]>
 1       2301 강남구 (957744.4 1947304)
 2       2302 강남구 (957953.8 1945252)
 3       2303 강남구 (957709.7 1945912)
 4       2304 강남구 (958979.2 1946017)
 5       2305 강남구 (958866.3 1946918)
 6       2306 강남구 (958358.2 1947640)
 7       2307 강남구 (959228.7 1947801)
 8       2308 강남구 (958967.5 1947878)
 9       2309 강남구   (960180 1946719)
10       2310 강남구 (959368.8 1947245)
# ℹ 1,521 more rows

각 grid에 대한 label 생성

grid_lab <- st_centroid(grid_50) %>% cbind(st_coordinates(.))
Warning: st_centroid assumes attributes are constant over geometries
grid_lab
Simple feature collection with 200 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 935961.1 ymin: 1938182 xmax: 971141.8 ymax: 1965471
Projected CRS: PCS_ITRF2000_TM
First 10 features:
   grid_id        X       Y                        .
1        1 935961.1 1938182 POINT (935961.1 1938182)
2        2 937812.7 1938182 POINT (937812.7 1938182)
3        3 939664.3 1938182 POINT (939664.3 1938182)
4        4 941515.9 1938182 POINT (941515.9 1938182)
5        5 943367.5 1938182 POINT (943367.5 1938182)
6        6 945219.1 1938182 POINT (945219.1 1938182)
7        7 947070.8 1938182 POINT (947070.8 1938182)
8        8 948922.4 1938182 POINT (948922.4 1938182)
9        9 950774.0 1938182   POINT (950774 1938182)
10      10 952625.6 1938182 POINT (952625.6 1938182)

서울시 경계구역과 grid 시각화

ggplot() +
  geom_sf(data = nc, fill = 'white', lwd = 0.05) +
  geom_sf(data = pts, color = 'red', size = 1.7) + 
  geom_sf(data = grid_50, fill = 'transparent', lwd = 0.3) +
  geom_text(data = grid_lab, aes(x = X, y = Y, label = grid_id), size = 2) +
  coord_sf(datum = NA)  +
  labs(x = "") +
  labs(y = "")

따릉이 대여소가 몇번째 grid에 해당하는지 찾기

pts %>% st_join(grid_50, join = st_intersects) %>% as.data.frame %>% head()
  대여소번호 자치구                 geometry grid_id
1       2301 강남구 POINT (957744.4 1947304)      73
2       2302 강남구 POINT (957953.8 1945252)      53
3       2303 강남구 POINT (957709.7 1945912)      73
4       2304 강남구 POINT (958979.2 1946017)      73
5       2305 강남구 POINT (958866.3 1946918)      73
6       2306 강남구 POINT (958358.2 1947640)      73

Citation

BibTeX citation:
@online{don2024,
  author = {Don, Don and Don, Don},
  title = {Split the Map into a Grid Using Sf},
  date = {2024-03-09},
  url = {https://dondonkim.netlify.app/posts/2021-06-12-split-the-map-into-a-grid-using-sf/sf2.html},
  langid = {en}
}
For attribution, please cite this work as:
Don, Don, and Don Don. 2024. “Split the Map into a Grid Using Sf.” March 9, 2024. https://dondonkim.netlify.app/posts/2021-06-12-split-the-map-into-a-grid-using-sf/sf2.html.