利用raster包进行栅格数据处理（1）

，欢迎也在微信上查看相关内容。

这系列开始介绍raster包对遥感数据的处理。与terra包相比，raster包虽然处理速度较慢，但功能更多、更稳定。语法上面，也各有有缺点，大家可以通过动过代码去体会。

     这次内容主要为以下部分：

•  栅格计算

•  裁剪

•  重投影和重采用

    了解下RasterStack对象

 它的本质是具有相同空间范围、分辨率和坐标参考系统(CRS)的栅格图层集合。

   数据下载链接：

https://raw.githubusercontent.com/joaofgoncalves/R_exercises_raster_tutorial/master/data/LT8_PNPG.zip

 library(raster)library(sp)library(rgdal)setwd("F:/test")#数据路径，lansat8fp <- list.files(path ="./data-raw",                 pattern = ".tif$",full.names = TRUE)fprst <- stack(fp)print(rst)#重命名names(rst) <-paste("b",1:7,sep="")#显示plot(rst)

栅格计算

以计算NDVI为例：NDVI = (NIR – Red) / (NIR + red)

#-----栅格计算---

#方法1,由于处理是在内存中同时进行

#须确保数据能够装入RAM

ndvi<-(rst[[5]]-rst[[4]])/(rst[[5]]+rst[[4]])

#方法2，利用calc函数.对于大型对象，calc将通过栅格块计算值，

#从而节省内存。

calcNDVI_1 <- function(x) return((x[[5]] - x[[4]]) / (x[[5]] + x[[4]]))

ndvi1 <- calc(rst, fun = calcNDVI_1)

#方法3，利用overlay函数允许组合两个(或更多)raster*对象。

#类似于calc。

calcNDVI_2 <- function(x, y) return((x - y) / (x + y))

ndvi2 <- overlay(x = rst[[5]], y = rst[[4]],

                 fun = calcNDVI_2)

# NDVI的值在-1至1之间

#ndvi值有<-1和>1的异常值

summary(ndvi)

ndvi[ndvi < -1] <- NA

ndvi[ndvi>1]<-NA

plot(ndvi, main="NDVI Peneda-Geres National Park",

     xlab = "X-coordinates", ylab = "Y-coordinates")

#提取NDVI>0.4

ndviMask<-ndvi>0.4

plot(ndviMask,main="NDVI mask")

#裁剪xmin <- 554615xmax <- 589015ymin <- 4618355ymax <- 4654705newExtent <- extent(xmin, xmax, ymin,ymax)cropRst <- crop(rst, newExtent)