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

지난번 포스팅에서는 레스터 데이터의 유형으로 (a) 셀 ID (Cell IDs), (b) 셀 값 (Cell Values) 에 대해서 알아보고, RasterLayer 예제 데이터를 R raster 패키지의 raster() 함수로 불러와서 여러 속성 정보를 살펴보았습니다. (rfriend.tistory.com/605)

이번 포스팅에서는 3가지의 레스터 클래스 (Raster Classes) 의 특장점과 언제 사용하면 좋을지에 대해서 소개하겠습니다.

(1) RasterLayer class
(2) RasterBrick class
(3) RasterStack class

(4) 언제 어떤 레스터 클래스를 사용하는 것이 좋은가?

[ 3가지의 레스터 클래스 (Raster Classes) ]

(1) RasterLayer Class

RasterLayer Class는 레스터 객체 중에서 가장 간단한 형태의 클래스로서, 단 한 개의 층으로 구성되어 있습니다. 

지난번 포스팅에서 소개했던 것처럼, RasterLayer Class 객체를 만드는 가장 쉬운 방법은 기존의 RasterLayer Class 객체 파일을 읽어오는 것입니다.  아래 예에서는 raster 패키지의 raster() 함수를 사용해서 spDataLarge 패키지에 내장되어 있는 srtm.tif 레스터 층 클래스 객체를 읽어와서 raster_layer 라는 이름의 단 한개의 층만을 가진 RasterLayer Class 객체를 만들어보겠습내다 .  nlayers() 함수로 층의 개수를 살펴보면 '1'개 인 것을 확인할 수 있습니다.

## ==========================================
## R Raster Classes
## : RasterLayer, RasterBrick and RasterStack
## ==========================================

library(raster)
library(spDataLarge)

## -- (1) RasterLayer
## : The RasterLayer class represents the simplest form of a raster object, 
##   and consists of only one layer.

## read-in a raster file from disk or from a server
raster_filepath = system.file("raster/srtm.tif", package = "spDataLarge")
raster_filepath
# [1] "/Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/srtm.tif"

raster_layer = raster(raster_filepath)


raster_layer
# 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)


## number of layers
nlayers(raster_layer)
# [1] 1

raster 패키지는 rgdal 패키지의 도움을 받아 수많은 드라이버를 지원합니다. 

아래는 raster::writeFormats() 로 확인해 본, raster 패키지가 지원하는 "쓰기 포맷 (Write Formats)" 들 입니다.

## write Formats
raster::writeFormats()
#      name        long_name                                
# [1,] "raster"    "R-raster"                               
# [2,] "SAGA"      "SAGA GIS"                               
# [3,] "IDRISI"    "IDRISI"                                 
# [4,] "IDRISIold" "IDRISI (img/doc)"                       
# [5,] "BIL"       "Band by Line"                           
# [6,] "BSQ"       "Band Sequential"                        
# [7,] "BIP"       "Band by Pixel"                          
# [8,] "ascii"     "Arc ASCII"                              
# [9,] "CDF"       "NetCDF"                                 
# [10,] "ADRG"      "ARC Digitized Raster Graphics"          
# [11,] "BMP"       "MS Windows Device Independent Bitmap"   
# [12,] "BT"        "VTP .bt (Binary Terrain) 1.3 Format"    
# [13,] "BYN"       "Natural Resources Canada's Geoid"       
# [14,] "CTable2"   "CTable2 Datum Grid Shift"               
# [15,] "EHdr"      "ESRI .hdr Labelled"                     
# [16,] "ELAS"      "ELAS"                                   
# [17,] "ENVI"      "ENVI .hdr Labelled"                     
# [18,] "ERS"       "ERMapper .ers Labelled"                 
# [19,] "GPKG"      "GeoPackage"                             
# [20,] "GS7BG"     "Golden Software 7 Binary Grid (.grd)"   
# [21,] "GSBG"      "Golden Software Binary Grid (.grd)"     
# [22,] "GTX"       "NOAA Vertical Datum .GTX"               
# [23,] "GTiff"     "GeoTIFF"                                
# [24,] "HDF4Image" "HDF4 Dataset"                           
# [25,] "HFA"       "Erdas Imagine Images (.img)"            
# [26,] "IDA"       "Image Data and Analysis"                
# [27,] "ILWIS"     "ILWIS Raster Map"                       
# [28,] "INGR"      "Intergraph Raster"                      
# [29,] "ISCE"      "ISCE raster"                            
# [30,] "ISIS2"     "USGS Astrogeology ISIS cube (Version 2)"
# [31,] "ISIS3"     "USGS Astrogeology ISIS cube (Version 3)"
# [32,] "KRO"       "KOLOR Raw"                              
# [33,] "LAN"       "Erdas .LAN/.GIS"                        
# [34,] "Leveller"  "Leveller heightfield"                   
# [35,] "MBTiles"   "MBTiles"                                
# [36,] "MRF"       "Meta Raster Format"                     
# [37,] "NGW"       "NextGIS Web"                            
# [38,] "NITF"      "National Imagery Transmission Format"   
# [39,] "NTv2"      "NTv2 Datum Grid Shift"                  
# [40,] "NWT_GRD"   "Northwood Numeric Grid Format .grd/.tab"
# [41,] "PAux"      "PCI .aux Labelled"                      
# [42,] "PCIDSK"    "PCIDSK Database File"                   
# [43,] "PCRaster"  "PCRaster Raster File"                   
# [44,] "PDF"       "Geospatial PDF"                         
# [45,] "PDS4"      "NASA Planetary Data System 4"           
# [46,] "PNM"       "Portable Pixmap Format (netpbm)"        
# [47,] "RMF"       "Raster Matrix Format"                   
# [48,] "ROI_PAC"   "ROI_PAC raster"                         
# [49,] "RRASTER"   "R Raster"                               
# [50,] "RST"       "Idrisi Raster A.1"                      
# [51,] "SAGA"      "SAGA GIS Binary Grid (.sdat, .sg-grd-z)"
# [52,] "SGI"       "SGI Image File Format 1.0"              
# [53,] "Terragen"  "Terragen heightfield"                   
# [54,] "VICAR"     "MIPL VICAR file"                        
# [55,] "netCDF"    "Network Common Data Format" 

