Raster07：Extract-NDVI-From-Rasters-in-R

最新推荐文章于 2021-11-16 11:07:46 发布

GeoSuper

最新推荐文章于 2021-11-16 11:07:46 发布

阅读量262

点赞数

分类专栏： R

R 专栏收录该内容

27 篇文章 41 订阅

订阅专栏

library(pacman)
p_load(raster, rgdal, ggplot2)
wd <- "G:/Rdata/neon_data/NEONDSLandsatNDVI/NEON-DS-Landsat-NDVI/"
setwd(wd)

#create list of NDVI file paths
all_HARV_NDVI <- list.files(paste0(wd, "HARV/2011/NDVI"),
                                   full.names = TRUE,
                                   pattern = ".tif$")
all_HARV_NDVI

#Create a time series raster stack
NDVI_HARV_stack <- stack(all_HARV_NDVI)

#apply scale factor
NDVI_HARV_stack <- NDVI_HARV_stack*0.0001

#calculate mean NDVI for each raster
avg_NDVI_HARV <- cellStats(NDVI_HARV_stack,mean)

#Convert output array to data.frame
avg_NDVI_HARV <- as.data.frame(avg_NDVI_HARV)
avg_NDVI_HARV

#                          avg_NDVI_HARV
# X005_HARV_ndvi_crop      0.365150
# X037_HARV_ndvi_crop      0.242645
# X085_HARV_ndvi_crop      0.251390
# X133_HARV_ndvi_crop      0.599300
# X181_HARV_ndvi_crop      0.878725
# X197_HARV_ndvi_crop      0.893250
# X213_HARV_ndvi_crop      0.878395
# X229_HARV_ndvi_crop      0.881505
# X245_HARV_ndvi_crop      0.850120
# X261_HARV_ndvi_crop      0.796360
# X277_HARV_ndvi_crop      0.033050
# X293_HARV_ndvi_crop      0.056895
# X309_HARV_ndvi_crop      0.541130


#view only the value in row 1, column 1 of the data frame

avg_NDVI_HARV[1,1]
# [1] 0.36515
avg_NDVI_HARV[2,1]
# [1] 0.242645

#view column name slot
names(avg_NDVI_HARV)
# [1] "avg_NDVI_HARV"

#rename the NDVI column
names(avg_NDVI_HARV) <- "meanNDVI"
names(avg_NDVI_HARV)
# [1] "meanNDVI"

#add a site column to our data
avg_NDVI_HARV$site <- "HARV"

#add a year column to our data
avg_NDVI_HARV$year <- "2011"

#view data
avg_NDVI_HARV

#                     meanNDVI site year
# X005_HARV_ndvi_crop 0.365150 HARV 2011
# X037_HARV_ndvi_crop 0.242645 HARV 2011
# X085_HARV_ndvi_crop 0.251390 HARV 2011
# X133_HARV_ndvi_crop 0.599300 HARV 2011
# X181_HARV_ndvi_crop 0.878725 HARV 2011
# X197_HARV_ndvi_crop 0.893250 HARV 2011
# X213_HARV_ndvi_crop 0.878395 HARV 2011
# X229_HARV_ndvi_crop 0.881505 HARV 2011
# X245_HARV_ndvi_crop 0.850120 HARV 2011
# X261_HARV_ndvi_crop 0.796360 HARV 2011
# X277_HARV_ndvi_crop 0.033050 HARV 2011
# X293_HARV_ndvi_crop 0.056895 HARV 2011
# X309_HARV_ndvi_crop 0.541130 HARV 2011

#Extract Julian Day from row.names
JulianDays <- gsub(pattern = "X|_HARV_ndvi_crop", #the pattern to find
               x = row.names(avg_NDVI_HARV), #the object containing the strings
               replacement = "")
# [1] "005" "037" "085" "133" "181" "197" "213" "229" "245" "261" "277" "293" "309"
#alternately you can include the above code on one single line
#JulianDays <- gsub("X|_HARV_NDVI_crop","",row.names(avg_NDVI_HARV))



avg_NDVI_HARV$JulianDay <- JulianDays
avg_NDVI_HARV

#what class is the new column
class(avg_NDVI_HARV$JulianDay)
# [1] "character"

# convert Julian Day to Date Class
#set the origin for the julian date(1 Jan 2011)
origin <- as.Date("2011-01-01")

#convert "JulianDay" from class character to integer
avg_NDVI_HARV$JulianDay <-as.integer(avg_NDVI_HARV$JulianDay)
class(avg_NDVI_HARV$JulianDay)
# [1] "integer"

