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

지난번 포스팅에서는 좌표계 (CRS, Coordinate Reference System) 로서 지리 좌표계(Geographic CRS)와 투영 좌표계(Projected CRS) 그리고 R에서 좌표계를 조회하고 설정하는 방법에 대해서 알아보았습니다. (rfriend.tistory.com/618)

이번 포스팅에서는 좌표계 (CRS) 정보 안에 들어있는 공간의 단위 (Spatial Units)에 대해서 알아보겠습니다. 지도를 제작하거나 볼 때 측정 단위 (measurement units)가 미터(meters) 인지 혹은 피트(feets) 인지 명시적으로 표현하고 정확하게 확인할 필요가 있습니다. 그래야지 벡터의 지리적 데이터나 레스터의 픽셀에서 측정되는 단위라는 맥락(context)를 알 수 있고, 실제 지표면과 지도 표현 간의 관계, 거리를 알 수 있고, 또 거리나 면적 등을 계산할 수 있기 때문입니다.

(1) 지리공간 벡터 데이터에서의 측정 단위 (Units in Vector data)

(2) 지리공간 레스터 데이터에서의 측정 단위 (Units in Raster data)

(1) 지리공간 벡터 데이터에서의 측정 단위 (Units in Vector data)

sf 객체의 지리공간 벡터 데이터는 단위에 대해서 native support 이여서, 다른 외부 모듈이나 확장 프로그램을 설치하지 않아도 sf 객체 내에 단위가 들어가 있습니다. 그래서 sf 객체 벡터 데이터에 대해서 연산을 하게 되면 units 패키지에 의해 정의된 "단위 속성(units arrtibute)"도 같이 반환해주어서 단위로 인한 혼란을 미연에 방지해줍니다. (대부분의 좌표계는 미터(meters)를 사용하지만, 일부는 피트(feets)를 사용하기 때문에 단위가 혼란스러울 수 있음. raster 패키지는 단위가 native support가 아님.)

아래 예제는 R의 spData 패키지에 들어있는 "world" 데이터셋에서 '대한민국 (South Korea, Republic of Korea)' 의 벡터 데이터를 가져와서, sf 패키지의 st_area() 함수로 면적을 계산해본 것입니다.

## ========================
## Units in Geospatial data
## ========================

## -- getting 'world' data from spData package
library(spData)
names(world)
# [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"      "area_km2"  "pop"       "lifeExp"  
# [10] "gdpPercap" "geom"

south_korea = world[world$name_long == "Republic of Korea", ]

south_korea
# Simple feature collection with 1 feature and 10 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: 126.1174 ymin: 34.39005 xmax: 129.4683 ymax: 38.61224
# geographic CRS: WGS 84
# # A tibble: 1 x 11
# iso_a2 name_long  continent region_un subregion type  area_km2    pop lifeExp gdpPercap             geom
# <chr>  <chr>      <chr> <chr>  <chr>     <chr>   <dbl>  <dbl>   <dbl>   <dbl> <MULTIPOLYGON [arc_degre>
# 1 KR   Republic ~ Asia   Asia   Eastern ~ Sove~ 99044. 5.07e7    81.7  33426. (((126.1748 37.74969, 12~

plot(south_korea[1])

sf 패키지의 st_area() 함수로 벡터 데이터의 면적으로 계산하면, 결과값의 뒤에 [m^2] 이라고 해서 2차원 공간 상의 "제곱미터" 단위가 같이 반환된 것을 알 수 있습니다.

## -- calculating the area of South Korea
library(sf)
install.packages("lwgeom")
library(lwgeom)

st_area(south_korea)
# 99044366830 [m^2]

attributes() 함수로 벡터 데이터에 대한 면적 계산 결과에 대한 속성 정보를 확인해 볼 수 있습니다. "units" 클래스에 "m^2" (제곱미터) 단위라고 나오는 군요.

## -- attributes of units
attributes(st_area(south_korea))
# $units
# $numerator
# [1] "m" "m"
# 
# $denominator
# character(0)
# 
# attr(,"class")
# [1] "symbolic_units"
# 
# $class
# [1] "units"