그리고 아래는 rgdal::gdalDrivers() 로 확인해 본, rgdal 패키지에서 지원하는 140 여개의 드라이버들입니다. 드라이버별로 생성(create), 복사(copy) 가능 여부는 아래에서 확인이 가능합니다.

## gdal Drivers
rgdal::gdalDrivers()
#                    name                                                long_name create  copy isRaster
# 1               AAIGrid                                      Arc/Info ASCII Grid  FALSE  TRUE     TRUE
# 2                  ACE2                                                     ACE2  FALSE FALSE     TRUE
# 3                  ADRG                            ARC Digitized Raster Graphics   TRUE FALSE     TRUE
# 4                   AIG                                     Arc/Info Binary Grid  FALSE FALSE     TRUE
# 5                   ARG                                Azavea Raster Grid format  FALSE  TRUE     TRUE
# 6                AirSAR                                AirSAR Polarimetric Image  FALSE FALSE     TRUE
# 7                   BAG                               Bathymetry Attributed Grid  FALSE  TRUE     TRUE
# 8                BIGGIF                       Graphics Interchange Format (.gif)  FALSE FALSE     TRUE
# 9                   BLX                                     Magellan topo (.blx)  FALSE  TRUE     TRUE
# 10                  BMP                     MS Windows Device Independent Bitmap   TRUE FALSE     TRUE
# 11                  BSB                              Maptech BSB Nautical Charts  FALSE FALSE     TRUE
# 12                   BT                      VTP .bt (Binary Terrain) 1.3 Format   TRUE FALSE     TRUE
# 13                  BYN                         Natural Resources Canada's Geoid   TRUE FALSE     TRUE
# 14                  CAD                                           AutoCAD Driver  FALSE FALSE     TRUE
# 15                 CALS                                            CALS (Type 1)  FALSE  TRUE     TRUE
# 16                 CEOS                                               CEOS Image  FALSE FALSE     TRUE
# 17                COASP                          DRDC COASP SAR Processor Raster  FALSE FALSE     TRUE
# 18                  COG                        Cloud optimized GeoTIFF generator  FALSE  TRUE     TRUE
# 19                COSAR               COSAR Annotated Binary Matrix (TerraSAR-X)  FALSE FALSE     TRUE
# 20                  CPG                                          Convair PolGASP  FALSE FALSE     TRUE
# 21                  CTG                           USGS LULC Composite Theme Grid  FALSE FALSE     TRUE
# 22              CTable2                                 CTable2 Datum Grid Shift   TRUE FALSE     TRUE
# 23                 DAAS          Airbus DS Intelligence Data As A Service driver  FALSE FALSE     TRUE
# 24              DERIVED               Derived datasets using VRT pixel functions  FALSE FALSE     TRUE
# 25                DIMAP                                               SPOT DIMAP  FALSE FALSE     TRUE
# 26                DIPEx                                                    DIPEx  FALSE FALSE     TRUE
# 27                 DOQ1                                     USGS DOQ (Old Style)  FALSE FALSE     TRUE
# 28                 DOQ2                                     USGS DOQ (New Style)  FALSE FALSE     TRUE
# 29                 DTED                                    DTED Elevation Raster  FALSE  TRUE     TRUE
# 30              E00GRID                                 Arc/Info Export E00 GRID  FALSE FALSE     TRUE
# 31              ECRGTOC                                          ECRG TOC format  FALSE FALSE     TRUE
# 32                EEDAI                              Earth Engine Data API Image  FALSE FALSE     TRUE
# 33                 EHdr                                       ESRI .hdr Labelled   TRUE  TRUE     TRUE
# 34                  EIR                                        Erdas Imagine Raw  FALSE FALSE     TRUE
# 35                 ELAS                                                     ELAS   TRUE FALSE     TRUE
# 36                 ENVI                                       ENVI .hdr Labelled   TRUE FALSE     TRUE
# 37                  ERS                                   ERMapper .ers Labelled   TRUE FALSE     TRUE
# 38                 ESAT                                     Envisat Image Format  FALSE FALSE     TRUE
# 39                 FAST                                        EOSAT FAST Format  FALSE FALSE     TRUE
# 40                  FIT                                                FIT Image  FALSE  TRUE     TRUE
# 41              FujiBAS                                   Fuji BAS Scanner Image  FALSE FALSE     TRUE
# 42                  GFF Ground-based SAR Applications Testbed File Format (.gff)  FALSE FALSE     TRUE
# 43                  GIF                       Graphics Interchange Format (.gif)  FALSE  TRUE     TRUE
# 44                  GMT                                   GMT NetCDF Grid Format  FALSE  TRUE     TRUE
# 45                 GPKG                                               GeoPackage   TRUE  TRUE     TRUE
# 46       GRASSASCIIGrid                                         GRASS ASCII Grid  FALSE FALSE     TRUE
# 47                 GRIB                             GRIdded Binary (.grb, .grb2)  FALSE  TRUE     TRUE
# 48                GS7BG                     Golden Software 7 Binary Grid (.grd)   TRUE  TRUE     TRUE
# 49                 GSAG                        Golden Software ASCII Grid (.grd)  FALSE  TRUE     TRUE
# 50                 GSBG                       Golden Software Binary Grid (.grd)   TRUE  TRUE     TRUE
# 51                  GSC                                              GSC Geogrid  FALSE FALSE     TRUE
# 52                  GTX                                 NOAA Vertical Datum .GTX   TRUE FALSE     TRUE
# 53                GTiff                                                  GeoTIFF   TRUE  TRUE     TRUE
# 54                  GXF                             GeoSoft Grid Exchange Format  FALSE FALSE     TRUE
# 55               GenBin                           Generic Binary (.hdr Labelled)  FALSE FALSE     TRUE
# 56                 HDF4                       Hierarchical Data Format Release 4  FALSE FALSE     TRUE
# 57            HDF4Image                                             HDF4 Dataset   TRUE FALSE     TRUE
# 58                 HDF5                       Hierarchical Data Format Release 5  FALSE FALSE     TRUE
# 59            HDF5Image                                             HDF5 Dataset  FALSE FALSE     TRUE
# 60                  HF2                               HF2/HFZ heightfield raster  FALSE  TRUE     TRUE
# 61                  HFA                              Erdas Imagine Images (.img)   TRUE  TRUE     TRUE
# 62                 HTTP                                    HTTP Fetching Wrapper  FALSE FALSE     TRUE
# 63                  IDA                                  Image Data and Analysis   TRUE FALSE     TRUE
# 64  IGNFHeightASCIIGrid                  IGN France height correction ASCII Grid  FALSE FALSE     TRUE
# 65                ILWIS                                         ILWIS Raster Map   TRUE  TRUE     TRUE
# 66                 INGR                                        Intergraph Raster   TRUE  TRUE     TRUE
# 67                 IRIS                             IRIS data (.PPI, .CAPPi etc)  FALSE FALSE     TRUE
# 68                 ISCE                                              ISCE raster   TRUE FALSE     TRUE
# 69                  ISG                      International Service for the Geoid  FALSE FALSE     TRUE
# 70                ISIS2                  USGS Astrogeology ISIS cube (Version 2)   TRUE FALSE     TRUE
# 71                ISIS3                  USGS Astrogeology ISIS cube (Version 3)   TRUE  TRUE     TRUE
# 72           JAXAPALSAR               JAXA PALSAR Product Reader (Level 1.1/1.5)  FALSE FALSE     TRUE
# 73                 JDEM                                      Japanese DEM (.mem)  FALSE FALSE     TRUE
# 74          JP2OpenJPEG               JPEG-2000 driver based on OpenJPEG library  FALSE  TRUE     TRUE
# 75                 JPEG                                                JPEG JFIF  FALSE  TRUE     TRUE
# 76      KMLSUPEROVERLAY                                        Kml Super Overlay  FALSE  TRUE     TRUE
# 77                  KRO                                                KOLOR Raw   TRUE FALSE     TRUE
# 78                  L1B                     NOAA Polar Orbiter Level 1b Data Set  FALSE FALSE     TRUE
# 79                  LAN                                          Erdas .LAN/.GIS   TRUE FALSE     TRUE
# 80                  LCP                        FARSITE v.4 Landscape File (.lcp)  FALSE  TRUE     TRUE
# 81               LOSLAS                        NADCON .los/.las Datum Grid Shift  FALSE FALSE     TRUE
# 82             Leveller                                     Leveller heightfield   TRUE FALSE     TRUE
# 83                  MAP                                         OziExplorer .MAP  FALSE FALSE     TRUE
# 84              MBTiles                                                  MBTiles   TRUE  TRUE     TRUE
# 85                  MEM                                         In Memory Raster   TRUE FALSE     TRUE
# 86                  MFF                                        Vexcel MFF Raster   TRUE  TRUE     TRUE
# 87                 MFF2                                 Vexcel MFF2 (HKV) Raster   TRUE  TRUE     TRUE
# 88                  MRF                                       Meta Raster Format   TRUE  TRUE     TRUE
# 89                 MSGN                           EUMETSAT Archive native (.nat)  FALSE FALSE     TRUE
# 90                  NDF                                        NLAPS Data Format  FALSE FALSE     TRUE
# 91             NGSGEOID                              NOAA NGS Geoid Height Grids  FALSE FALSE     TRUE
# 92                  NGW                                              NextGIS Web   TRUE  TRUE     TRUE
# 93                 NITF                     National Imagery Transmission Format   TRUE  TRUE     TRUE
# 94                 NTv1                                    NTv1 Datum Grid Shift  FALSE FALSE     TRUE
# 95                 NTv2                                    NTv2 Datum Grid Shift   TRUE FALSE     TRUE
# 96              NWT_GRC               Northwood Classified Grid Format .grc/.tab  FALSE FALSE     TRUE
# 97              NWT_GRD                  Northwood Numeric Grid Format .grd/.tab   TRUE  TRUE     TRUE
# 98                  OZI                                   OziExplorer Image File  FALSE FALSE     TRUE
# 99                 PAux                                        PCI .aux Labelled   TRUE FALSE     TRUE
# 100              PCIDSK                                     PCIDSK Database File   TRUE FALSE     TRUE
# 101            PCRaster                                     PCRaster Raster File   TRUE  TRUE     TRUE
# 102                 PDF                                           Geospatial PDF   TRUE  TRUE     TRUE
# 103                 PDS                               NASA Planetary Data System  FALSE FALSE     TRUE
# 104                PDS4                             NASA Planetary Data System 4   TRUE  TRUE     TRUE
# 105            PLMOSAIC                                  Planet Labs Mosaics API  FALSE FALSE     TRUE
# 106            PLSCENES                                   Planet Labs Scenes API  FALSE FALSE     TRUE
# 107                 PNG                                Portable Network Graphics  FALSE  TRUE     TRUE
# 108                 PNM                          Portable Pixmap Format (netpbm)   TRUE FALSE     TRUE
# 109                 PRF                                      Racurs PHOTOMOD PRF  FALSE FALSE     TRUE
# 110       PostGISRaster                                    PostGIS Raster driver  FALSE  TRUE     TRUE
# 111                   R                                      R Object Data Store  FALSE  TRUE     TRUE
# 112                 RDA                   DigitalGlobe Raster Data Access driver  FALSE FALSE     TRUE
# 113                 RIK                                  Swedish Grid RIK (.rik)  FALSE FALSE     TRUE
# 114                 RMF                                     Raster Matrix Format   TRUE FALSE     TRUE
# 115             ROI_PAC                                           ROI_PAC raster   TRUE FALSE     TRUE
# 116              RPFTOC                         Raster Product Format TOC format  FALSE FALSE     TRUE
# 117             RRASTER                                                 R Raster   TRUE  TRUE     TRUE
# 118                 RS2                                   RadarSat 2 XML Product  FALSE FALSE     TRUE
# 119                 RST                                        Idrisi Raster A.1   TRUE  TRUE     TRUE
# 120          Rasterlite                                               Rasterlite  FALSE  TRUE     TRUE
# 121                SAFE                              Sentinel-1 SAR SAFE Product  FALSE FALSE     TRUE
# 122                SAGA                  SAGA GIS Binary Grid (.sdat, .sg-grd-z)   TRUE  TRUE     TRUE
# 123            SAR_CEOS                                           CEOS SAR Image  FALSE FALSE     TRUE
# 124                SDTS                                              SDTS Raster  FALSE FALSE     TRUE
# 125           SENTINEL2                                               Sentinel 2  FALSE FALSE     TRUE
# 126                 SGI                                SGI Image File Format 1.0   TRUE FALSE     TRUE
# 127              SIGDEM                       Scaled Integer Gridded DEM .sigdem  FALSE  TRUE     TRUE
# 128              SNODAS                            Snow Data Assimilation System  FALSE FALSE     TRUE
# 129                 SRP                      Standard Raster Product (ASRP/USRP)  FALSE FALSE     TRUE
# 130             SRTMHGT                                      SRTMHGT File Format  FALSE  TRUE     TRUE
# 131                 TIL                                          EarthWatch .TIL  FALSE FALSE     TRUE
# 132                 TSX                                       TerraSAR-X Product  FALSE FALSE     TRUE
# 133            Terragen                                     Terragen heightfield   TRUE FALSE     TRUE
# 134             USGSDEM                       USGS Optional ASCII DEM (and CDED)  FALSE  TRUE     TRUE
# 135               VICAR                                          MIPL VICAR file   TRUE  TRUE     TRUE
# 136                 VRT                                           Virtual Raster   TRUE  TRUE     TRUE
# 137                 WCS                                 OGC Web Coverage Service  FALSE FALSE     TRUE
# 138                WEBP                                                     WEBP  FALSE  TRUE     TRUE
# 139                 WMS                                      OGC Web Map Service  FALSE  TRUE     TRUE
# 140                WMTS                                 OGC Web Map Tile Service  FALSE  TRUE     TRUE
# 141                 XPM                                        X11 PixMap Format  FALSE  TRUE     TRUE
# 142                 XYZ                                        ASCII Gridded XYZ  FALSE  TRUE     TRUE
# 143                ZMap                                           ZMap Plus Grid  FALSE  TRUE     TRUE
# 144              netCDF                               Network Common Data Format   TRUE  TRUE     TRUE

