Raster 00: Intro to Raster Data in R

R中的栅格数据介绍：

1、栅格的定义

栅格数据是以二维矩阵的形式来表示空间地物或现象分布的数据组织方式．每个矩阵单位称为一个栅格单元（cell）。栅格的每个数据表示地物或现象的属性数据。因此栅格数据有属性明显，定位隐含的特点。其中栅格数据按数据类型又可以分为连续型的栅格数据（高程、降水、温度）和分类型的栅格数据（分类图（一个像素代表一种地物））。

2、在R中加载并处理栅格

查看栅格数据属性

GeoTIFF包含一组嵌入式标签，其中包含有关栅格的元数据。

library(pacman)
p_load(raster,rgdal,rasterVis,ggplot2)

wd <- "G:/Rdata/neon_data/NEONDSAirborneRemoteSensing/"
setwd(wd)
#view attributes before opening file
GDALinfo(paste0(wd,"HARV/DSM/HARV_dsmCrop.tif"))
# rows        1367 
# columns     1697 
# bands       1 
# lower left origin.x        731453 
# lower left origin.y        4712471 
# res.x       1 
# res.y       1 
# ysign       -1 
# oblique.x   0 
# oblique.y   0 
# driver      GTiff 
# projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
# file        G:/Rdata/neon_data/NEONDSAirborneRemoteSensing/HARV/DSM/HARV_dsmCrop.tif 
# apparent band summary:
#         GDType hasNoDataValue NoDataValue blockSize1 blockSize2
# 1 Float64           TRUE       -9999          1       1697
# apparent band statistics:
#         Bmin   Bmax    Bmean      Bsd
# 1 305.07 416.07 359.8531 17.83169
# Metadata:
#         AREA_OR_POINT=Area 

加载栅格

#load raster into R
DSM_HARV <- raster(paste0(wd,"HARV/DSM/HARV_dsmCrop.tif"))

#view raster structure
DSM_HARV

# class      : RasterLayer
# dimensions : 1367, 1697, 2319799  (nrow, ncol, ncell)
# resolution : 1, 1  (x, y)
# extent     : 731453, 733150, 4712471, 4713838  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
# source     : G:/Rdata/neon_data/NEONDSAirborneRemoteSensing/HARV/DSM/HARV_dsmCrop.tif
# names      : HARV_dsmCrop
# values     : 305.07, 416.07  (min, max)

#plot raster
plot(DSM_HARV,main="NEON Digital Surface Model\nHarvard Forest")

GeoTIFF：

在各种地理信息系统、摄影测量与遥感等应用中, 都要求图像具有地理编码信息, 例如图像所在的坐标系、比例尺、图像上点的坐标、经纬度、长度单位及角度单位等等。对于存储和读取这些信息，纯TIFF 格式的图像文件是很难做到的，而是GeoTIFF 作为TIFF 的一种扩展, 在TIFF 的基础上定义了 一些GeoTag (地理标签)，从而对各种坐标系统、椭球基准、投影信息等进行定义和存储，使图像数据和地理数据存储在同一图像文件中, 这样就为广大用户制作和使用带有地理信息的图像提供了方便的途径。

Coordinate Reference System：

A comprehensive online library of CRS information.
QGIS Documentation - CRS Overview.
Choosing the Right Map Projection.
NCEAS Overview of CRS in R.
坐标参考系统或CRS告诉R栅格在地理空间中的位置。 它还告诉R应该使用什么方法在地理空间中“展平”或投影栅格。

myCRS <- crs(DSM_HARV)
myCRS
# CRS arguments:
#   +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs

范围：

分辨率

计算栅格最小最大值

#calculate Raster Min and Max Values
DSM_HARV <- setMinMax(DSM_HARV)
minValue(DSM_HARV)
# [1] 305.07
maxValue(DSM_HARV)
# [1] 416.07

NoData Values

#NoData Values in rasters
RGB_stack <- stack(paste0(wd, "HARV/RGB_Imagery/HARV_RGB_Ortho.tif"))
#create an RGB image from the raster stack
par(col.axis="white", col.lab="white", tck=0)
plotRGB(RGB_stack, r=1, g=2, b=3,
        axes=TRUE, main="Raster With NoData Values\nRendered in Black")

给NoData 赋值为NA值

func <- function(x){
        x[rowSums(x == 0) == 3, ] <- NA
        x
}
newRGBImage <- calc(RGB_stack,func)
par(col.axis="white", col.lab="white", tck=0)
plotRGB(newRGBImage, r=1, g=2, b=3,
     axes = TRUE, main = "Raster with No Data Values\nNoDataValues = NA"
     )

坏值/错误值/异常值

即超出数据集适用范围的值（例如，0-1（0-10000）范围内反射率数据超过1（10000）的值；NDVI数据中超过[-1，1]范围的值）。
那么如何利用R语言找出这些在数据收集或者处理的过程中出现的异常值呢？
创建直方图

但是直方图默认值处理100000个像素值，而本栅格数据有2319799个像素值
所以：

#create histogram that includes with all pixels values in the raster
hist(DSM_HARV,
     maxpixels = ncell(DSM_HARV),
     main="Distribution of DSM Values\n All Pixels Values Included
     NEON Harvard Forest Field site",
     xlab="DSM Elevation Value(m)",
     ylab="Frequency",
     col="wheat4")

多波段栅格数据

#view number of bands
nlayers(DSM_HARV)
# [1] 1