R语言实现ndvi计算

Wuren

基于terra包，先放代码：
library(pacman)
p_load(sp,terra,fs,sf,tidyverse)
b1<-rast(".../LC08_L2SP_122043_20221224_20230103_02_T1_SR_B1.TIF")
> b2<-rast(".../LC08_L2SP_122043_20221224_20230103_02_T1_SR_B2.TIF")
> b3<-rast(".../LC08_L2SP_122043_20221224_20230103_02_T1_SR_B3.TIF")
> b4<-rast(".../LC08_L2SP_122043_20221224_20230103_02_T1_SR_B4.TIF")
> b5<-rast(".../LC08_L2SP_122043_20221224_20230103_02_T1_SR_B5.TIF")
> b6<-rast(".../LC08_L2SP_122043_20221224_20230103_02_T1_SR_B6.TIF")
> b7<-rast(".../LC08_L2SP_122043_20221224_20230103_02_T1_SR_B7.TIF")
pre<-c(b1,b2,b3,b4,b5,b6,b7)
names(pre)<-c('coa','blue','green','red','nir','swir1','swir2')##基础数据准备，加载影像，命名波段，实测landset8/9和sentinel2都可以这样读取。我这下载的数据是预先去云和校准了的，R里怎么搞些我没找到包。

fg<-vect(".../xian.shp")
pre_fg<-crop(pre_crop,fg,mask=T)##按掩膜提取，已经预先把.shp的坐标系转为WGS 1984 UTM Zone 49N。

mndvi<-function(nir,red){
ndvi<-(nir-red)/(nor+red)
return(ndvi)}##定义ndvi计算方法。

pre_fg_ndvi<-mndvi(pre_fg$nir,pre_fg$red)##计算ndvi

quantile(as.array(pre_fg_ndvi),probs=seq(0,1,0.1),T)##提取栅格的0.1分位数，要将栅格转为数组
hist(pre_fg_ndvi,breaks=40)

rr<-c(-0.1,0.24,1,0.24,0.32,2,0.32,0.6,3)
rm<-matrix(rr,ncol=3,byrow=T)
fg_ndvi<-classify(pre_fg_ndvi,rm,include.lowest=TRUE)##栅格重分类,-0.1~0.24定义为1，依此类推，区间瞎选的。

writeRaster(fg_ndvi,".../fg_ndvi.tif",overwrite=T)##输出栅格

AEFvgb

感谢分享

毛毛虫、

谢谢分享