RasterLayer 클래스 객체를 raster() 함수를 사용해서 처음부터 직접 만들 수도 있습니다.

아래의 예는 8개의 행과 8개의 열을 가진, 총 64 개의 셀(or 픽셀)을 가진 RasterLayer 클래스를 직접 만들어본 것입니다. 이때 레스터 객체의 좌표 참조 시스템(CRS, Coordinates Reference System)은 WGS84 가 기본 설정값이며, 이는 해상도(resolution)의 단위가 도 (in degrees) 라는 의미입니다. 아래 예에서는 res = 0.5 로서 해상도를 0.5도로 설정해주었습니다. 각 셀의 값은 왼쪽 상단부터 시작하여, 행 방향(row-wise)으로 왼쪽에서 오른쪽으로 채워나가게 됩니다. 아래 예는 총 64개 셀의 각 셀의 값을 1~64 의 정수값을 좌측 상단에서 시작하여 우측으로 한 줄씩 채워나가도록 했으며, 제일 아래에 8*8 셀들의 값과 오른쪽에 plot()으로 시각화한 결과를 확인해보시기 바랍니다.

## -- Rasters can also be created from scratch using the raster() function.
## RasterLayer object, The CRS is the default of raster objects: WGS84.
## raster of 64 cells (pixels) = 8 rows * 8 cols (row-wise)
my_raster = raster(nrows = 8, ncols = 8, res = 0.5, 
                   xmn = -2.0, xmx = 2.0, ymn = -2.0, ymx = 2.0,
                   vals = 1:64)


