R, Python 분석과 프로그래밍의 친구 (by R Friend)

지리적 레스터 데이터 (geographic raster data) 는 일반적으로 (a) 레스터 헤더(raster header)와 (b) 행렬 (matrix)로 구성이 됩니다.

- (a) 레스터 헤더 (raster header): 좌표 참조 시스템(CRS, Coordinate Reference System), 시작점(the origin)과 범위 (the extent)를 정의함
- (b) 행렬 (matrix): 동일한 크기의 픽셀 또는 셀(pixel, or cell)을 표현. 픽셀 ID(pixel IDs)와 픽셀 값(pixel values).

보통 행렬의 시작점(the origin, or starting point)은 행렬의 좌측 하단 구석에 위치한 좌표를 의미하는데요, R의 raster 패키지의 레스터 헤더에서는 시작점의 기본값으로 좌측 상단 구석에 위치한 좌표를 시작점으로 사용하므로 주의가 필요합니다.

레스터 헤더에서 범위 (the extent) sms 행의 개수, 열의 개수, 셀의 크기 해상도로 정의합니다. 각각의 셀에 접근하거나 수정하려면 시작점(the origin)으로 부터의 셀 ID를 사용하거나 또는 명시적으로 행과 열을 지정하면 됩니다.

레스터 데이터의 행렬 표현법은 각 셀의 네 개 구석의 좌표를 명시적으로 저장하지 않으며, 대신 시작점(the origin)만 저장하고 나머지는 시작점으로부터의 행과 열의 ID를 가지고 각 셀에 접근하는 방식이므로 벡터 데이터 표현과 비교해서 상대적으로 효율적이고 속도가 빠릅니다. (예: 벡터 사각형 폴리곤의 경우 각 셀별로 5개 점의 좌표를 저장해야 함.) 하지만 벡터 데이터의 각 도형별로 여러개의 값을 가질 수 있는 반면에, 레스터 데이터의 경우 각 셀별로 하나의 값만을 가질 수 있습니다.

[레스터 데이터 유형 (1) 셀 ID (Cell IDs), (2) 셀 값 (Cell Values) ]

지리공간 벡터 데이터 처리 및 분석에 sf 패키지를 사용했었는데요, 레스터 데이터 처리 및 분석은 R raster 패키지를 사용합니다. 

spDataLarge 패키지에 내장되어 있는 미국 유타지역의 Zion 국립공원 지역의 레스터 샘플 데이터를 raster() 함수를 사용해서 불러오고, 레스터 데이터의 속성정보들을 살펴보겠습니다. 혹시 raster, spDataLarge, rgdal 패키지를 사전에 설치하지 않았다면 install.packages() 로 패키지를 먼저 설치해주세요.

raster() 함수로 불러온 레스터 데이터셋 이름 'srtm_raster' 을 입력해주면 클래스(class) 차원(dimentions), 해상도(resolution), 범위(extent), 좌표 참조 시스템 (Coordinates Reference System), 출처(Source), 이름(names), 최소/최대 값(min, max values) 속성 정보를 볼 수 있습니다.

# =================================
# R GeoSpatial Data Analysis 
# Raster Data using raster package
# =================================

library(raster)
library(spDataLarge)

install.packages("rgdal")
library(rgdal)

## raster dataset from the spDataLarge package, 
## a few raster objects and one vector object covering an area of the Zion National Park (Utah, USA).
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
srtm_raster = raster(raster_filepath)
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
srtm_raster
# class      : RasterLayer 
# dimensions : 457, 465, 212505  (nrow, ncol, ncell)
# resolution : 0.0008333333, 0.0008333333  (x, y)
# extent     : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/srtm.tif 
# names      : srtm 
# values     : 1024, 2892  (min, max)

이들 레스터 데이테의 속성 정보에 함수를 사용해서 일일이 접근이 가능합니다. 먼저, 좌표 참조 시스템(Coordinate Reference System, CRS)는 crs() 함수를 통해 확인할 수 있습니다.

## crs() : coordinate reference system
crs(srtm_raster)
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs

dim() 함수는 레스터 데이터의 차원인 행의 개수(number of rows), 열의 개수 (number of columns), 층의 개수 (number of layers)를 한꺼번에 볼 수 있게 해줍니다. 

