用R语言计算Landsat影像NDVI的实验

加载包；

library(pacman)
p_load(raster,sp,rgdal,rgeos,tidyverse)

加载影像，并以真彩色显示；

img <- stack('Anhui/LC08_L2SP_121038_20201024_20201105_02_T1/hfLandsat8.tif')
plotRGB(img,r=3,g=2,b=1,axes=TRUE,stretch='lin',main="Landsat True Colour Composite")
nlayers(img)#查看图层数
res(img)#查看空间分辨率
img
# class      : RasterStack 
# dimensions : 5884, 4034, 23736056, 6  (nrow, ncol, ncell, nlayers)
# resolution : 30, 30  (x, y)
# extent     : 466845, 587865, 3424335, 3600855  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs 
# names      : hfLandsat8.1, hfLandsat8.2, hfLandsat8.3, hfLandsat8.4, hfLandsat8.5, hfLandsat8.6 
# min values :       -31788,       -32711,       -32103,       -32762,       -32745,       -32748 
# max values :        32502,        32294,        32749,        32746,        32749,        32741 

把raster转为data frame

img.df <- as.data.frame(img)#把raster转为dataframe
str(img.df)#查看数据结构
# 'data.frame':	23736056 obs. of  6 variables:
# $ hfLandsat8.1: num  NA NA NA NA NA NA NA NA NA NA ...
# $ hfLandsat8.2: num  NA NA NA NA NA NA NA NA NA NA ...
# $ hfLandsat8.3: num  NA NA NA NA NA NA NA NA NA NA ...
# $ hfLandsat8.4: num  NA NA NA NA NA NA NA NA NA NA ...
# $ hfLandsat8.5: num  NA NA NA NA NA NA NA NA NA NA ...
# $ hfLandsat8.6: num  NA NA NA NA NA NA NA NA NA NA ...
colnames(img.df)<-c('Band1','Band2','Band3','Band4','Band5','Band6')#重命名dataframe的列
head(img.df)#查看前几行
# Band1 Band2 Band3 Band4 Band5 Band6
# 1    NA    NA    NA    NA    NA    NA
# 2    NA    NA    NA    NA    NA    NA
# 3    NA    NA    NA    NA    NA    NA
# 4    NA    NA    NA    NA    NA    NA
# 5    NA    NA    NA    NA    NA    NA
# 6    NA    NA    NA    NA    NA    NA

提取经纬度到dataframe

SpatialPoints(img,CRS("+init=epsg:32650"))%>%
  spTransform(., CRS("+init=epsg:4326"))->lonlat1
head(lonlat1)

lonlat1.df <-as.data.frame(lonlat1,na.rm=TRUE)
head(lonlat1.df)
colnames(lonlat1.df)<-c('Lon','Lat')
head(lonlat1.df)
df<-cbind(img.df,lonlat1.df)
head(df)
# Band1 Band2 Band3 Band4 Band5 Band6      Lon      Lat
# 1    NA    NA    NA    NA    NA    NA 116.6470 32.54444
# 2    NA    NA    NA    NA    NA    NA 116.6474 32.54444
# 3    NA    NA    NA    NA    NA    NA 116.6477 32.54444
# 4    NA    NA    NA    NA    NA    NA 116.6480 32.54444
# 5    NA    NA    NA    NA    NA    NA 116.6483 32.54444
# 6    NA    NA    NA    NA    NA    NA 116.6486 32.54444

计算NDVI

NDVICalc <- function(x){
  ndvi<-(x[[4]]-x[[3]])/(x[[3]]+x[[4]])
  return(ndvi)
}
NDVI <- NDVICalc(img)
plot(NDVI)

很明显上面计算的NDVI中有异常值，我们NDVI应该是在-1至1之间的；怎么办呢？我也不知道。。。Google下吧！
果然找到了解决方法,
先看下NDVI的直方图，和NDVI信息（values : -31382, 109.7899 (min, max)）

hist(NDVI,main="NDVI_hist",xlab="NDVI_VALUE",ylab="Frequency",col="wheat1")
NDVI
# class      : RasterLayer 
# dimensions : 5884, 4034, 23736056  (nrow, ncol, ncell)
# resolution : 30, 30  (x, y)
# extent     : 466845, 587865, 3424335, 3600855  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs 
# source     : memory
# names      : layer 
# values     : -31382, 109.7899  (min, max)

利用阈值法把-1,1范围外的赋为NA，然后再看下NDVI的原信息
这回，values : -0.8023534, 0.8727696 (min, max)，搞定。

NDVI[NDVI>1|NDVI<=-1]<-NA
NDVI
hist(NDVI)
# class      : RasterLayer 
# dimensions : 5884, 4034, 23736056  (nrow, ncol, ncell)
# resolution : 30, 30  (x, y)
# extent     : 466845, 587865, 3424335, 3600855  (xmin, xmax, ymin, ymax)
# crs        : +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs 
# source     : memory
# names      : layer 
# values     : -0.8023534, 0.8727696  (min, max)

plot(NDVI)

把NDVI存储为dataframe，并统计每列

NDVI.df <- as.data.frame(NDVI)
str(NDVI.df)
colnames(NDVI.df)<-"NDVI"
head(NDVI.df)
allvars <- cbind(df,NDVI.df)
head(allvars)
allvars
allvars <-na.exclude(allvars)
summary(allvars)
# Band1            Band2            Band3            Band4            Band5            Band6       
# Min.   :-31788   Min.   :-32711   Min.   :-32103   Min.   :-29310   Min.   :-32745   Min.   :-32748  
# 1st Qu.:  8786   1st Qu.:  9546   1st Qu.:  9382   1st Qu.: 13686   1st Qu.: 12085   1st Qu.:  9993  
# Median :  9167   Median : 10040   Median : 10089   Median : 15298   Median : 13600   Median : 11164  
# Mean   :  9289   Mean   : 10169   Mean   : 10216   Mean   : 14641   Mean   : 13252   Mean   : 11190  
# 3rd Qu.:  9660   3rd Qu.: 10667   3rd Qu.: 10949   3rd Qu.: 16374   3rd Qu.: 14905   3rd Qu.: 12459  
# Max.   : 32502   Max.   : 32294   Max.   : 29784   Max.   : 32746   Max.   : 32749   Max.   : 32741  
# Lon             Lat             NDVI        
# Min.   :116.6   Min.   :30.95   Min.   :-0.8024  
# 1st Qu.:117.1   1st Qu.:31.53   1st Qu.: 0.1255  
# Median :117.3   Median :31.77   Median : 0.1897  
# Mean   :117.3   Mean   :31.76   Mean   : 0.1693  
# 3rd Qu.:117.5   3rd Qu.:32.01   3rd Qu.: 0.2488  
# Max.   :117.9   Max.   :32.54   Max.   : 0.8728