my_raster
# class      : RasterLayer 
# dimensions : 8, 8, 64  (nrow, ncol, ncell)
# resolution : 0.5, 0.5  (x, y)
# extent     : -2, 2, -2, 2  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# source     : memory
# names      : layer 
# values     : 1, 64  (min, max)


## plotting
plot(my_raster, main = "my raster (64 cells = 8 rows * 8 cols)")

(2) RasterBrick Class

위 (1) 의 RasterLayer 클래스가 단지 1개 층 (only 1 layer) 으로만 구성되는 가장 간단한 형태의 레스터 클래스라고 소개했는데요, (2) RasterBrick 클래스와 (3) RasterStack 클래스는 여러개의 층(multiple layers)을 가질 수 있습니다.

특히, RasterBrick 클래스는 단일 다중 스펙트럼 위성 파일 (a single multispectral satellite file) 이나 또는 메모리의 단일 다층 객체 (a single multilayer object in memory) 의 형태로 다층의 레스터 객체를 구성합니다.

아래의 예는 raster 패키지의 brick() 함수를 사용해서 spDataLarge 패키지에 들어있는 landsat.tif 의 다층 레스터 파일을 RasterBrick 클래스 객체로 불러온 것입니다. 

nlayers() 함수로 총 몇 개의 층이 있는지 확인해 보니 4개의 층이 하나의 파일에 들어있네요.

plot() 함수로 RasterBrick 클래스 객체의 4개 층을 모두 시각화보면 아래와 같습니다.

## -- (2) RasterBrick
## : A RasterBrick consists of multiple layers, 
##  which typically correspond to a single multispectral satellite file 
##  or a single multilayer object in memory. 
multi_raster_file = system.file("raster/landsat.tif", package = "spDataLarge")
raster_brick = brick(multi_raster_file)

