The data are from the Fluox, and are already processed for the clear sky days, also the database are filtered for the DOY>=128, the bellow code was apply to be abble of uploading this data to github because the size of the first file was larger than 100Mb:
The time serie consist in observations on Alfafa crop during DOY 144 and 193
#
#namefile <- list.files("data-raw/spectra/", pattern="JB-009-ESA_F")
#
# ori_tab <- readr::read_csv(paste0("data-raw/spectra/",namefile),col_names = TRUE)
# col_Names <- as.character(ori_tab$F_final1)
#
# trans_tab <- t(ori_tab)
# colnames(trans_tab) <- col_Names
# trans_tab <- as.data.frame(trans_tab)
#
# trans_tab <- trans_tab[-1, ]
# trans_tab <- trans_tab[,-2]
# colnames(trans_tab)[1] <- "DOY"
#
# trans_tab <- trans_tab |>
# dplyr::filter(
# DOY > 127.9
# )
# namefile_rds <- sub(".csv",".rds",namefile)
# readr::write_rds(trans_tab, paste0("data/",namefile_rds))Now we are keeping going from here with the .rds file
namefile <- list.files("data-raw/", pattern = "_F_")
ori_tab <- readr::read_rds(paste0("data-raw/",namefile))
names(ori_tab) <- c("DOY",
round(
as.numeric(names(ori_tab)[2:length(ori_tab)]), 2)
)Reading the Index file
ind_file <- readr::read_csv("data-raw/table_all_index.csv") |>
dplyr::filter(DOYdayfrac > 127.9)First we are gonna make the time table and made a mean for every 30 m
time_table <- cbind(ind_file,ori_tab)
time_table <- time_table |>
dplyr::filter(DOY >= 144 & DOY <=193)
time_table <- time_table |>
dplyr::filter(lubridate::hour(time_table$UTC_datetime) > 8)
min_int<-"30 min" #step
sttime <- time_table$UTC_datetime[1] # start time of the serie
endtime <- time_table$UTC_datetime[8344] #final time of the serie
## Making the agreggation
timetoagg<-seq.POSIXt(from = as.POSIXct(sttime,tz = "UTC"),
to = as.POSIXct(endtime,tz = "UTC"),
by = min_int)
time_agg <- aggregate(x = time_table,
by = list(cut(as.POSIXct(time_table$UTC_datetime,tz = "UTC"),
timetoagg)), FUN = mean)Now we need to load the GPP data and find the days and hours that bolth have measurments
gppFile <- list.files("data-raw/", pattern = "aa_")
gpp_table <- readr::read_csv(paste0("data-raw/", gppFile))If we try to use the function dplyr::semi_join() we are not gonna get
any result, because the data frames don’t have any column to compare, so
we must create in this case
time_agg <- time_agg |> dplyr::mutate(day = time_agg$DOY%/%1,
Hr1 = lubridate::hour(as.character(Group.1)),
Mn1 = lubridate::minute(as.character(Group.1)))
time_agg <- time_agg[,-7] # this line is because they have the same same in the GPP table, but with no match
gpp_table <- gpp_table |> dplyr::mutate(day=DOY%/%1)
# Now we can use the semi_join
comp_f_gpp <- dplyr::semi_join(time_agg,gpp_table)
comp_gpp_f <- dplyr::semi_join(gpp_table,time_agg)
comp_gpp_f <- comp_gpp_f[,-c(6,176)]
df_final <- cbind(comp_f_gpp,comp_gpp_f)Now we are gonna do a integral in the F spectra, but first we have to do
some transformations in the time_agg.
spectra_Table <- time_agg[,-c(689:691)]
spectra_Table <- as.data.frame(t(spectra_Table))
c_spectra <- paste0(lubridate::hour(as.character(spectra_Table[1,])),":",
lubridate::minute(as.character(spectra_Table[1, ])))
colnames(spectra_Table) <- c_spectra
spectra_Table <- spectra_Table[-c(1:6), ]
rownames(spectra_Table) <- seq(1:682)
spectra <- names(ori_tab[0, c(2:683)])
spectra_Table <- cbind(spectra,spectra_Table)
colnames(spectra_Table) <- make.unique(names(spectra_Table))
spectra_Table[,1] <-as.numeric(spectra_Table[,1])Integrating function:
integral <- function(x,y){
approxfun(x,y)
}
result <- vector("list", length = ncol(spectra_Table))
for (i in 1:(ncol(spectra_Table))){
int <- integral(spectra_Table$spectra, spectra_Table[,i])
result[[i]] <- integrate(int, min(spectra_Table$spectra),
max(spectra_Table$spectra))
}
for(i in 1:length(result)){
if(i==1){
fint <- result[[i]]$value
}else{
fint_a <- result[[i]]$value
fint <- rbind(fint,fint_a)
fint <- as.data.frame(fint)
}
}
row.names(fint) <- seq(1:208)
fint <- as.data.frame(fint[-1,]) # the first observasation is the integral of wavelenght in function of it self
colnames(fint) <- "Fint"
## Now we are gonna add the integral in the final dataset
df_final <- cbind(df_final,fint)Before going plot the data, first we need another transformations to do a correlation plot for all the Fluoresce spectra
fint_t <- vector("numeric", length = 208)
fint_t[2:208] <- t(fint)
spectra_Table[683, ] <- fint_t
gpp <- vector("numeric", length = 208)
gpp[2:208] <- comp_gpp_f$GPP_DT_U95
spectra_Table[684, ] <- gpp
for(i in 2:ncol(spectra_Table)){
if(i==2){
tab <- as.numeric(spectra_Table[,i])
}else{
tab_a <- as.numeric(spectra_Table[,i])
tab <- cbind(tab,tab_a)
}
}
tab <- as.data.frame(tab)
colnames(tab) <- comp_gpp_f$TIMESTAMP
tab <- as.data.frame(t(tab))
colnames(tab) <- spectra_Table$spectra
colnames(tab)[683] <- "fint"
colnames(tab)[684] <- "gpp"
##
pear_resul <- vector("list", length = ncol(tab))
for(i in 1:length(pear_resul)){
pear_resul[[i]] <- cor.test(tab$gpp, tab[,i])
}
for(i in 1:length(pear_resul)){
if(i==1){
correl <- pear_resul[[i]]$estimate
}else{
correl_a <- pear_resul[[i]]$estimate
correl <- rbind(correl,correl_a)
correl <- as.data.frame(correl)
}
}
rownames(correl) <- colnames(tab)
correl[,2] <- spectra_Table$spectraPloting correlation in the Fluoresce spectra, without F integral.
F integral r = 0.65
correl[-c(683:684),] |>
ggplot2::ggplot(ggplot2::aes(x=as.numeric(V2), y = cor))+
ggplot2::geom_point()+
ggplot2::xlab(label="Fluorescence Spectra in nm")
Now we can subset the data frame for the most importants wavelength,
685nm,705nm, 740nm, 780nm, 757nm, 771nm and F integral,
and with the GPP data include for the Day Time (DT) observations
df_final_sub <- df_final |>
dplyr::select(
Group.1, day, Hr1, L750_Wm2nm1sr1, SZA, PAR_Wm2,
`685.09`,`705.1`,`779.96` ,`740.02`,`757.11`, `771.06`, Fint,
GPP_DT, GPP_DT_U95
) |>
dplyr::glimpse()## Rows: 207
## Columns: 15
## $ Group.1 <fct> 2018-05-25 09:00:02, 2018-05-25 09:30:02, 2018-05-25 10~
## $ day <dbl> 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, ~
## $ Hr1 <int> 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 1~
## $ L750_Wm2nm1sr1 <dbl> 0.29679919, 0.31430288, 0.32663515, 0.33453771, 0.33726~
## $ SZA <dbl> 32.92583, 28.54777, 24.98358, 22.66306, 21.92419, 22.93~
## $ PAR_Wm2 <dbl> 369.93567, 391.71895, 406.27653, 415.26607, 417.01186, ~
## $ `685.09` <dbl> 1.2308437, 1.3439059, 1.4053925, 1.4662777, 1.4739017, ~
## $ `705.1` <dbl> 1.5574706, 1.6609758, 1.7503835, 1.8034656, 1.8040760, ~
## $ `779.96` <dbl> 1.5097738, 1.5652948, 1.7294112, 1.7645806, 1.8263655, ~
## $ `740.02` <dbl> 4.334771, 4.471682, 4.948808, 5.048612, 5.246750, 5.469~
## $ `757.11` <dbl> 3.0777047, 3.1789153, 3.5247521, 3.5942228, 3.7366047, ~
## $ `771.06` <dbl> 1.9703019, 2.0396786, 2.2565206, 2.3006356, 2.3854447, ~
## $ Fint <dbl> 255.14553, 266.40490, 290.57136, 297.53778, 305.94678, ~
## $ GPP_DT <dbl> 24.2, 25.1, 25.9, 26.4, 26.7, 24.8, 24.7, 25.4, 26.1, 2~
## $ GPP_DT_U95 <dbl> 25.1, 26.2, 27.1, 27.7, 28.1, 26.6, 26.6, 27.0, 27.4, 2~
df_final_sub[,-c(1:3)] |> #excluding the time for the Pearson correlation
cor(use = "p") |>
corrplot::corrplot()
time_agg |>
tidyr::pivot_longer(
cols = "669.98":"779.96",
names_to = "wavelength",
values_to = "Fluorescence") |>
ggplot2::ggplot(ggplot2::aes(x=as.numeric(wavelength), y=Fluorescence, color=DOYdayfrac)) +
ggplot2::geom_point() +
ggplot2::geom_line()+
ggplot2::facet_wrap(~day)+
ggplot2::xlab(label="Wavelength") 
df_final_sub |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= Fint))+
ggplot2::geom_jitter(ggplot2::aes(colour=as.factor(day)))+
ggplot2::geom_smooth(method = "lm")+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~")))
Controlling by the hours between 11:00 AM and 01:00 PM
df_final_sub |>
dplyr::filter(Hr1 > 10 & Hr1 < 14) |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= Fint))+
ggplot2::geom_jitter(ggplot2::aes(colour=as.factor(day)))+
ggplot2::geom_smooth(method = "lm")+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~")))
With the 780 wavelength
df_final_sub |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= `779.96`))+
ggplot2::geom_jitter(ggplot2::aes(colour=as.factor(day)))+
ggplot2::geom_smooth(method = "lm")+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~")))
df_final_sub |>
dplyr::filter(Hr1 > 10 & Hr1 < 14) |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= `779.96`))+
ggplot2::geom_jitter(ggplot2::aes(colour=as.factor(day)))+
ggplot2::geom_smooth(method = "lm")+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~")))
With the other two selected waves
df_final_sub |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= `740.02`))+
ggplot2::geom_point()+
ggplot2::geom_smooth(method = "lm")+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~")))
df_final_sub |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= `685.09`))+
ggplot2::geom_point()+
ggplot2::geom_smooth(method = "lm")+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~")))
df_final_sub |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= Fint))+
ggplot2::geom_point()+
ggplot2::geom_smooth(method = "lm")+
ggplot2::facet_wrap(~day, scales = "free")+
ggplot2::theme_classic()+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..rr.label..)))
df_final_sub |>
ggplot2::ggplot(ggplot2::aes(x=GPP_DT_U95, y= `779.96`))+
ggplot2::geom_point()+
ggplot2::geom_smooth(method = "lm")+
ggplot2::facet_wrap(~day, scales = "free")+
ggplot2::theme_classic()+
ggpubr::stat_regline_equation(ggplot2::aes(
label = paste(..rr.label..)))
Detection of clounds using the incoming radiance at 750 nm
for(i in c(145,146,161,171,177,181,182,183,187,188,189,190,191)){
print(time_agg |>
dplyr::filter(DOYdayfrac%/%1 == i) |>
ggplot2::ggplot(ggplot2::aes(x=DOYdayfrac, y= L750_Wm2nm1sr1))+
ggplot2::geom_jitter()+
ggplot2::geom_smooth())
}
Here we are gonna see how the reflectance of the canopy are changing, and for this we are gonna use the NDVI and NDVI red-edge
rfile <- list.files("data-raw/", pattern = "_R_")
rtab <- readr::read_rds(paste0("data-raw/",rfile))
rspec <- as.data.frame(t(rtab))
c_name <- rspec[1,]
colnames(rspec) <- c_name
rspec <- rspec[-1, ]Using only the bands that are necessary
rnir <- rspec[c(616:682),]
rnir_spec <- as.numeric(rownames(rnir))
rnir <- cbind(rnir_spec,rnir)
rred <- rspec[c(1:49),]
rred_spec <- as.numeric(rownames(rred))
rred <- cbind(rred_spec,rred)
redge <- rspec[c(207:267),]
redge_spec <- as.numeric(rownames(redge))
redge <- cbind(redge_spec,redge)doing the integral for these bands
rnir_result <- vector("list", length = ncol(rnir))
rred_result <- vector("list", length = ncol(rred))
redge_result <- vector("list", length = ncol(redge))
for (i in 1:(ncol(rnir))){
int <- integral(rnir$rnir_spec, rnir[,i])
rnir_result[[i]] <- integrate(int, min(rnir$rnir_spec), max(rnir$rnir_spec))
}
for(i in 1:length(rnir_result)){
if(i==1){
rnir_int <- rnir_result[[i]]$value
}else{
rnir_int_a <- rnir_result[[i]]$value
rnir_int <- rbind(rnir_int,rnir_int_a)
rnir_int <- as.data.frame(rnir_int)
}
}
rnir_int <- rnir_int[-1,]
rnir_int <- as.data.frame(rnir_int)
for (i in 1:(ncol(rred))){
int <- integral(rred$rred_spec, rred[,i])
rred_result[[i]] <- integrate(int, min(rred$rred_spec), max(rred$rred_spec))
}
for(i in 1:length(rred_result)){
if(i==1){
rred_int <- rred_result[[i]]$value
}else{
rred_int_a <- rred_result[[i]]$value
rred_int <- rbind(rred_int,rred_int_a)
rred_int <- as.data.frame(rred_int)
}
}
rred_int <- rred_int[-1,]
rred_int <- as.data.frame(rred_int)
for (i in 1:(ncol(redge))){
int <- integral(redge$redge_spec, redge[,i])
redge_result[[i]] <- integrate(int, min(redge$redge_spec), max(redge$redge_spec))
}
for(i in 1:length(redge_result)){
if(i==1){
redge_int <- redge_result[[i]]$value
}else{
redge_int_a <- redge_result[[i]]$value
redge_int <- rbind(redge_int,redge_int_a)
redge_int <- as.data.frame(redge_int)
}
}
redge_int <- redge_int[-1,]
redge_int <- as.data.frame(redge_int)Now the formula for the both ndvi
ndvi <- (rnir_int - rred_int)/(rnir_int+rred_int)
colnames(ndvi) <- "NDVI"
ndvi_r <- (rnir_int- redge_int)/(rnir_int + redge_int)
colnames(ndvi_r) <- "NDVI Red"
rtab <- cbind(rtab,ndvi,ndvi_r)
mean_ndvi <- rtab |>
dplyr::mutate(day = DOY%/%1) |>
dplyr::group_by(day) |>
dplyr::summarise(ndvi = mean(NDVI))
mean_ndvi_r <- rtab |>
dplyr::mutate(day = DOY%/%1) |>
dplyr::group_by(day) |>
dplyr::summarise(ndvi_r = mean(`NDVI Red`))We also are gonna need the NDVI aggregated, so now we are gonna do the same process with this data
rtab <- cbind(ind_file,rtab)
rtab <- rtab |>
dplyr::filter(lubridate::hour(rtab$UTC_datetime) > 8)
ndvi_agg <- aggregate(x = rtab,
by = list(cut(as.POSIXct(rtab$UTC_datetime,tz = "UTC"),
timetoagg)), FUN = mean)
ndvi_agg <- ndvi_agg |> dplyr::mutate(day = DOY%/%1,
Hr1 = lubridate::hour(as.character(Group.1)),
Mn1 = lubridate::minute(as.character(Group.1)))
ndvi_agg <- ndvi_agg[,-7]
comp_ndvi_gpp <- dplyr::semi_join(ndvi_agg,gpp_table)## Joining, by = c("day", "Hr1", "Mn1")
Calculatin the NIRv for latter calculate the aPAR
t_rnir <- data.frame(t(rnir))
names(t_rnir) <- round(as.numeric(t_rnir[1,]),2)
t_rnir <- t_rnir[-1,]
t_rnir <- cbind(ind_file,t_rnir)
t_rnir <- t_rnir |>
dplyr::filter(lubridate::hour(UTC_datetime) > 8 & DOYdayfrac >= 144 & DOYdayfrac <=193)
min_int<-"30 min"
sttime <- t_rnir$UTC_datetime[1]
endtime <- t_rnir$UTC_datetime[8344]
timetoagg<-seq.POSIXt(from = as.POSIXct(sttime,tz = "UTC"),
to = as.POSIXct(endtime,tz = "UTC"),
by = min_int)
nir_agg <- aggregate(x = t_rnir,
by = list(cut(as.POSIXct(t_rnir$UTC_datetime,tz = "UTC"),
timetoagg)), FUN = mean)
nir_agg <- nir_agg |>
dplyr::select(
Group.1,
`775`
)
names(nir_agg)[2] <- "R775"
nir_agg <- dplyr::semi_join(nir_agg,comp_ndvi_gpp)## Joining, by = "Group.1"
And Finaly the plots to see the changes
mean_ndvi |>
ggplot2::ggplot(ggplot2::aes(x=day, y=ndvi))+
ggplot2::geom_jitter()+
ggplot2::geom_point(data = subset(mean_ndvi, ndvi< 0.9),
ggplot2::aes(x=day,y=ndvi)
, colour="red")
mean_ndvi_r |>
ggplot2::ggplot(ggplot2::aes(x=day, y=ndvi_r))+
ggplot2::geom_jitter()+
ggplot2::geom_point(data = subset(mean_ndvi_r, ndvi_r< 0.59),
ggplot2::aes(x=day,y=ndvi_r)
, colour="red")
the following figure show to us, the days in the data-set that have problems and the most provable cause, this days we are gonna exclud for our analyses.
NDVI <- comp_ndvi_gpp$NDVI
NDVI_R <- comp_ndvi_gpp$`NDVI Red`
NIR = NDVI*nir_agg$R775
ndvi_df <- data.frame(NDVI,NDVI_R, NIR)
df_final_sub <- cbind(df_final_sub,ndvi_df)
df_tab <- df_final_sub |>
dplyr::filter(day != 171 & day !=177 & day !=188 & day != 189
& day !=199 & day != 200 & day != 201)The calculos for the aPAR is based on the Cesana, I., 2021 thesis
(not available yet)
## Rows: 143
## Columns: 23
## $ Group.1 <fct> 2018-05-25 09:00:02, 2018-05-25 09:30:02, 2018-05-25 10~
## $ day <dbl> 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, 145, ~
## $ Hr1 <int> 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 1~
## $ L750_Wm2nm1sr1 <dbl> 0.29679919, 0.31430288, 0.32663515, 0.33453771, 0.33726~
## $ SZA <dbl> 32.92583, 28.54777, 24.98358, 22.66306, 21.92419, 22.93~
## $ PAR_Wm2 <dbl> 369.93567, 391.71895, 406.27653, 415.26607, 417.01186, ~
## $ `685.09` <dbl> 1.2308437, 1.3439059, 1.4053925, 1.4662777, 1.4739017, ~
## $ `705.1` <dbl> 1.5574706, 1.6609758, 1.7503835, 1.8034656, 1.8040760, ~
## $ `779.96` <dbl> 1.5097738, 1.5652948, 1.7294112, 1.7645806, 1.8263655, ~
## $ `740.02` <dbl> 4.334771, 4.471682, 4.948808, 5.048612, 5.246750, 5.469~
## $ `757.11` <dbl> 3.0777047, 3.1789153, 3.5247521, 3.5942228, 3.7366047, ~
## $ `771.06` <dbl> 1.9703019, 2.0396786, 2.2565206, 2.3006356, 2.3854447, ~
## $ Fint <dbl> 255.14553, 266.40490, 290.57136, 297.53778, 305.94678, ~
## $ GPP_DT <dbl> 24.2, 25.1, 25.9, 26.4, 26.7, 24.8, 24.7, 25.4, 26.1, 2~
## $ GPP_DT_U95 <dbl> 25.1, 26.2, 27.1, 27.7, 28.1, 26.6, 26.6, 27.0, 27.4, 2~
## $ NDVI <dbl> 0.9571886, 0.9563613, 0.9564403, 0.9557606, 0.9551864, ~
## $ NDVI_R <dbl> 0.6397319, 0.6398108, 0.6416534, 0.6404618, 0.6407936, ~
## $ NIR <dbl> 0.5368435, 0.5222212, 0.5051232, 0.5054422, 0.5082543, ~
## $ k1 <dbl> 229.4349, 191.0885, 158.6574, 137.0217, 130.0536, 139.5~
## $ k2 <dbl> 2.581094, 2.535827, 2.503354, 2.484389, 2.478717, 2.486~
## $ aPAR <dbl> 225.5134, 185.3108, 151.2731, 130.2763, 123.8156, 133.1~
## $ Fy <dbl> 1.1313985, 1.4376113, 1.9208402, 2.2838981, 2.4709875, ~
## $ LUE <dbl> 0.11130159, 0.14138410, 0.17914625, 0.21262503, 0.22695~
Saving the data in a .rds file
readr::write_rds(df_tab,"data/base.rds")