ncol() 은 열의 개수, nrow() 는 행의 개수, nlayers() 는 층의 개수를 개별적으로 접근할 수 있는 함수입니다.

## dim() : number of rows, columns and layers
dim(srtm_raster)
# [1] 457 465   1


nrow(srtm_raster)
# [1] 457

ncol(srtm_raster)
# [1] 465

nlayers(srtm_raster)
# [1] 1

ncell() 함수는 레스터 데이터의 전체 셀의 개수 (number of cells, or pixels) 를 알려주며, 이는 (셀의 개수 = 행의 개수 X 열의 개수) 를 의미합니다.

res() 함수는 레스터 데이터의 공간 해상도 (the raster's spatial resolution)를 나타냅니다.

## ncell() : number of cells (pixels) = number of rows * number of columns
ncell(srtm_raster)
# [1] 212505


## res() : the raster's spatial resolution
res(srtm_raster)
# [1] 0.0008333333 0.0008333333

extent() 함수는 공간의 x와 y 좌표의 최소, 최대값의 범위를 나타냅니다. 

반면에 summary() 함수는 각 셀의 특성 값의 분위수값과 평균의 요약 통계량을 나타냅니다.

## extent() : spatial extent
extent(srtm_raster)
# class      : Extent 
# xmin       : -113.2396 
# xmax       : -112.8521 
# ymin       : 37.13208 
# ymax       : 37.51292


## summary() : Summary of the values of a Raster* object (quartiles and mean)
summary(srtm_raster)
# srtm
# Min.    1024
# 1st Qu. 1535
# Median  1837
# 3rd Qu. 2115
# Max.    2892
# NA's       0
# Warning message:
# In .local(object, ...) :
# summary is an estimate based on a sample of 1e+05 cells (47.06% of all cells)

inMemory() 함수는 레스터 데이터가 메모리에 저장되어 있는지 (기본 설정), 아니면  디스크에 저장되어 있는지를 확인할 때 사용합니다. (블리언 반환)

## inmemory() : whether the raster data is stored in memory (the default) or on disk
inMemory(srtm_raster)
# [1] FALSE

레스터 데이터를 시각화하는데는 여러개의 패키지와 함수가 있는데요, 가장 기본적으로 쉽고 빠르게 시각화하는 방법은 plot() 함수를 사용하는 것입니다.

## plotting
plot(srtm_raster, main = 'basic raster plot')

[ Reference ]

* Raster data: https://geocompr.robinlovelace.net/spatial-class.html#raster-data

이번 포스팅이 많은 도움이 되었기를 바랍니다.

행복한 데이터 과학자 되세요!  :-)