## 4 layers
raster_brick
# class      : RasterBrick 
# dimensions : 1428, 1128, 1610784, 4  (nrow, ncol, ncell, nlayers)
# resolution : 30, 30  (x, y)
# extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
# source     : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/landsat.tif 
# names      : landsat.1, landsat.2, landsat.3, landsat.4 
# min values :      7550,      6404,      5678,      5252 
# max values :     19071,     22051,     25780,     31961


## nlayers() : number of layers in a Raster* object
nlayers(raster_brick)
# [1] 4


## plotting RasterBrick object with 4 layers
plot(raster_brick)

(3) RasterStack Class

세번째 레스터 클래스인 RasterStack 클래스도 다 층 (multi-layers) 레스터 객체로 구성이 됩니다. 같은 범위와 해상도를 가진 여러개의 RasterLayer 클래스 객체들을 리스트로 묶어서 RasterStack 클래스 객체를 만들 수 있습니다.

이때, 위 (2)번의 RasterBrick 클래스가 동일한 복수개의 RasterLayer 층으로 구성되는 반면에, 이번 (3)번의 RasterStack 클래스는 여러개의 RasterLayer와 RasterBrick 클래스 객체가 혼합되어서 구성할 수 있다는 점입니다.  연산 속도면에서 보면 일반적으로 RasterBrick 클래스가 RasterStack 클래스보다 빠릅니다.

RasterBrick 클래스와 RasterStack 객체에 대한 연산은 보통은 RasterBrack 클래스 객체를 반환합니다.

아래 예에서는 (a) raster(raster_brick, layer = 1) 함수를 사용해서 위의 (2)번에서 불러왔던 RasterBrick 클래스 객체의 1번째 층만 가져다가 raster_on_disk 라는 이름으로 레스터 객체를 하나 만들고, (b) raster() 함수로 동일한 범위와 해상도, 좌표 참조 시스템(CRS)를 가지고 난수로 셀의 값을 채운 raster_in_memory 라는 이름의 메모리에 있는 RasterLayer 클래스 객체를 만들었습니다. 다음에 stac() 함수로 raster_stack = stack(raster_in_memory, raster_on_disk) 처럼 (a) + (b) 하여 쌓아서 raster_stack 라는 이름의 RasterStack 클래스 객체를 만들었습니다. 

마지막으로 plot() 함수로 RasterStack 클래스 객체에 쌓여 있는 2개의 객체를 시각화해보았습니다. (raster_in_memory 는 난수를 발생시켜 셀 값을 채웠기 때문에 시각화했을 때 아무런 패턴이 없습니다.)

## -- (3) RasterStack
## : RasterStack allows you to connect several raster objects 
##   stored in different files or multiple objects in memory.

raster_on_disk = raster(raster_brick, layer = 1)

raster_in_memory = raster(xmn = 301905, xmx = 335745,
                          ymn = 4111245, ymx = 4154085, 
                          res = 30)
values(raster_in_memory) = sample(seq_len(ncell(raster_in_memory)))
crs(raster_in_memory) = crs(raster_on_disk)


## RasterStack is a list of RasterLayer objects with the same extent and resolution.
raster_stack = stack(raster_in_memory, raster_on_disk)
raster_stack
# dimensions : 1428, 1128, 1610784, 2  (nrow, ncol, ncell, nlayers)
# resolution : 30, 30  (x, y)
# extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
# names      :   layer, landsat.1 
# min values :       1,      7550 
# max values : 1610784,     19071


## plotting RasterStack (raster_in_memory + raster_on_disk)
plot(raster_stack)

(4) 언제 어떤 레스터 클래스를 사용하는 것이 좋은가?

위의 (1), (2), (3)에서 소개한 RasterLayer 클래스, RasterBrick 클래스, RasterStack 클래스의 특징과 서로 간의 차이점을 보면 알겠지만 다시 한번 정리하자면요, Raster* 클래스 중에서 무엇을 사용할 지는 투입 데이터 (input data)의 특징이 무엇이냐에 따라 달라집니다.

하나의 다층 레스터 파일이나 객체(a single multilayer file or object)를 처리하는 것이라면 RasterBrick 이 적합합니다.

반면에, 여러개의 레스터 파일들(many files)이나 여러 종류의 레스터 클래스 (many Raster*) 를 한꺼번에 연결해서 연산하고 처리해야 하는 경우라면 RasterStack Class 가 적합하겠습니다.

[ Reference ]

* Geographic data in R : geocompr.robinlovelace.net/spatial-class.html

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

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

