空间数据分析与r语言实践_R语言空间数据分析（十）：简单栅格计算—

作者：黄天元，复旦大学博士在读，热爱数据科学与开源工具（R），致力于利用数据科学迅速积累行业经验优势和科学知识发现，涉猎内容包括但不限于信息计量、机器学习、数据可视化、应用统计建模、知识图谱等，著有《R语言高效数据处理指南》（《R语言数据高效处理指南》(黄天元)【摘要书评试读】- 京东图书）。知乎专栏：R语言数据挖掘。邮箱：huang.tian-yuan@qq.com.欢迎合作交流。

HopeR：R语言空间数据分析（零）：总目录zhuanlan.zhihu.com

NDVI的全称是Normalized Difference Vegetation Index，即归一化植被指数，是反映农作物长势和营养信息的重要参数之一，计算方法是(NIR-R)/(NIR+R)。具体介绍可见姚立坤：NDVI - Normalized Differnce Vegetation Index。在R中，我们可以通过定义函数来计算它。

NDVIfun <- function(NIR, Red) {
  NDVI <- (NIR - Red) / (NIR + Red)
  return(NDVI)
}

同时，数据的导入与之前类似（请参考HopeR：R语言空间数据分析（九）：栅格数据探索性分析），过程如下：

library(pacman)
p_load(sp,raster,rgeos,rgdal,rasterVis,raster,fs,sf,tidyverse)

dir("data/LC08_L1TP_203023_20190513_20190521_01_T1/mask",full.names = T) %>% 
  str_subset(pattern = "_B[1-7]_") %>% 
  stack() -> manc

# Name the Bands based on where they sample the electromagentic spectrum
names(manc) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2') 
manc_rgb <- stack(manc$red, manc$green, manc$blue)

然后，调用函数直接求得NDVI：

ndvi <- NDVIfun(manc$NIR, manc$red)

可以简单对其进行绘图：

ndvi %>%
  plot(.,col = rev(terrain.colors(10)), main = "Landsat-NDVI")

其中，col参数可以设置颜色，我们一共设置了10级，并用rev来反转颜色条带。不妨来看一下它具体是什么数值：

同时，我们还可以调查它的取值直方图：

ndvi %>%
  hist(., breaks = 40, main = "NDVI Histogram", xlim = c(-.3,.8))

上面的代码中，xlim可以控制横轴的取值，breaks代表分为40个区间，main则可以声明题目。根据不懂格点的NDVI数值，我们可以判断哪些地方更可能有植被覆盖。比如，我们定义大于3/4分位数的地方，是植被覆盖的地方，那么可以以这个为阈值开始进行植被判定。

veg <- ndvi %>%
  reclassify(., cbind(-Inf, quantile(ndvi,3/4), NA))

veg %>%
  plot(.,main = 'Possible Veg cover')

在reclassify函数中，第二个参数告诉计算机，从负无穷到ndvi的3/4分位数的数字，都转化为NA，那么剩下的就是植被覆盖的地方。

我们可以把曼彻斯特的真彩色图像与这个图层叠加起来进行观察：

manc_rgb %>%
  plotRGB(.,axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

veg %>%
  plot(., add=TRUE, legend=FALSE)

空间数据分析与r语言实践_R语言空间数据分析（十）：简单栅格计算——NDVI