지난번 포스팅에서는 지리공간 데이터의 포맷으로서 벡터 데이터(vector data)와 레스터 데이터(raster data) 모델에 대해서 간단하게 소개(https://rfriend.tistory.com/589)를 하였습니다.

이번 포스팅에서는 R의 sf 패키지를 활용해서 spData 에 내장되어 있는 세계지도(world map)를 시각화해보겠습니다.

(1) 지리공간 데이터 처리 및 분석을 위한 R 패키지 설치

(2) 요약통계량 summary() 및 부분 집합 subset 하기

(3) sf 패키지의 plot()으로 기본 지도 시각화하기

(4) 다른 지도 층을 추가하기 (add plots as layers)

(5) 원과 텍스트를 덮어쓰우기 (overlaying circles and texts over the plot)

R의 sf 패키지는 지리공간 벡터 데이터(vector data) 분석을 위한 패키지로서, 이전의 sp 패키지를 대체합니다. 그리고 GEOS와 GDAL과 R의 인터페이스를 제공하며 rgeos, rgdal 패키지를 대신합니다.

R의 raster 패키지는 지리공간 레스터 데이터(raster data)를 처리 및 분석하는데 사용합니다.

그리고 spData, spDataLarge 패키지에는 지리공간 데이터 샘플을 내장하고 있습니다. spDataLarge 패키지는 내장 데이터의 크기가 너무 크기 때문에 R CRAN 사이트에 패키지가 올라가 있지 않습니다. 따라서 CRAN에서 설치하려고 하면 에러가 나며, 아래처럼 Github 에서 type = "source" 로 해서 설치를 해주어야 합니다.

(만약 spData 패키지의 개발자 버전을 Github에서 바로 받아서 설치를 하려면 devtools::install_github("nowosad/spData") 로 하면 됩니다.)

참고로, spData 패키지에는 37개의 지리공간 데이터셋이 내장되어 있으며, https://nowosad.github.io/spData/  깃헙 페이지에 가면 전체 리스트를 확인할 수 있습니다. 아래의 데이터셋 리스트 중에서 'world' (World country polygons) 데이터셋을 가지고 sf 패키지를 사용해서 요약통계량 및 시각화를 예를 들어보겠습니다.

spData 패키지에 내장된 'world' 벡터 데이터에는 세계 대륙의 국가별로 10개의 집계데이터 칼럼과 마지막에 geom (지리공간 Polygon 리스트) 으로 구성되어 있습니다.

sf 패키지는 base R의 기본 함수를 그래도 사용할 수 있습니다. names()로 지리공간 데이터의 칼럼를 리스트업 할 수도 있고, summary() 함수로 요약통계량을 계산할 수 있습니다.

summary() 함수로 전세계 국가별 인구("pop")와 예상수명("lifeExp")에 대해 177개 Multipolygon 별로 최소/Q1/중앙값/Q3/최대값/결측값의 요약통계량을 계산해 보겠습니다.

world 데이터셋의 제일 마지막 'geom' 칼럼에는 177개 국가에 대한 WGS 84 좌표 시스템의MultiPolygon 정보 (좌표 리스트)가 들어있습니다.

행과 열의 부분집합을 subset 하는 방법도 기본 R 과 동일합니다. 아래 예에서는 1행~2행, 1열~3열까지의 데이터셋을 subset 해보겠습니다.

 world 세계지도 데이터 칼럼 중 7열~10열까지 4개 칼럼 별로 177 개 MultiPolygon별 값에 따라 색깔을 달리해서 R의 plot() 함수로 시각화를 해보겠습니다. ('area_km2', 'pop', 'lifeExp', 'gdpPercap' 칼럼의 값에 따라서 'geom' 의 177개 MultiPolygon 별 색깔이 달라짐.)

단 plot() 함수로 시각화를 하면 매우 간단하게 지도를 시각화할 수 있지만, 상호작용하는 동적인 지도(interactive map)가 아니라 그냥 정적인 지도만 그려지게 됩니다. (interactive map 은 다음번 포스팅에서 소개 예정)

위처럼 여러개의 칼럼에 대해 한꺼번에 여러개의 지도를 시각화하면 색깔이 무엇을 의미하는지 범례가 표시되지 않습니다.

반면에, 아래처럼 특정 한개의 변수에 대해서만 plot() 으로 시각화를 하면 색깔이 의미하는 바에 대한 범례(legend)를 색깔 막대 형태로 표시해줍니다.

plot() 으로 지리공간 데이터를 시각화할 때 col 매개변수로 특정 색깔을 지정해줄 수도 있습니다. 대신, 특정 색깔을 col 로 지정해주면 위의 지도 시각화처럼 오른쪽에 색깔 막대 범례는 필요가 없으므로, 비록 1개 변수만 시각화한다고 해도 범례는 없습니다.

위의 (3)번 처럼 지도를 먼저 한번 시각화를 하고, add = TRUE 매개변수를 사용하면 나중에 그 위에 다른 지도를 겹쳐서, 즉 층을 추가하여 지도를 덮어쓰기로 그릴 수 있습니다.

단, 첫번째 지도 그래프에 키(key)가 있을 경우에는 reset = FALSE 매개변수를 꼭 설정해준 다음에, 이후에 다음번 plot(add = TRUE)를 사용해주어야 합니다.

cex 매개변수를 사용하면 시각화한 지도 위에 원을 덮어쓰울 수 있습니다.

아래 예에서는 세계지도 위에 국가별 인구에 비례해서 원을 덮어쓰워보겠습니다. (pop 변수에 제곱근을 취하고 10000 으로 나누어준 것은 시각화한 지도 크기에 맞추어서 너무 크지도, 작지도 않은 원을 덮어쓰우기 위해 숫자를 조정한 것입니다.)

이번 포스팅이 많은 도움이 되었기를 바랍니다.

행복한 데이터 과학자 되세요!  :-)