지리적 레스터 데이터 (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

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

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

지난번 포스팅에서는 R의 sf 패키지를 사용하여 "sf 클래스 (The sf class)"를 만드는 방법을 소개하였습니다.

이번 포스팅에서는 지리공간 벡터 데이터(GeoSpatial Vector data)를 분석할 때 사용하는 "sf 패키지"와 "sp 패키지"의 관계, "sf 클래스"와 "sp 클래스" 간 변환에 대해서 소개하겠습니다.

(1) R의 sf 패키지와 sp 패키지의 관계/ 역사

(2) sf 클래스를 sp 클래스로 전환하기: as(sf_class, Class="Spatial")

(3) sp 클래스를 sf 클래스로 전환하기: st_as_sf(sp_class)

[ R 지리공간 데이터 분석 패키지: sf package와 sp package 클래스 간 전환 ]

sf 패키지는 {sp 패키지 + R 인터페이스 with (GDAL, GEOS, PROJ) + s2 패키지} 의 패키지 기능들을 승계한 종합 패키지로서, 지리공간 벡터 데이터를 단순 지리특성 (Simple features) 으로 인코딩할 수 있게 해줍니다. [1]

>> sp 패키지 : Classes and Methods for Spatial Data [2], [3]

(sf 패키지의 부모격이 되는) sp 패키지는 2005년 Pebesma 와 Bivand 가 개발해서 처음으로 공개하였습니다. 2005년 전에는 지리적인 좌표가 일반적으로 여느 숫자와 다름없이 처리되었다면, sp 패키지가 나오고 부터는 지리공간 점(Point), 선(LineString), 면/다각형(Polygon), 격자(Grid) 와 속성(Attributes)을 지원하는 클래스(Classes)와 메소드(Methods) 로 처리되는 방식으로 바뀌게 됩니다.

sp 패키지는 테두리 상자(bounding box), 좌표 참조 시스템(CRS, Coordinate Reference System), 속성(Attributes) 등의 정보를 S4 클래스 시스템을 사용해서 "Spatial" 객체안의 슬롯에 저장을 합니다. 이를 통해 지리공간 데이터에 대한 연산 작업을 할 수 있게 해줍니다. 또한, sp 패키지는 지리공간 데이터에 대해 summary()나 plot() 함수와 같은 R 에 내장된 함수도 사용할 수 있게 해줍니다.

>> sf 패키지 : Simple Features for R

sf 패키지는 Edzer Pebesma, Roger Bivand 등이 2016년 10월에 최초로 오픈소스로 공개하였으며, R로 단순 지리특성 기하 (Simple Feature Geometry) 형태로 지리 벡터 데이터를 인코딩하는 표준화된 방법을 지원합니다. sf 패키지는 sp 패키지의 기능을 승계하였으며, 이에 더해 지리공간 데이터를 읽고 쓰는 'GDAL', 지리적 연산을 할 때 사용하는 'GEOS', 지도의 투영 변환(projection conversions)과 데이터 변환(datum transformations)을 위한 'PROJ' 와 R과의 인터페이스를 제공합니다. 그리고 선택적으로 지리적 좌표에 대한 구면 기하 연산 (spherical geometry operations) 을 위해 's2' 패키지를 사용합니다.

단순 지리특성 모델 (Simple Features Model)을 지원하는 "sf 패키지"를 사용하면 좋은 점들로는[4],

• 지리공간 벡터 데이터를 빠르게 읽고 쓸 수 있음
• 지리공간 벡터 데이터 시각화 성능의 고도화
  (tmap, leaflet, mapview 지리공간 데이터 시각화 패키지가 sf 클래스 지원)
  
• 대부분의 연산에서 sf 객체는 DataFrame 처럼 처리가 가능함
• sf 함수들은 '%>%' 연산자 (chain operator) 와 함께 사용할 수 있고, R의 tidyverse 패키지들과도 잘 작동함
  (sp 패키지도 spdplyr 패키지를 설치하면 dplyr의  %>% 체인 연산자와 기능을 사용할 수 있음)
  
• sf 함수이름은 'st_' 로 시작하여 상대적으로 일관성이 있고 직관적임
  

등을 들 수 있습니다.

sf 패키지의 장점이 이렇게 많으므로 지리공간 벡터 데이터를 처리하고 분석하고 시각화하는데 있어 sf 패키지를 사용하지 않을 이유가 없습니다!

예제로 spData 패키지에 들어있는 'world' 라는 이름의 "sf", "data.frame" 클래스의 데이터셋을 사용하겠습니다.

"sf" 클래스를 "sp" 클래스로 변환(converting the sf class into the sp class) 할 때는 as(sf_class, Class = "Spatial") 함수와 매개변수를 사용합니다.(sp 패키지는 S4의 Spatial 클래스를 사용)

반대로, "sp" 클래스를 "sf" 클래스로 변환(converting the sp class into sf class) 할 때는 sf 패키지의 st_as_sf() 함수를 사용합니다.

[ Reference ]

[1] sf package: https://cran.r-project.org/web/packages/sf/index.html

[2] sp package: https://cran.r-project.org/web/packages/sp/index.html

[3] The history of R-spatial: https://geocompr.robinlovelace.net/intro.html

[4] Why simple features?: https://geocompr.robinlovelace.net/spatial-class.html

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

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

지난번 포스팅에서는 R의 sf 패키지를 사용해서

- 벡터 데이터의 기하 유형 (Geometry types of vector data)을 단순 지리특성 기하로 정의하기

- 벡터 데이터의 여러개 단순 지리특성 기하(sfg)를 하나의 단순 지리특성 칼럼(sfc) 객체로 합치기

- st_sfc() 함수로 좌표 참조 시스템(CRS, Coordinate Reference Systems) 을 설정하고 확인하기

에 대해서 소개하였습니다.

이번 포스팅에서는 앞서의 포스팅들에서 소개했던 내용을 모두 포괄하는 sf패키지의 "sf 클래스 (The sf class)" 를 만드는 방법을 소개하겠습니다.

(1) sf 패키지의 "sf 클래스(the sf class)"는 무엇인가?

(2) sf 클래스 만들기: st_sf()

(3) sf 클래스의 클래스, 기하유형, CRS 확인하기: class(), st_geometry_type(), st_crs)

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-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 을 설정해주겠습니다.

(2-2) 비기하 data.frame 속성 (data.frame Non-geometry Attributes)

이번에는 서울의 COEX, GBC(현대자동차 Global Business Center), JansilStadium(잠실 종합 운동장) 의 3개 POI (Point of Interest)별 면적(area)과 가격(price) 속성(Attributes) 정보를 가지는 data.frame을 만들어보겠습니다. (아래 속성 데이터 값은 모두 그냥 가짜로 입력한 예시 값입니다.)

(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) 를 만들어보겠습니다.

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

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

지난번 포스팅에서는 R의 sf패키지의 st_sfc() 함수를 사용해서 여러개의 기본 지리특성 기하(sfg, Simple feature geometry)들을 하나의 단순 지리특성 칼럼(sfc, Simple feature columns) 객체로 합치는 방법을 소개하였습니다.

지리공간 데이터 분석을 할 때 제일 먼저 해야 하고 또 확인해야 하는 것이 좌표계, 좌표 참조 시스템(Coordinate Reference Systems, CRS)을 설정하고 확인하는 것입니다.

좌표 참조 시스템(CRS) (또는 좌표계, 공간 참조 시스템, Spatial Reference System, SRS) 은 좌표 기반으로 국소, 지역, 또는 전지구적으로 지리상의 객체를 위치시킬 때 사용되는 참조 시스템입니다. 좌표 참조 시스템은 특정 지도 투영법 (map projection)을 정의하고, 서로 다른 좌표 참조 시스템 간의 변환을 할 수 있도록 해줍니다.