위의 대한민국 면적 계산의 단위가 "제곱미터 [m^2]" 이다보니 결과값의 자리수가 너무 길게 나와서 보기에 어렵습니다.  그래서 보기에 편하도록 면적 계산의 단위를 "제곱킬로미터 [km^2]"로 변경하려면 units 패키지의 set_units(sf_object, units) 함수로 단위를 설정해주면 변경이 됩니다.

만약 set_units() 함수로 단위를 재설정하는게 아니라, 기존의 면적 단위인 '제곱미터(m^2)' 로 계산된 결과값을 1,000,000 으로 나누게 되면 결과값은 맞더라도 단위가 '제곱미터(m^2)' 로 그대로여서, 우리가 원하던 단위인 '제곱킬로미터(km^2)' 가 아니게 되므로 주의가 필요합니다. 즉, 단위를 변경해서 계산하고 결과값을 얻고 싶다면 units 패키지의 st_units() 함수로 단위를 명확하게 설정해주는게 필요합니다.

## -- cautious; the units is [m^2], not [km^2]
st_area(south_korea) / 1000000
# 99044.37 [m^2]


##-- set_units() : setting the units
library(units)
units::set_units(st_area(south_korea), km^2)
# 99044.37 [km^2]

(2) 지리공간 레스터 데이터에서의 측정 단위 (Units in Raster data)

앞에서 소개한 벡터 데이터를 다루는 sf 패키지는 단위가 native support 여서 조회나 계산 결과를 반환할 때 단위(units)를 속성으로 반환한다고 하였습니다. 하지만 레스터 데이터를 다루는 raster 패키지는 단위에 대해서 native support 가 아직까지는 아니므로, 단위에 대해서 혼란스러울 수 있으므로 조심해야 합니다.

아래의 예는 spDataLarge 패키지에 들어있는 strm.tif 파일을 raster() 함수로 읽어왔습니다. 이 데이터는 st_crs() 함수로 좌표계를 확인해보면 "WGS84 투영"을 사용하므로, 십진수 각도 (decimal degrees as units) 를 단위로 사용합니다. res() 함수로 해상도를 확인해보면, 단지 숫자형 벡터 (numeric vector) 만 반환할 뿐, 단위에 대한 속성 정보는 없습니다! (no units attributes).

## -- units in raster data
library(raster)
library(spDataLarge)

raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge") 
new_raster = raster(raster_filepath)


## -- getting CRS
st_crs(new_raster)
# Coordinate Reference System:
#   User input: +proj=longlat +datum=WGS84 +no_defs 
# wkt:
#   GEOGCRS["unknown",
#           DATUM["World Geodetic System 1984",
#                 ELLIPSOID["WGS 84",6378137,298.257223563,
#                           LENGTHUNIT["metre",1]],
#                 ID["EPSG",6326]],
#           PRIMEM["Greenwich",0,
#                  ANGLEUNIT["degree",0.0174532925199433],
#                  ID["EPSG",8901]],
#           CS[ellipsoidal,2],
#           AXIS["longitude",east,
#                ORDER[1],
#                ANGLEUNIT["degree",0.0174532925199433,
#                          ID["EPSG",9122]]],
#           AXIS["latitude",north,
#                ORDER[2],
#                ANGLEUNIT["degree",0.0174532925199433,
#                          ID["EPSG",9122]]]]


## -- WGS84 projection with decimal degrees as units
## -- cautious; simply returns numeric vector without any unit.
res(new_raster)
# [1] 0.0008333333 0.0008333333

만약 우리가 아래의 예처럼 UTM 투영을 사용한다면, 이에 따라서 단위가 바뀔 텐데요, res()로 해상도를 살펴보면 역시 단지 숫자형 벡터만 반환할 뿐, 단위 속성 정보는 없습니다.  따라서 R에서 raster 패키지로 레스터 데이터를 다룰 때는 단위에 대해서 특별히 조심할 필요가 있겠습니다.

## -- if we used the UTM projection, the units would change.
repr = projectRaster(new_raster, crs = "+init=epsg:26912")

## -- no units attributes, just only returns numeric vector
res(repr)
# [1] 73.8 92.5

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

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