library(GOESDiurnalNDVI)
#> Warning: replacing previous import 'lubridate::here' by 'plyr::here' when
#> loading 'GOESDiurnalNDVI'Please cite the paper if you use this package: Wheeler, K.I.; Dietze, M.C. A Statistical Model for Estimating Midday NDVI from the Geostationary Operational Environmental Satellite (GOES) 16 and 17. Remote Sens. 2019, 11, 2507.
This is an example of how to use the diurnal fit method explained in Wheeler and Dietze (in press, Remote Sensing) to estimate daily midday, maximum NDVI values from GOES-R (GOES-16 and GOES-17) data. The diurnal model is fit using a Bayesian framework, but it is our intent that as long as the right packages are installed (specifically the packages rjags and runjags), those with limited Bayesian knowledge will be able to use this R package to apply the method. In summary, the model is comprised of two main parts. (1) It assumes that without any noise (e.g., from clouds that were not filtered out by a cloud mask), the “true” NDVI values fall along an inverted “u” shape (negative exponential decrease –actually increases– in the morning and a negative exponential increase in the afternoon – actually decreases).(2) It includes an error model that accounts for the negative bias in the noise by calculating the probability that each observation is clear or cloudy and then the amount that the atmosphere transmits if it is “cloudy.” For a full description, see the paper (open access available).
In this document, we will lay out the general steps to fit the model to one example day (2017-08-21 at Russell Sage, USA). We have included a subset of the data in this folder to allow for the model to fit. We will walk you through the steps from requesting the data through fitting the model and interpretating the outputs (these are also explained in the paper).
1.1 Go to the website and create an account: https://www.avl.class.noaa.gov/saa/products/search?datatype_family=GRABIPRD
1.1 Fill in the desired fields (note: you can only request so much data at a time) 1.1.1 Radiance channel 2 (red): Datatype: “ABI L1b Radiances Data” Satellite: G16 (This package has been specifically built for G16 satellite data and has yet to be tested on G17) ABI Mode: M3 Product Type: Radiances ABI Channel: CO2 ABI Scan Sector: CONUS (This package has been specifically built for the continental US and has yet to be tested on other scan sectors)
1.1.2 Radiance channel 3 (near-infrared): Datatype: “ABI L1b Radiances Data” Satellite: G16 (This package has been specifically built for G16 satellite data and has yet to be tested on G17) ABI Mode: M3 Product Type: Radiances ABI Channel: CO3 ABI Scan Sector: CONUS (This package has been specifically built for the continental US and has yet to be tested on other scan sectors)
1.1.3 Clear Sky Mask (ACM): Datatype: “ABI L2+ Cloud and Moisture Imagery Data”" Satellite: G16 (This package has been specifically built for G16 satellite data and has yet to be tested on G17) ABI Mode: M3 Product Type: Clear Sky Mask ABI Channel:
library(devtools)
#> Loading required package: usethis
#install_github("k-wheeler/NEFI_pheno/GOESDiurnalNDVI")
library(GOESDiurnalNDVI)###Site Characteristics:
siteName <- "RussellSage"
lat <- 32.457
long <- -91.9743
year <- 2017
TZ <- 5 #Time zone
day <- 234 #Day of year
dataPath <- "GOES_data" #Folder where the data is located
savePath <- paste(getwd(),"/",sep="")
print(savePath)
#> [1] "/Users/Kathryn/Documents/PhD_Research/NEFI_pheno/GOESDiurnalNDVI/vignettes/"
siteData <- cbind(siteName,lat,long,TZ)
fileName <- paste(savePath,"GOES_NDVI_Diurnal",siteName,"_",year,day,".csv",sep="")
if(!file.exists(fileName)){ #Note the file of the calculated NDVI's is included on github
calculateNDVI_GOES_MAIN(day=day,siteData=siteData,year=year,TZ=TZ,
dataPath=dataPath,TZ_name="America/New_York",
savePath=savePath)
}fileName <- paste(savePath,"GOES_NDVI_Diurnal",siteName,"_",year,day,".csv",sep="")
GOESdat <- read.csv(fileName,header=FALSE)
plot(as.numeric(GOESdat[3,]),as.numeric(GOESdat[2,]),pch=20,xlim=c(0,20),ylim=c(0,1),
ylab="NDVI",xlab="Time (Hour)")
data <- list(x=as.numeric(GOESdat[3,]),y=as.numeric(GOESdat[2,]))
modelFitFileName <- paste(savePath,siteName,"_",year,day,"_varBurn.RData",sep="")
if(!file.exists(modelFitFileName)){
j.model <- createDiurnalModel(siteName=siteName,data=data)
}if(!file.exists(modelFitFileName)){
var.burn <- runMCMC_Model(j.model=j.model,variableNames=c("a","c","k","prec"),
baseNum=20000,iterSize =10000) #The baseNum and iterSize can be increased/decreased to make the code run faster if you know it will converge easier
save(var.burn,file=modelFitFileName)
}load(modelFitFileName)
plotCI(siteName=siteName,year=year,day=day,savePath=savePath)
out.mat <- data.frame(as.matrix(var.burn))
c <- out.mat$c
c.quantiles <- quantile(c,c(0.025,0.50,0.975))
print(paste("Midday NDVI Quantiles (0.025,0.5,0.975):",as.character(c.quantiles[1]),
as.character(c.quantiles[2]),as.character(c.quantiles[3])))
#> [1] "Midday NDVI Quantiles (0.025,0.5,0.975): 0.627971566237793 0.654016384016397 0.681470456327048"
c.mean <- mean(c)
print(paste("Midday NDVI Mean Estimate:",as.character(c.mean)))
#> [1] "Midday NDVI Mean Estimate: 0.654185756548796"