위도, 경도 좌표의 기준이 어떻게 되는지, 누가 언제 정의한 좌표 체계인지를 정의해 놓은 것입니다.

오픈소스 PostGIS 에서는 Spatial Reference System Table에 SRID 별로 'auth_name', 'auth_srid', 'srtext', 'proj4text' 등의 정보가 저장되어 있는데요, R에서도 좌표계, 좌표 참조 시스템(CRS)이 있어서 좌표계 간에 변환을 할 수 있게 해줍니다.

참고로, 전세계를 한번에 나타내야 할 때 많이 사용하는 전지구 좌표계로는 WGS84(EPSG:4326), Bessel 1841(EPSG:4004, EPSG:4162), GRS80(EPSG)4019, EPGS:4737), Google Mercator(EPSG:3857) 등이 있습니다.

[전지구 좌표계]
*WGS84 경위도: GPS가 사용하는 좌표계EPSG:4326, EPSG:4166 (Korean 1995)+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs 
*Bessel 1841 경위도: 한국과 일본에 잘 맞는 지역타원체를 사용한 좌표계EPSG:4004, EPSG:4162 (Korean 1985)+proj=longlat +ellps=bessel +no_defs +towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43
*GRS80 경위도: WGS84와 거의 유사EPSG:4019, EPSG:4737 (Korean 2000)+proj=longlat +ellps=GRS80 +no_defs
*Google Mercator: 구글지도/빙지도/야후지도/OSM 등 에서 사용중인 좌표계EPSG:3857(공식), 

EPSG:900913(통칭)

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

출처: https://www.osgeo.kr/17 [OSGeo(Open Source GeoSpatial) 한국어 지부 - OSGeo Korean Chapter]

한국 주요 좌표계 EPSG 코드 및 proj4 인자 정리는 https://www.osgeo.kr/17 를 참고하세요.

이번 포스팅에서는 R의 sf패키지를 사용해서 단순 지리특성 칼럼(sfc, Simple feature columns)에서

(1) 특정 좌표 참조 시스템(CRS, Coordinate Reference Systems) 지정해주기

(2) 좌표 참조 시스템(CRS, Coordinate Reference Systems) 확인: st_crs()

(3) raw proj4string 을 사용해서 좌표 참조 시스템 지정해주기

를 소개하겠습니다.

기본 지리특성 칼럼(sfc, Simple feature columns)의 좌표 참조 시스템의 기본 값은 NA (Not Available) 입니다.

st_sfc() 함수의 crs 매개변수를 사용해서 특정 좌표 참조 시스템(CRS)을 설정해줄 수 있습니다.

가령, 아래 예에서는 EPSG 코드 4326 번으로 좌표 참조 시스템을 설정하여 보았습니다.

* EPSG 코드의 EPSG는 European Petroleum Survey Group 의 앞글자를 딴 줄임말입니다. EPSG 그룹은 좌표 시스템 정보 데이터베이스와 지도 투영법과 데이터에 대한 표준 문서 자료를 발행하는 곳입니다.

* WGS84 는 World Geodetic System 1984 의 앞글자를 딴 줄임말입니다. WGS84 는 전지구 좌표계로서, 우리가 매일 사용하는 GPS(Global Positioning System), 지도 제작, 측지학 등에서 사용합니다.

좌표 참조 시스템 확인은 st_crs() 함수를 사용합니다.

raw proj4string 을 직접 입력해서 좌표 참조 시스템을 지정해주려면 st_sfc() 함수의 crs 매개변수에 직접 따옴표 안에 proj, datum 을 입력해주면 됩니다. 아래 예는 WGS84를 raw proj4string로 해서 좌표 참조 시스템을 등록한 것입니다.

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

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

지난번 포스팅에서는 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) 객체로 만들기

R의 sf 패키지를 사용할 것이므로, 먼저 library(sf) 로 패키지를 불러옵니다.

다음으로, 예제로 사용할 점(Point) 두개를 st_point() 함수를 사용해서 만들었습니다.

그리고, st_sfc() 함수를 사용하여, 단순 지리특성 기하(Simple feature geometry, sfg)인 두개의 점을

--> 단순 지리특성 칼럼(Simple feature columns, sfc) 의 하나의 객체로 합쳤습니다.

st_geometry_type() 함수를 사용하면 기하 유형을 확인해 볼 수 있습니다. 동일한 기하 유형인 두 개 점을 합쳤으므로 하나의 객체 집합에 두 개의 점(Points)가 들어있습니다.

이번에는 동일한 기하 형태인 두개의 단순 지리특성 기하 면 (2 Simple feature geometry polygons) 을 st_polygon() 함수로 만들어보겠습니다.

그리고 이를 st_sfc() 함수를 사용해서 하나의 단순 지리특성 칼럼(Simple feature column) 객체로 합쳐보겠습니다.

st_geometry_type()함수로 기하 유형을 확인해보면 2개 면의 집합인 하나의 객체임을 확인할 수 있습니다.

이번에는 (1)번에서 만들었던 점(Point)과 (2)번에서 만들었던 면(Polygon)을 st_sfc() 함수를 사용해서 하나의 단순 지리특성 칼럼 객체로 합쳐보겠습니다. 

이처럼 서로 다른 기하 형태도 st_sfc() 함수로 하나의 단순 지리특성 칼럼 객체로 합칠 수 있습니다.

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

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

이번 포스팅에서는 R의 leaflet 패키지와 Open Street Map을 활용하여 RStudio 또는 Jupyter Notebook R kernel IDE 에서 Interactive Map 을 그리는데 있어,

(1) 기본 상호작용하는 지도 생성하기 (create a basic interactive map)

(2) 지도 스타일 추가하기 (add map styles)

(3) 지도 좌표에 표시 풍선 추가하기 (add marker popup at latitude/longitude coordinates)

