레스터 객체에 대해 summary(raster_obj) 함수를 사용하면 레스터 객체의 픽셀 속성값에 대한 최소값(Minimum value), 1사분위수(1st Quantile), 중위값(Median), 3사분위수(3rd Quantile), 최대값(Maximum value), 결측값 개수(NA's) 에 대한 기술통계량(descriptive statistics)을 확인할 수 있습니다.
## summary() function for raster objects
summary(elev)
# layer
# Min. 1.00
# 1st Qu. 9.75
# Median 18.50
# 3rd Qu. 27.25
# Max. 36.00
# NA's 0.00
레스터 객체의 픽셀 속성값에 대해 cellStats(raster_obj, summary_function) 함수를 사용하여 평균(mean), 분산(variance), 표준편차(standard deviation), 사분위수(quantile), 합계(summation) 를 구할 수 있습니다.
RasterBlack Calss, RasterStack Class와 같이 여러개의 층을 가지는 레스터 클래스 객체 (multi-layered raster classes objects) 에 대해서 앞서 소개한 요약 통계량을 구하는 함수를 적용하면 각 층별로 따로 따로 요약 통계량이 계산됩니다.
(2) 레스터 객체 시각화 하기 (visualizing raster objects)
## -- visualization of raster objects
## plot()
plot(elev, main = "raster objects")
레스터 객체의 픽셀 속성값에 대해 히스토그램(histogram for raster objects)을 그리려면 hist() 함수를 사용하면 됩니다.
## histogram
hist(elev, main = "histogram for raster objects")
레스터 객체의 픽셀 속성값에 대해서 raster 패키지의 raster::density() 함수를 사용하여 추정 밀도 곡선(smoothed density estimates curve) 을 그릴 수 있습니다.
## raster::density() : density plot (smoothed density estimates)
density(elev, main = "density plot for raster objects")
레스터 객체의 픽셀 속성값에 대해 boxplot() 함수를 사용하면 상자그림(box plot for raster objects)을 그릴 수 있습니다.
## boxplot
boxplot(elev, main = "box plot for raster objects")
만약 레스터 객체에 대해 시각화를 하려고 하는데 안된다면 values() 또는 getValues() 함수를 사용해서 레스터 객체로 부터 픽셀의 속성값을 반환받은 결과에 대해서 시각화를 하면 됩니다.
지난번 포스팅에서는 3가지 종류의 레스터 클래스 객체로서 RasterLayer 클래스, RasterBrick 클래스, RasterStack 클래스가 있고, 각 특성은 무엇인지, 그리고 R raster 패키지를 사용해서 각 레스터 클래스 객체를 어떻게 만들 수 있는지 살펴보았습니다. (rfriend.tistory.com/617)
또 이전 포스팅 rfriend.tistory.com/596 에서 좌표계가 무엇이고, 벡터 데이터를 다룰 때 사용하는 sf 패키지를 사용해서 좌표계를 지정하고 확인하는 방법을 소개한 적이 있습니다.
이번 포스팅에서는 지리 공간 데이터 분석에서 가장 기본이 되고 또 처음에 확인을 해보아야 하는 좌표계, 좌표 참조 시스템(CRS, Coordinate Reference Systems)에 대해서 좀더 자세하게 소개하겠습니다. (이하, CRS는 "좌표계"로 부르겠습니다.)
(1) 좌표계, 좌표 참조 시스템 (CRS, Coordinate Reference Systems)은 무엇인가?
좌표계 (Coordinate Reference Systems) 는 지리공간 데이터를 저장/처리/분석하는데 있어서 지구의 표면을 어떻게 데이터의 공간 원소들로 표현할 지를 정의합니다. 좌표계는 좌표를 기반으로 국소, 지역, 또는 전지구적으로 지리상의 객체를 위치시킬 때 사용하는 참조 시스템입니다. 좌표계는 특정 지도 투영법(map projection)을 정의하고, 서로 다른 좌표계 간의 변환 (CRS transformtations)을 할 수 있도록 해줍니다.
좌표계는 지리 좌표계 (Geographic CRSs) 이거나 또는 투영 좌표계 (Projected CRSs) 인데요, 아래에 각각에 대해서 하나씩 설명하겠습니다.
지리 좌표계는 지구의 표면을 위도(Latitude)와 경도(Longitude) 의 두 개의 값을 사용해서 위치를 정의합니다. 위도(Latitude)는 적도면(Equator plane)으로 부터 남과 북 방향으로의 거리를 각도로 나타낸 것입니다. 경도(Longitude)는 본초 자오선 면(the Prime Meridian plane)으로 부터 동-서 방향으로의 거리를 각도로 나타낸 것입니다. 즉, 지리 좌표계에서의 거리는 미터(meters)가 아니라 각도(anglular distance)를 사용해서 측정합니다.
[ 위도와 경도 (Latitude and Longitude) ]
Latitude and Longitude (*Credit : Illinois State University), https://journeynorth.org/tm/LongitudeIntro.html
지리 좌표계에서의 지구 표면은 구면의 표면(Spherical surface) 또는 타원의 표면(Ellipsoidal surface)으로 표현됩니다.
구면 모델 (Spherical models) 은 지구가 주어진 반경(radius)으로 부터 완벽한 구면이라고 가정합니다. 구면 모델은 단순하고 계산이 빠르다는 장점이 있지만, 부정확하기 때문에 많이 이용되지는 않습니다. (실제의 지구 표면은 아래의 그림처럼 울퉁불퉁하며, 완벽한 구면이 아니기 때문입니다.)
타원 모델(Ellipsoidal models)은 적도 반지름(the equatorial radius)과 극 반지름(the polar radius)의 두 개의 모수를 사용해서 지구 표면을 정의합니다. 실제 지구는 아래의 그림처럼 위아래가 눌려서 옆으로 찌그러져 있으므로 타원 모델로 지구 표면을 표현하는 것이 구면 모델보다 더 적합합니다. (적도 반지름이 극 반지름보다 약 11.5 km 더 길이가 깁니다.)
[ 지구의 표면, 지오이드(GEOID), 타원(Ellipse), and 구(Sphere) ]
타원면(Ellipsoids)은 좌표계(CRSs)의 광의의 구성요소 중의 하나인 'the datum' 의 한 부분입니다. 'the datum'은 PROJ CRS 라이브러리 안의 'ellps' 모수로 무슨 타원면을 사용할지에 대한 정보와, 데카르트 좌표(Cartesian coordinates)와 지구 표면 위치 간의 정확한 관계에 대한 정보를 포함하고 있습니다. 이러한 추가적인 세부사항들은 'proj4string' 표기의 'towgs84' 매개변수에 저장이 되어 있습니다.
이러한 타원면의 the datum 정보들은 지구 표면의 지역적인 변이를 허용합니다. (가령, 큰 산맥이 있다면 지역 좌표계에 반영할 수 있음.)
the datum 에는 지역적 데이터(local datum)와 지심에서 바라 본 데이터(geocentric datum)의 두 종류가 있습니다. 'NAD83'과 같은 지역적 데이터(local datum)에서는 타원 표면이 특정 위치의 표면과 정렬이 되어 이동이 되어 있습니다. 반면에, 'WGS84'와 같은 지심에서 바라 본 데이터 (geocentric datum) 에서는 지구 중력의 중심이 중심 위치가 되며, 투사의 정확도가 특정 위치에 최적화되어 있지는 않습니다.
투영 좌표계(Projected Coordinate Reference Systems)는 암묵적으로 평평한 표면 위의 데카르트 좌표를 기반으로 합니다. 투영 좌표계는 원점(origin), x축과 y축, 그리고 미터(meters)와 같은 선형 측정 단위를 가지고 있습니다. 모든 투영 좌표계는 지리 좌표계(geographic CRSs)를 기반으로 하며, 3차원의 지구 표면을 투영 좌표계의 x축(동-서)와 y축(남-북)의 값으로 변환하는 지도 투사에 의존합니다.
평평한 평면을 가정하고 투영 변환을 하기 때문에 면적, 방향, 거리, 모양 등에 대한 왜곡이 발생하게 됩니다. 투영된 좌표계는 이들 중 단지 하나 혹은 두개의 특성 만을 보존할 수 있습니다. 투영 후 좌표계가 보존하는 특성이 무엇이냐에 따라서 보통 이름이 지어지는데요, 가령, 등면적 보존 면적 (equal-area preserve area), 방위각 보존 방향 (azimuthal preserve direction), 등거리 보존 거리 (equidistant preserve distance), 등각 보존 국소 형태 (conformal preserve local shape) 등과 같이 부릅니다.
투영의 유형과 원리(projection types and principles)에는 크게 평면(planar), 원뿔(conic), 실린더(cylindrical) 의 3가지 그룹이 있습니다. (아래 그림 참조)
(a) 평면 투영 (planar projection) 은 데이터를 지구와 접촉하는 하나의 점이나 접촉하는 선과 함께 평평한 표면에 투사를 합니다. 평면 투영법은 일반적으로 극 지방(polar regions)을 매핑할 때 사용됩니다. 투영의 기준(시작)이 되는 표준 점(standard point)으로 부터 멀어질 수록 데이터 왜곡이 심해집니다.
(b) 원뿔 투영 (conic projection)에서는 지구의 표면이 원뿔과 접촉하는 하나의 선이나 두개의 선들로 원뿔에 투영이 됩니다. 원뿔과 만나는 선인 표준 선(standard line)에서 왜곡이 가장 작고, 표준 선으로부터 멀어질 수록 왜곡이 심해집니다. 따라서 중간 위도 지역(mid-latitude areas)의 투영법으로 가장 적합합니다.
(c) 원통 투영 (cylindrical projection) 은 지구의 표면을 원통과 접촉하는 하나의 선이나 두개의 선들을 원통에 투사합니다. 원통 투영은 보통 지구 전체(entire world)를 매핑할 때 많이 사용합니다.
[ 지도 투영의 3가지 유형(Projection types)과 각 투영 유형별 왜곡(Distortion) ]
R에서 좌표계를 표현할 때는 (a) epsg 코드 (epsg code)나 또는 (b) proj4string 정의 (proj4string definition)를 사용합니다. 이들 두 가지의 좌표계 표현 방법별로 장단점이 있습니다.
(a) epsg 코드는 더 짧고, 따라서 기억하기 쉽습니다. 그리고 epsg 코드는 단 하나의 잘 정의된 좌표계를 지칭합니다. (EPSG: European Petroleum Survey Group, 지도 투영과 datums 에 대한 좌표계 정보 데이터베이스를 제공함)
반면에, (b) proj4string 정의는 투사 유형이나 datum, 타원체 등의 다른 모수들을 구체화할 수 있는 유연성이 있다는 장점이 있는 대신에, 사용자가 구체화를 해야 하므로 길고 복잡하며 기억하기 어렵습니다.
지리공간 R 패키지들은 PROJ 라이브러리를 사용하여 광범위한 좌표계(CRSs)를 지원합니다. R의 rgdal::make_EPSG() 함수를 사용하면 R이 지원하는 투영 좌표계의 리스트를 data.frame 형태로 반환받아 확인할 수 있습니다. 아래에 조회해본 결과, 현재 6,609개의 좌표계를 지원하고 있으며, EPSG 코드 (code), 설명 (note), proj4string 정의 (prj4), 투영 방법 (prj_method) 의 4개 칼럼으로 좌표계 정보를 제공합니다. 가령 'EPSG코드가 4326' 인 행을 조회해보면 note에 'WGS 84' (World Geodetic System 1984) 라고 설명이 나오는군요.
(4-2) 레스터 지리 데이터의 좌표계 확인: raster 패키지의 projection() 함수
(4-1) 벡터 데이터의 좌표계 (CRSs in Vector geographic data)
벡터 지리 데이터에 대해서는 sf 패키지의 st_crs() 함수를 사용해서 좌표계를 확인할 수 있습니다.
아래의 예에서는 spDataLarge 패키지에 들어있는 Zion 국립 공원의 경계를 다각형면(Polygon)으로 나타내는 zion.gpkg 벡터 데이터를 st_read() 함수로 불러와서, st_crs() 함수로 좌표계를 조회해 본 것입니다. 좌표계가 'User input: UTM Zone 12, Northern Hemisphere' 라고 나오는군요.
만약 좌표계가 비어있거나 잘못 입력되어 있는 경우 st_set_crs(vector_object, EPSG code) 구문으로 좌표계를 설정할 수 있습니다. 단, 아래의 경고 문구처럼 "st_set_crs() 함수는 좌표계를 변경하는 것이 투영 데이터를 변환하는 것은 아니며, 투영 데이터 변환을 하려면 st_transform() 함수를 이용 (st_crs<- : replacing crs does not reproject data; use st_transform for that)"해야 합니다.
## -- st_set_crs() : setting a CRS (coordinate reference system)
vector_zion_2 = st_set_crs(vector_zion, 4326) # set CRS with EPSG 4326 code
# Warning message:
# st_crs<- : replacing crs does not reproject data; use st_transform for that
[ 벡터 데이터: 지리 좌표계 (WGS 84; 왼쪽) vs. 투영 좌표계 (NAD83/ UTM zone 12N; 오른쪽) ]
레스터 모델의 객체에 대해서는 raster 패키지의 projection() 함수를 사용해서 좌표계를 확인하거나 설정할 수 있습니다.
아래 예는 R spDataLarge 패키지에 포함되어 있는 strm.tif 파일을 raster() 함수를 사용해 레스터 데이터로 읽어와서, projection() 함수로 좌표계를 확인해본 것입니다.
## -- raster::projection() : get or set CRS in raster* objects
library(raster)
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
new_raster = raster(raster_filepath)
projection(new_raster) # get CRS in raster objects
# [1] "+proj=longlat +datum=WGS84 +no_defs"
레스터 데이터에 대해서 좌표계를 새로 설정할 때도 역시 projection() 함수를 사용합니다. 벡터 데이터의 경우 좌표계를 설정할 때 'EPSG 코드'나 'Proj4string 정의' 모두 사용 가능한 반면에, 레스터 데이터는 'Proj4string 정의'만 사용할 수 있습니다.
## -- set CRS for raster objects
projection(new_raster) = "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs"
# Warning message:
# In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
# Discarded datum Unknown based on GRS80 ellipsoid in Proj4 definition
projection(new_raster)
# [1] "+proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs"
지리적 레스터 데이터 (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() 함수를 통해 확인할 수 있습니다.
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')
이번 포스팅에서는 앞서의 포스팅들에서 소개했던 내용을 모두 포괄하는 sf패키지의 "sf 클래스 (The sf class)" 를 만드는 방법을 소개하겠습니다.
(1) sf 패키지의 "sf 클래스(the sf class)"는 무엇인가?
(2) sf 클래스 만들기: st_sf()
(3) sf 클래스의 클래스, 기하유형, CRS 확인하기: class(), st_geometry_type(), st_crs)
(1) sf 패키지의 "sf 클래스(the sf class)"는 무엇인가?
R의 sf 패키지의 "sf 클래스(the sf class)"는 위의 그림에서 설명하는 바와 같이,
>> 지리공간 벡터 데이터 (GeoSpatial Vector Data)의 단순 지리특성 기하 (sfg, Simple feature geometry) 객체들을 리스트로 하나의 객체로 합친 단순 지리특성 칼럼 (sfc, Simple feature columns) 과,
>> 지리공간 벡터의 단순 지리특성 기하별로 속성 값 (예: 이름, 용도, 면적, 가격, 인구수, 온도, 습도 등)들을 모아놓은 data.frame 데이터 구조를
하나의 "sf 클래스" ( = sfc + data.frame) 로 합쳐 놓은 것입니다.
실제 지리공간 벡터 데이터를 분석한다고 했을 때는 sfc 객체의 기하 (Geometry) 정보와 data.frame 의 속성 (Non-geometric Attributes) 정보가 모두 있어야지, 기하 유형(점, 선, 면/다각형 등) 별로 속성 값을 분석할 수 있겠지요.
R의 sf 패키지로 지리공간 벡터 데이터 분석을 한다고 했을 때 새로운 용어가 자꾸 나오다 보니 위계 체계 (hierarchy), 자료 유형과 구조가 혼란스러울 수 있는데요, 위의 그림의 위계 체계, Input/Output 관계를 참고하시기 바랍니다.
(2) sf 클래스 만들기: st_sf()
(2-1) 기하 객체 (Geometry object)
아래의 예에서는 먼저, library(sf) 로 sf 패키지를 로딩하고, 지리공간 벡터 데이터(Vector data)의 점 기하 유형 (Point geometry types) 들로 이루어진 단순 지리특성 기하 (sfg, Simple feature geometry) 객체를 만들어 보겠습니다.
다음으로 st_sfc() 함수를 사용하여 앞서 만든 3개의 점(points) 단순 지리특성 기하(sfg) 객체를 하나의 단순 지리특성 칼럼 (sfc, Simple feature columns) 객체로 합치고, 이때 좌표 참조 시스템(CRS, Coordinate Reference Systems) 은 EPSG 코드 정의를 사용해서 WGS84 (World Geodetic System 1984) 인 crs = 4326 을 설정해주겠습니다.
이번에는 서울의 COEX, GBC(현대자동차 Global Business Center), JansilStadium(잠실 종합 운동장) 의 3개 POI (Point of Interest)별 면적(area)과 가격(price) 속성(Attributes) 정보를 가지는 data.frame을 만들어보겠습니다. (아래 속성 데이터 값은 모두 그냥 가짜로 입력한 예시 값입니다.)
## (3) data.frame object with non-geometric attributes seoul_df_attrib = data.frame( poi = c("COEX", "GBC", "JamsilStadium"), # Point of Interest area = c(20, 25, 15), price = c(80, 100, 50) )
seoul_df_attrib # poi area price # 1 COEX 20 80 # 2 GBC 25 100 # 3 JamsilStadium 15 50
(2-3) sf 클래스 = sfc 기하 객체(sfc Geometry) + data.frame 속성 (Non-geometry Attributes)
마지막으로, st_sf() 함수를 사용하여 위의 (2-1)에서 정의한 벡터 데이터의 기하 정보와, (2-2)에서 정의한 비기하 속성값 data.frame (Non-geometry Attributes) 을 합쳐서 sf 클래스 (the sf class) 를 만들어보겠습니다.
sf 클래스의 클래스를 class() 함수를 사용해서 확인해보면, 기하(Geometry) 정보를 가지고 있는 "sf" 클래스 객체와 속성 정보 값을 가지고 있는 "data.frame" 클래스의 두개 클래스 객체로 구성되어 있음을 알 수 있습니다.
## check the class class(seoul_sf) # [1] "sf" "data.frame"
st_geometry_type() 함수를 사용해서 기하 유형을 확인해보면, 이 예제에서는 점(Point) 3개로 구성되어 있음을 알 수 있습니다.
## check the geometry types st_geometry_type(seoul_sf) # [1] POINT POINT POINT # 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING MULTIPOLYGON ... TRIANGLE
좌표 참조 시스템(CRS, Coordinate Reference Systems)은 st_crs() 함수로 확인할 수 있는데요, 위의 (1)번에서 단순 지리특성 칼럼(sfc, Simple feature columns)을 생성할 때 crs = 4326 으로 설정을 해주었기 때문에 User input: EPSG:4326 으로 등록되어 있음을 확인할 수 있습니다.
지난번 포스팅에서는 R의 sf 패키지를 사용해서 지리공간 벡터 데이터의 기하 유형(geometry types of vector data)로서 단순 지리특성 기하(Simple feature geometry)인 점(Point), 선(LineString), 면(다각형) (Polygon), 다중점(MultiPoint), 다중선(MultiLineString), 다중면(MultiPolygon), 기하집합(GeometryCollection)에 대해서 소개하였습니다 .
이번 포스팅에서는 앞서 소개했던 벡터 데이터의 '단순 지리특성 기하 (Simple feature geometry, sfg)' 들을 sf 패키지의 st_sfc() 함수를 사용하여 '단순 지리특성 칼럼 (Simple feature column, sfc)' 으로 합치는 방법을 소개하겠습니다.
여기서 단순 지리특성 칼럼(Simple feature column, sfc)은 단순 지리특성 기하들의 리스트(a list of sfg) 로서, 두 개의 지리특성(features)를 하나의 칼럼 객체로 합쳐놓은 것입니다. (한국말로 feature, sfc, sfg 용어를 번역하기가 쉽지 않네요.)
-- 동일한 단순 지리특성 기하 유형 합치기
(1) 두개의 단순 지리특성 기하 점(2 sfg points)를 한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
(2) 두개의 단순 지리특성 기하 면(2 sfg polygons)를 한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
-- 서로 다른 단순 지리특성 기하 유형 합치기
(3) 단순 지리특성 기하 점과 면을 합쳐서 한개의 단순 지리특성 칼럼(1 sfc) 객체로 만들기
(1) 두개의 단순 지리특성 기하 점(2 sfg points)를 st_sfc() 함수로 한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
R의 sf 패키지를 사용할 것이므로, 먼저 library(sf) 로 패키지를 불러옵니다.
다음으로, 예제로 사용할 점(Point) 두개를 st_point() 함수를 사용해서 만들었습니다.
그리고, st_sfc() 함수를 사용하여, 단순 지리특성 기하(Simple feature geometry, sfg)인 두개의 점을
--> 단순 지리특성 칼럼(Simple feature columns, sfc) 의 하나의 객체로 합쳤습니다.
library(sf)
## Simple feature geometry (sfg) ## Point point_sfg_1 = st_point(c(1, 2)) point_sfg_1 # POINT (1 2)
point_sfg_2 = st_point(c(3, 4)) point_sfg_2 # POINT (3 4)
## Simple feature columns (sfc) ## st_sfc() function combine two simple features into one object with two features. ## sfc represents the geometry column in sf data frames. points_sfc = st_sfc(point_sfg_1, point_sfg_2)
points_sfc # Geometry set for 2 features # geometry type: POINT # dimension: XY # bbox: xmin: 1 ymin: 2 xmax: 3 ymax: 4 # CRS: NA # POINT (1 2) # POINT (3 4)
st_geometry_type() 함수를 사용하면 기하 유형을 확인해 볼 수 있습니다. 동일한 기하 유형인 두 개 점을 합쳤으므로 하나의 객체 집합에 두 개의 점(Points)가 들어있습니다.
## check the geometry type of sfc
st_geometry_type(points_sfc) # [1] POINT POINT # 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING MULTIPOLYGON ... TRIANGLE
(2) 두개의 단순 지리특성 기하 면(2 sfg polygons)를 st_sfc() 함수로
한개의 단순 지리특성 칼럼(1 sfc) 객체로 합치기
이번에는 동일한 기하 형태인 두개의 단순 지리특성 기하 면 (2 Simple feature geometry polygons) 을 st_polygon() 함수로 만들어보겠습니다.
그리고 이를 st_sfc() 함수를 사용해서 하나의 단순 지리특성 칼럼(Simple feature column) 객체로 합쳐보겠습니다.
## Polygon using list of matrices of points polygon_list_1 = list(rbind(c(1, 1), c(1, 3), c(3, 5), c(3, 1), c(1, 1))) polygon_sfg_1 = st_polygon(polygon_list_1) polygon_sfg_1 # POLYGON ((1 1, 1 3, 3 5, 3 1, 1 1))
이번 포스팅에서는 지리공간 데이터 중에서 실제 세계를 도형으로 표현(graphical representation) 하는 벡터 데이터 모델 (Vector data model)에서 사용하는 기하 유형 (geometry types)에 대해서 소개하겠습니다.
(* 참고로 레스터 데이터는 실제 세계를 픽셀로 이루어진 그리드 안에 지리특성 값을 가지는 형태라고 하였습니다.)
벡터 데이터의 기하 유형에는 총 17개가 있는데요, 이번 포스팅에서는 그중에서도 많이 사용되는 아래의 7개의 대표적인 유형에 대해서만 설명을 하겠습니다.
(1) 점 (Point)
(2) 선 (LineString)
(3) 면 (다각형) (Polygon)
(4) 다중점 (MultiPoint)
(5) 다중선 (MultiLineString)
(6) 다중면 (MultiPolygon)
(7) 기하집합 (GeometryCollection)
총 17개의 기하 유형 중에서 위에서 설명한 7개 이외에 나머지 유형에 대해서는 PostGIS 사이트를 참고하시기 바랍니다. ( http://postgis.net/docs/using_postgis_dbmanagement.html )
[ 벡터 데이터의 기하 유형 (Geometry Types of Vector data) ]
(1) 점 (Point), 다중점 (MultiPoint)
점은 2D (위도, 경도), 3D (위도, 경도, 해발고도), 또는 4D (위도, 경도, 고도, 측정 정확도) 공간 상의 하나의 위치 좌표에 해당합니다.
예: Point (3 3)
점 (Point)은 0 차원으로서, 각 기하도형의 가장 기본이 되는 형태입니다. 여러개의 점을 이어주면 선이 되고, 선을 닫아주면 면이 되는 식이니깐요. 이렇게 점 --> 선 --> 면 혹은 점 --> 다중점, 선 --> 다중선, 면 --> 다중면 방향으로 기하도형을 변환하는 것을 구성(composition) 한다고 하며, 그 반대 방향으로의 변환을 해체(Decompostion) 한다고 합니다.
점은 우물, 랜드마크, 관심지역(Point of Interest) 등의 지리특성을 정의할 때 사용합니다.
R에서 점(point)을 만들 때 st_point() 함수를 사용합니다. st_point(dim="XYM") 에서 dim 은 차원(dimension)을 의미하고, M은 보통 Measurement 의 M (측정 정확도, measurement accuracy)을 의미합니다. 점의 2D, 3D, 4D 좌표는 c() 로 묶어줍니다.(concatenate)
참고로, 접두사 'ST_' 는 PostGIS의 함수에서도 동일하게 사용되며, 'standard Spatial Type (ST)' 의 앞글자를 딴 것입니다.
install.packages("sf")
library(sf)
##---------------- ## point
## XY point st_point(c(3, 3)) # POINT (3 3)
## XYZ point st_point(c(3, 3, 2)) # POINT Z (3 3 2)
## XYM point st_point(c(3, 3, 1), dim = "XYM") # POINT M (3 3 1)
## XYZM point st_point(c(3, 3, 2, 1)) # POINT ZM (3 3 2 1)
(2) 선 (LineString)
선(LineString) 은 점들을 직선으로 연결한, 여러개 점들의 연속적인 집합입니다. 선은 아래의 예처럼 몇 개 점들의 집합으로 표현할 수 있습니다(이들 점을 순서대로 직선으로 연결).
예: LineString (1 1, 1 3, 3 5)
선은 1차원(dimension of 1)이며 길이(length)가 있습니다. (반면, 점은 0 차원이고 길이가 없음)
선은 길, 철도, 강, 전선 등과 같이 선형 특성을 가진 지형지물을 정의하는데 사용합니다.
R에서 선(LineString)은 st_linestring() 함수를 사용하여 만들 수 있습니다. 선을 이루는 여러개 점들의 좌표는 rbind()로 행렬을 만들어서 사용합니다.
면(Polygon)은 2차원의 표면을 가지는 형태로 점들이 연결된 것으로서, 안에 구멍(hole)이 있을 수도 있고 또는 없을 수도 있습니다. (위의 그림 예에서는 면 안에 구멍이 없음). 면이 되기 위해서는 시작점과 끝점이 동일해야 하며, 이를 닫혀있어야 한다(closed, not open)고 말합니다. 아래 예의 경우 면의 시작점과 끝점이 (x=1, y=1) 로서 동일하며, 즉 닫혀있습니다.
예: Polygon ((1 1, 1 3, 3 5, 3 1, 1 1))
면(Polygon)은 육지, 저수지, 건물 등과 같이 공간적인 외곽 (spatial extent)과 면적(area)을 가진 지리지형을 나타내는데 사용할 수 있습니다.
R의 st_polygon() 함수를 사용해서 면(Polygon)을 만들 수 있습니다. 면을 만들 때는 점 행렬의 리스트(list of matrices of points)를 입력값으로 사용합니다.
## Polygon using list of matrices of points polygon_list = list(rbind(c(1, 1), c(1, 3), c(3, 5), c(3, 1), c(1, 1)))
다중점(MultiPoint)은 점을 여러개 모아놓은 것(a collection of points)입니다.
예: MultiPoint (1 1, 3 3, 3 5)
R에서 점은 st_point() 함수로 만들었다면, 다중점(MultiPoint)은 st_multipoint() 로 만들 수 있습니다. 이때 st_point() 에서는 점의 좌표를 c()로 묶어주기만 했다면, 다중점(MultiPoint)을 만들 때는 rbind()로 묶어서 점의 행렬(matrices)을 만들어서 st_multipoint() 함수를 사용합니다.
## MultiPoint using matrices of points multipoint_matrix = rbind(c(1, 1), c(3, 3), c(3, 5))
R에서 면(Polygon)은 st_polygon() 함수로 만들었다면, 다중면(MultiPolygon)은 st_multipolygon() 함수로 만들 수 있습니다. 면(Polygon)을 만들 때 점 행렬의 리스트를 st_polygon() 함수에 넣어주었다면, 다중면(MultiPolygon)을 만들 때는 여러개의 면이 모인 것이므로 점 행렬 리스트들의 리스트 (list of lists)를 st_multipolygon() 함수에 넣어줍니다.
기하(도형)집합(GeometryCollection)은 위에서 설명한 점(Point), 선(LineString), 면(다각형)(Polygon), 다중점(MultiPoint), 다중선(MultiLineString), 다중면(MultiPolygon) 중에서 여러개의 기하 유형을 모아놓은 것이며, 어떤 조합(any combination of geometries)도 가능합니다.
아래 예에서는 점(Point), 다중선(MultiLineString), 면(Polygon)이 혼합된 조합의 기하집합(GeometryCollection) 입니다.
R의 st_geometrycollection() 함수를 사용해서 기하집합(GeometryCollection)을 만들 수 있습니다. 이때 각 기하 유형을 리스트(list)로 묶어서 넣어주면 됩니다.
## GeometryCollection using list of combination of geometries point_coord = c(3, 3) multilinestring_list = list(rbind(c(1, 1), c(1, 3), c(3, 5)), rbind(c(2, 2), c(5, 2))) polygon_list = list(rbind(c(4, 3), c(4, 5), c(5, 5), c(5, 3), c(4, 3)))
지리공간 데이터 (GeoSpatial data)를 처리하고 분석하는데 있어서 첫번째 관문이자 큰 도전사항 중에 하나가 지리공간 데이터 포맷이 매우 다양하다는 것입니다.
아래에 다양한 지리공간 데이터(various GeoSpatial data foramts)의 리스트를 소개하고, 특히 이중에서 점, 선, 다각형으로 구성된 벡터 데이터 포맷의 이미지 시각화를 예시로 보였습니다. 지리공간 데이터 포맷이 상당히 많지요?
이들 지리공간 데이터 포맷별로 데이터를 DB나 R로 불러오기 (importing)할 때 사용하는 DB utility tools 이나 R의 package가 달라지다 보니 번거롭고 또 어려운 점이 있습니다.
[ 다양한 지리공간 데이터 포맷 (various GeoSpatial data formats) ]
R을 활용한 지리공간 데이터의 처리 및 분석, 시각화를 본격적으로 들어가기 전에 먼저, 이들 지리공간 데이터 포맷들 중에서 특히 벡터 데이터(Vector data)와 레스트 데이터 (Raster data) 모델에 대해서 이들이 무엇이고, 어떻게 활용이 되며, 무슨 R 패키지를 사용해서 분석할 수 있는지에 대해서 알아보겠습니다.
[ 지리공간벡터 데이터(Vector data) vs. 레스터 데이터 (Raster data) ]
(1) 지리공간 벡터 데이터 (Vector data)
벡터 데이터에는 KML(.kml or .kmz), GML, GeoJSON, Shapefile (.shp), WKT 등의 데이터 포맷이 있습니다.
GeoJSON 데이터 포맷은 이름에서 짐작할 수 있듯이 JSON 기반으로 간단한 지리공간 데이터와 그 외 일반 데이터를 저장합니다. GeoJSON 데이터는 인터넷으로 지리공간 & 일반 데이터를 공유하는데 역시 많이 사용됩니다.
Shapefile 데이터 포맷은 GIS (Geographic Information System) 소프트웨어를 위한 지리공간 벡터 데이터입니다. Shapefile 은 GIS 의 국제적인 제공사인 Esri(Environmental Systems Research Institute)에서 개발하고 관리하며, GIS 소프트웨어 간 상호운용성(interoperability)를 보장합니다.
WKT 데이터 포맷은 Well-Known Text 의 약자로서, 벡터 지리공간 데이터를 표현하는데 텍스트 마크업 언어(Text Markup Language)를 사용합니다. WKB (Well-Known Bianry)는 WKT와 같은 정보를 저장하는데 있어 이진(binary) 포맷을 사용해 보다 간소하고 컴퓨터가 처리하기에 편리하도록 하며, 대신 사람이 읽을 수는 없습니다.
벡터 데이터는 실제 세상을 그래픽으로 재표현(graphical representation of the real world)한 것으로서, 점, 선, 다각형(points, lines, polygons) 유형의 그래픽을 이용합니다. 벡터 데이터는 지구 표면의 객체나 특징을 일반화하여 표현하는데 사용됩니다.
벡터 데이터는 별개로 분리되고, 경계가 잘 정의되어 있어서 보통 높은 수준의 정밀도 (high level of precision) 을 가지고 있습니다. 이런 이유로 벡터 데이터는 사회 과학 (social sciences) 분야에서 많이 사용됩니다.
R 의 sf 패키지 (spatial data frame) 를 사용하여 벡터 데이터를 불러오고, 처리 및 분석, 시각화를 할 수 있습니다. (다음 포스팅에서 소개) sf 패키지는 이전의 sp 패키지, rgeos, rgdal 패키지를 모두 아우르고 있고, GEOS, GDAL, PROJ 와 R 의 interface를 제공해주어서, R로 지리공간 벡터 데이터를 다루는데 있어 매우 편리하고 강력합니다.
[ 강과 도심 지역을 나타낸 벡터 데이터(vector data)와 레스터 데이터(raster data) 비교 ]
레스터 데이터(Raster data)에는 ESRI Grid, GeoTIFF, JPEG 2000, NITF 등이 데이터 포맷이 있습니다.
레스터 데이터는 픽셀의 격자(grid of pixels) 로 지구의 표면을 표현합니다. 각 픽셀 안에는 색, 측정 단위 등과 같이 질문의 요소에 대한 정보를 전달하는 값이 있습니다.
레스터 데이터는 인공위성이나 항공장비에서 지구 표면을 향해 위에서 아래로 수직으로 찍은 사진으로 생각하면 이해하기가 쉽습니다.(예: NASA에서 제공하는, 인공위성에서 찍은 지구의 야간 사진 등) 이 지구표면을 수직으로 찍은 사진을 픽셀의 격자로 나누어서 각 픽셀(pixel, cell)에 지리특성정보 값을 가지고 있는 것입니다.
레스터 데이터의 품질은 사진을 찍었던 장비의 해상도의 한계나, 활용하고자 하는 분야의 목적에 따라서 다양합니다. 레스터 데이터는 많은 환경관련 과학 분야 (environmental sciences)에 많이 사용되고 있습니다.
R의 raster 패키지를 사용하면 R에서 레스터 데이터를 처리할 수 있습니다.
위에서 각각 소개한 벡터 데이터와 레스터 데이터는 상호 간에 변환(converting from vector to raster data, from raster to vector data)이 가능하며, 하나의 분석 목적에 두 유형의 데이터 포맷이 동시에 사용되기도 합니다.
다음번 포스팅에서는 R의 spData 패키지에 내장되어 있는 지리공간 벡터 데이터 모델(Vector data model)을 가지고 sf 패키지로 시각화하는 간단한 예를 소개하겠습니다.