#create a date column:-1 added because origin is the 1st
#if not -1,01/01/2011 + 5 = 06/01/2011  which is Julian day6, not 5
avg_NDVI_HARV$Date <- origin + (avg_NDVI_HARV$JulianDay-1)
avg_NDVI_HARV
class(avg_NDVI_HARV$Date)
# [1] "Date"

#plot NDVI
ggplot(avg_NDVI_HARV,aes(JulianDay,meanNDVI),na.rm = TRUE)+
  geom_point(size = 4,colour = "PeachPuff4")+
  ggtitle("Landsat Derived NDVI-2011\nNEON Harvard Forest Field Site"
          )+
  xlab("Julian Days")+ylab("Mean NDVI")+
  theme(text=element_text(size = 20))

在这里插入图片描述

#load the SJER data
all_NDVI_SJER <- list.files(paste0(wd,"SJER/2011/NDVI"),
                            full.names = TRUE,
                            pattern = ".tif$")

NDVI_SJER_stack <- stack(all_NDVI_SJER)
avg_NDVI_SJER <- cellStats(NDVI_SJER_stack,mean)
avg_NDVI_SJER <- as.data.frame(avg_NDVI_SJER)
avg_NDVI_SJER
names(avg_NDVI_SJER) <-"MeanNDVI"
avg_NDVI_SJERa
avg_NDVI_SJER$site <- "Site"
avg_NDVI_SJER$year <- "year"
JulianDay1 <- gsub("X|_SJER_ndvi_crop","",row.names(avg_NDVI_SJER))
class(JulianDay1)
JulianDay1 <- as.integer(JulianDay1)
class(JulianDay1)

avg_NDVI_SJER$JulianDay <- JulianDay1
class(avg_NDVI_SJER$JulianDay)


Date <- origin + (avg_NDVI_SJER$JulianDay-1)
avg_NDVI_SJER$Date <- Date
avg_NDVI_SJER
ggplot(avg_NDVI_SJER,aes(JulianDay,MeanNDVI),na.rm = TRUE)+
  geom_point(size = 4,colour = "SpringGreen4")+
  ggtitle("Landsat Derived NDVI -2011\nNEON SJER Field Site")+
  xlab("Julian Days")+ylab("Mean NDVI")+
  theme(text = element_text(size = 20))

在这里插入图片描述

#Compare NDVI from Two Different Sites in one Plot
#Merge Data Frames
avg_NDVI_HARV
avg_NDVI_SJER
NDVI_HARV_SJER <- rbind(avg_NDVI_HARV,avg_NDVI_SJER)
NDVI_HARV_SJER
# meanNDVI site year JulianDay       Date
# X005_HARV_ndvi_crop    0.365150 HARV 2011         5 2011-01-05
# X037_HARV_ndvi_crop    0.242645 HARV 2011        37 2011-02-06
# X085_HARV_ndvi_crop    0.251390 HARV 2011        85 2011-03-26
# X133_HARV_ndvi_crop    0.599300 HARV 2011       133 2011-05-13
# X181_HARV_ndvi_crop    0.878725 HARV 2011       181 2011-06-30
# X197_HARV_ndvi_crop    0.893250 HARV 2011       197 2011-07-16
# X213_HARV_ndvi_crop    0.878395 HARV 2011       213 2011-08-01
# X229_HARV_ndvi_crop    0.881505 HARV 2011       229 2011-08-17
# X245_HARV_ndvi_crop    0.850120 HARV 2011       245 2011-09-02
# X261_HARV_ndvi_crop    0.796360 HARV 2011       261 2011-09-18
# X277_HARV_ndvi_crop    0.033050 HARV 2011       277 2011-10-04
# X293_HARV_ndvi_crop    0.056895 HARV 2011       293 2011-10-20
# X309_HARV_ndvi_crop    0.541130 HARV 2011       309 2011-11-05
# X014_SJER_ndvi_crop  482.160000 SJER 2011        14 2011-01-14
# X046_SJER_ndvi_crop 5297.800000 SJER 2011        46 2011-02-15
# X062_SJER_ndvi_crop 5543.680000 SJER 2011        62 2011-03-03
# X094_SJER_ndvi_crop 6010.960000 SJER 2011        94 2011-04-04
# X110_SJER_ndvi_crop 5558.360000 SJER 2011       110 2011-04-20
# X126_SJER_ndvi_crop 5383.360000 SJER 2011       126 2011-05-06
# X142_SJER_ndvi_crop 4008.680000 SJER 2011       142 2011-05-22
# X158_SJER_ndvi_crop 4452.400000 SJER 2011       158 2011-06-07
# X174_SJER_ndvi_crop 4074.920000 SJER 2011       174 2011-06-23
# X190_SJER_ndvi_crop 4015.240000 SJER 2011       190 2011-07-09
# X206_SJER_ndvi_crop 3742.600000 SJER 2011       206 2011-07-25
# X222_SJER_ndvi_crop 3670.520000 SJER 2011       222 2011-08-10
# X238_SJER_ndvi_crop 3558.960000 SJER 2011       238 2011-08-26
# X254_SJER_ndvi_crop  270.920000 SJER 2011       254 2011-09-11
# X270_SJER_ndvi_crop 3493.080000 SJER 2011       270 2011-09-27
# X286_SJER_ndvi_crop 4390.680000 SJER 2011       286 2011-10-13
# X302_SJER_ndvi_crop 4683.840000 SJER 2011       302 2011-10-29