(4) 상호작용하는 동적 지도를 웹 페이지로 저장하기 (save an interactive map as web page)

하는 방법을 소개하겠습니다.

아래의 포스팅 예시는 RStudio IDE 에서 진행하였습니다. (Jupyter Notebook R kernel 에서도 가능)

먼저 R의 leaflet 패키지를 설치하고 불러오기를 해보겠습니다.

다음으로 leaflet 패키지의 leaflet() 함수를 사용하여, 기본 OpenStreetMap 타일을 추가하고 (add default OpenStreetMap map tiles), 대한민국의 위도(latitude), 경도(longitude) 위치에 zoom = 6 의 축소/확대 비율로 설정해서 동적인 지도(interactive map)를 그려보겠습니다.

마우스 커서로 좌측 상단의 '+'를 클릭하면 지도가 'zoom in' 확대되어 더 상세하게 볼 수 있으며, 반대로 좌측 상단의 '-' 를 클릭하면 지도가 'zoom out' 축소되어 더 넓은 지역을 볼 수 있게 됩니다. (마우스 휠을 사용해서 zoom in, zoom out 할 수도 있습니다.)

애초에 기본 지도를 생성할 때 지도의 zoom 수준을 setView(zoom = xx) 에서 숫자로 설정해줄 수도 있습니다. 바로 위에서 생성했던 지도가 zoom = 6 이었다면, 아래 지도는 zoom = 8, zoom = 16 으로 기본 지도를 생성하게 해보겠습니다. 아래에 보면, 지도가 좀더 zoom in 으로 확대되어 더 상세하게 볼 수 있습니다.

(물론 interactive map 이기 때문에 마우스를 이용해서 언제든지 자유롭게 지도를 동적으로 zoom in, zoom out 할 수 있습니다.)

addProviderTiles("Tile Name Here") 를 이용하여 외부 지도 타일을 추가할 수 있습니다.

아래는 NASA에서 제공하는 우주의 인공위성에서 바라본 2012년 밤의 지도 타일("NASAGIBS.ViirsEarthAtNight2012") 입니다. 별도의 파일을 미리 다운로드 하거나 할 필요는 없구요, 그냥 타일 이름을 addProviderTiles() 메소드 안에 큰 따옴표로 해서 적어주기만 하면 됩니다. (물론, 외부와 통신할 수 있는 인터넷망에 접속해 있어야 합니다.)

한국의 경우 남한은 밤에 환한 반면에, 북한은 평양 빼고는 어둡네요. 38선 경계도 불빛으로 환하게 구분이 되는 모습을 보노라니 약간 서글프네요.

아래는 zoom = 16 으로 zoom in 해서 삼성동 COEX 부근을 확대해본 것인데요, 여기에 addProviderTiles("Esri.WorldImagery") 타일을 추가해 보았습니다. 항공사진 모습의 타일이네요.

위의 두 예 외에 leaflet() 에 다른 외부 타일 제공자의 타일 지도를 추가하고 싶으면 아래의 url 의 Leaflet-prodivers preview 에 방문해서 원하는 형태의 지도 스타일을 찾았다면 그 이름을 addProviderTiles("Tile Name Here") 메소드 안에 타일 이름을 큰 따옴표 안에 써주기만 하면 됩니다.

https://leaflet-extras.github.io/leaflet-providers/preview/

addMarkers() 메소드를 사용하면 위도(latitude), 경도(longitude) 좌표 위치에 풍선 모양의 표식과 커서를 클릭했을 때 팝업으로 나타나는 설명을 추가할 수 있습니다.

아래 예제에서는 대한민국 삼성동의 COEX Mall (코엑스 몰), GBC (현대자동차 Global Business Center), Samsil-Stadium (올림픽 주 경기장) 의 3개 주요 POI (Point Of Interest) 지점에 말풍선 표식을 추가해보았습니다. 마우스 커서를 말풍선 표식에 가져가서 클릭을 해보면 popup 설명이 튀어나옵니다.

(아래 예제는 스크린 캡쳐를 한 이미지 파일이어서 interactive map 기능이 안됩니다. RStudio 나 Jupyter Notebook R kernel 에서 interactive map 기능 사용 가능합니다.)

위에서 만든 Interactive Map 을 HTML 웹 페이지로 저장하여 내보내기를 할 수 있습니다.

RStudio의 우측 하단의 Viewer 화면에서

   > Export

      > Save as Web Page

메뉴를 선택한 후에

아래와 같은 'Save As Web Page' 저장하기 팝업 화면이 나오면 HTML 파일의 이름과 저장할 경로를 지정해주면 됩니다.  나중에 이 HTML 파일을 웹 페이지에 올려서 바로 Interactive Map 을 바로 웹 페이지에 올릴 수 있습니다.

2020년 마지막 포스팅이 되겠네요.

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

행복한 데이터 과학자 되세요.

그리고 2020년 한해 마무리 잘하시고, 2021년 새해에는 건강하시고 복 많이 받으세요.

지난번 포스팅에서는 지리공간 데이터의 포맷으로서 벡터 데이터(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 으로 나누어준 것은 시각화한 지도 크기에 맞추어서 너무 크지도, 작지도 않은 원을 덮어쓰우기 위해 숫자를 조정한 것입니다.)

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

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

지리공간 데이터 (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) ]

벡터 데이터에는 KML(.kml or .kmz), GML, GeoJSON, Shapefile (.shp), WKT 등의 데이터 포맷이 있습니다.

KML (Keyhole Markup Language), GML (Geography Markup Language) 데이터 포맷은 XML 기반으로 지리공간 데이터를 저장합니다. KML은 OGC(Open Geospatil Consortium)의 공식 표준입니다. KML과 GML 데이터 포맷은 non-GIS 사용자들과 인터넷을 통해 쉽게 지리공간 데이터를 공유하는데 많이 사용됩니다.

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) 비교 ]

* source: https://blog.rmotr.com/spatial-data-with-python-lets-begin-e29b5c41ead3

레스터 데이터(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 패키지로 시각화하는 간단한 예를 소개하겠습니다.

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

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