#plot NDVI values for both sites
ggplot(NDVI_HARV_SJER,aes(JulianDay,meanNDVI,color = site))+
  geom_point(size = 4, aes(group=site))+
  geom_line(aes(group = site))+
  ggtitle("Landsat Derived NDVI-2011\n Harv vs SanJ\nNEON field Sites")+
  xlab("Julian Day")+ylab("Mean NDVI")+
  scale_colour_manual(values=c("PeachPuff4","SpringGreen4"))+
  theme(text=element_text(size = 20))

在这里插入图片描述

#Remove outlier Date
#open up RGB Image
all_rgb <-list.files(paste0(wd,"HARV/2011/RGB"),
                     full.names = TRUE,
                     pattern = ".tif$")

#create a layout
par(mfrow = c(4,4))


for (aimg in all_rgb){
  RGB_stack <- stack(aimg)
  plotRGB(RGB_stack,stretch="lin")

在这里插入图片描述

par(mfrow = c(4,4))
for (afile in all_rgb_SJER){
  RGB_SJER_stack <- stack(afile)
  plotRGB(RGB_SJER_stack,stretch="lin")
}

在这里插入图片描述

#retain only rows with meanNDVI>0.1
avg_NDVI_HARV_clean <- subset(avg_NDVI_HARV,meanNDVI>0.1)
avg_NDVI_HARV$meanNDVI<0.1
avg_NDVI_HARV_clean$meanNDVI<0.1
#plot without questionable data
ggplot(avg_NDVI_HARV_clean,aes(JulianDay,meanNDVI))+
  geom_point(size = 4,colour = "SpringGreen4")+
  ggtitle("Landsat Derived NDVI-2011\n NEON Harvard Forest Field Site")+
  xlab("Julian Days")+ylab("Mean NDVI")+
  theme(text = element_text(size=20))

在这里插入图片描述
导出为csv

#create a new df to prevent changes to avg_NDVI_HARV
NDVI_HARV_toWrite <- avg_NDVI_HARV_clean

#drop the row.names column
row.names(NDVI_HARV_toWrite)<-NULL
NDVI_HARV_toWrite
# 1  0.365150 HARV 2011         5 2011-01-05
# 2  0.242645 HARV 2011        37 2011-02-06
# 3  0.251390 HARV 2011        85 2011-03-26
# 4  0.599300 HARV 2011       133 2011-05-13
# 5  0.878725 HARV 2011       181 2011-06-30
# 6  0.893250 HARV 2011       197 2011-07-16
# 7  0.878395 HARV 2011       213 2011-08-01
# 8  0.881505 HARV 2011       229 2011-08-17
# 9  0.850120 HARV 2011       245 2011-09-02
# 10 0.796360 HARV 2011       261 2011-09-18
# 11 0.541130 HARV 2011       309 2011-11-05

#create a csv of meanNDVI
write.csv(NDVI_HARV_toWrite,file = paste0(wd,"meanNDVI_HARV_2011.csv"))

GeoSuper

关注

0
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
Raster07：Extract-NDVI-From-Rasters-in-R

library(pacman)p_load(raster, rgdal, ggplot2)wd <- "G:/Rdata/neon_data/NEONDSLandsatNDVI/NEON-DS-Landsat-NDVI/"setwd(wd)#create list of NDVI file pathsall_HARV_NDVI <- list.files(paste0(wd, "HARV/2011/NDVI"),
复制链接

扫一扫

专栏目录