LUE variations

Evaluate light use efficiency.

Daily data

Gather data and create one nice data frame with daily data.

if (!exists("out_eval_ORG"))  load("./calib_results/out_eval_ORG.Rdata")
if (!exists("out_eval_FULL")) load("./calib_results/out_eval_FULL.Rdata")
load("~/eval_pmodel/data/obs_eval_NT_WITHFORCING.Rdata")
mc <- 12.0107  # molecular mass of carbon

## As done for publication:
ddf <- out_eval_FULL$gpp$fluxnet2015$data$ddf %>% 
  left_join(dplyr::select(obs_eval_NT$ddf, date, sitename, fapar, ppfd_fluxnet2015, temp_day_fluxnet2015, vpd_day_fluxnet2015, prec_fluxnet2015, patm_fluxnet2015), by = c("sitename", "date")) %>% 
  mutate(setup_pmodel = "FULL") %>% 
  
  ## Calculate APAR (using interpolated PPFD and fAPAR for very few points)
  mutate(ppfd_fluxnet2015 = myapprox(ppfd_fluxnet2015),
         fapar            = myapprox(fapar)) %>%
  mutate(apar = ppfd_fluxnet2015 * fapar) %>% 

  ## Remove very low APAR 
  dplyr::filter(apar > 0.1) %>% 
  rename(gpp_obs = obs, gpp_mod = mod) %>% 
  
  ## Calculate LUE (division by molecular mass to get LUE in units of mol/mol!)
  mutate(lue_obs = gpp_obs / (mc * apar),
         lue_mod = gpp_mod / (mc * apar) ) %>% 
  
  ## Remove outliers (outside 1.5 times the inter-quartile-range)
  mutate(lue_obs = remove_outliers(lue_obs),
         lue_mod = remove_outliers(lue_mod))

Look at distribution of values and comparison to obs-derived LUE with daily data.

ddf %>% 
  filter(setup_pmodel == "FULL") %>% 
  drop_na(gpp_obs, gpp_mod) %>% 
  mutate(gpp_mod = gpp_mod * 1e3, gpp_obs = gpp_obs * 1e3) %>% 
  gather(modobs, gpp, c(gpp_obs, gpp_mod)) %>% 
  ggplot(aes(x = modobs, y=gpp, fill=modobs)) +
  geom_violin() +
  theme_classic() +
  scale_fill_brewer(palette="Paired", 
                    breaks=c("gpp_obs", "gpp_mod"), 
                    labels=c("Observed", "Modelled")) +
  labs(fill="", x="Source", y=expression(paste("gpp (mmol mol" ^{-1}, ")"))) + 
  theme(legend.position="top")

ddf %>% 
  filter(setup_pmodel == "FULL") %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  gather(modobs, lue, c(lue_obs, lue_mod)) %>% 
  ggplot(aes(x = modobs, y=lue, fill=modobs)) +
  geom_violin() +
  theme_classic() +
  scale_fill_brewer(palette="Paired", 
                    breaks=c("lue_obs", "lue_mod"), 
                    labels=c("Observed", "Modelled")) +
  labs(fill="", x="Vegetation type", y=expression(paste("LUE (mmol mol" ^{-1}, ")"))) + 
  theme(legend.position="top")

ddf %>% 
  filter(setup_pmodel == "FULL") %>% 
  ggplot() +
  geom_histogram(
    aes(x = apar, y = ..density..), 
    color = "black", alpha = 0.3, binwidth = 0.5, 
    position="identity")

gg <- ddf %>% analyse_modobs2("lue_obs", "lue_mod", type = "heat")

## Loading required package: LSD

## Loading required package: ggthemes

gg$gg

Ok. They are on the same order.

Compare coefficient of variations in LUE model.

cv <- function(vec, ...){
  sd(vec, ...) / mean(vec, ...)
}
df_cv <- ddf %>% 
  summarise(gpp_obs = cv(gpp_obs, na.rm=TRUE),
            ppfd_fluxnet2015 = cv(ppfd_fluxnet2015, na.rm=TRUE),
            fapar = cv(fapar, na.rm=TRUE),
            lue_obs = cv(lue_obs, na.rm=TRUE))
print(df_cv)

## # A tibble: 1 x 4
##   gpp_obs ppfd_fluxnet2015 fapar lue_obs
##     <dbl>            <dbl> <dbl>   <dbl>
## 1    1.08            0.601 0.500   0.883

Monthly

Aggregate to monthly values.

Sum GPP and APAR and get LUE from sums ==> works much better? But is dominated by large values.

mdf_nice <- ddf %>%
  mutate(moy = month(date), year = year(date)) %>%
  group_by(sitename, year, moy) %>%
  summarise( apar = sum(apar)) %>%
  dplyr::filter(apar > 30) %>%
  right_join(out_eval_FULL$gpp$fluxnet2015$data$mdf, by=c("sitename", "year")) %>%
  rename(gpp_obs = obs, gpp_mod = mod) %>%
  mutate(lue_obs = gpp_obs / (mc * apar),
         lue_mod = gpp_mod / (mc * apar) ) %>%
  ungroup()  %>%
  left_join(rsofun::metainfo_Tier1_sites_kgclimate_fluxnet2015 %>% dplyr::select(sitename, classid),
            by="sitename")

## Warning: replacing previous import 'dplyr::intersect' by
## 'lubridate::intersect' when loading 'rsofun'

## Warning: replacing previous import 'dplyr::union' by 'lubridate::union'
## when loading 'rsofun'

## Warning: replacing previous import 'dplyr::setdiff' by 'lubridate::setdiff'
## when loading 'rsofun'

## Warning: replacing previous import 'purrr::list_along' by
## 'rlang::list_along' when loading 'rsofun'

## Warning: replacing previous import 'purrr::invoke' by 'rlang::invoke' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_raw' by
## 'rlang::flatten_raw' when loading 'rsofun'

## Warning: replacing previous import 'purrr::modify' by 'rlang::modify' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::as_function' by
## 'rlang::as_function' when loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_dbl' by
## 'rlang::flatten_dbl' when loading 'rsofun'

## Warning: replacing previous import 'Metrics::ll' by 'rlang::ll' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_lgl' by
## 'rlang::flatten_lgl' when loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_int' by
## 'rlang::flatten_int' when loading 'rsofun'

## Warning: replacing previous import 'purrr::%@%' by 'rlang::%@%' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten_chr' by
## 'rlang::flatten_chr' when loading 'rsofun'

## Warning: replacing previous import 'purrr::splice' by 'rlang::splice' when
## loading 'rsofun'

## Warning: replacing previous import 'purrr::flatten' by 'rlang::flatten'
## when loading 'rsofun'

## Warning: replacing previous import 'purrr::prepend' by 'rlang::prepend'
## when loading 'rsofun'

## Warning: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

## Warning: WARNING:

## Warning: Setting option to specify location of SOFUN (Fortran):

## Warning: rsofun.dir.sofun = '~/sofun/'

## Warning: Change this option by:

## Warning: options( list( rsofun.dir.sofun='string_path_where_sofun_is' ) )

## Warning: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

mdf_nice %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  gather(modobs, lue, c(lue_obs, lue_mod)) %>% 
  ggplot(aes(x = modobs, y=lue, fill=modobs)) +
  geom_violin() +
  theme_classic() +
  scale_fill_brewer(palette="Paired", 
                    breaks=c("lue_obs", "lue_mod"), 
                    labels=c("Observed", "Modelled")) +
  labs(fill="", x="Vegetation type", y=expression(paste("LUE (mmol mol" ^{-1}, ")"))) + 
  theme(legend.position="top")

out_mdf <- mdf_nice %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs", type = "heat")

gg_modobs_mdf_nice <- out_mdf$gg +
  labs(x = expression(paste("Modelled LUE (mmol mol" ^{-1}, ")")),
       y = expression(paste("Observed LUE (mmol mol" ^{-1}, ")")))

gg_modobs_mdf_nice

Mean across daily values ==> works worse? This is now in the paper.

mdf <- ddf %>%
  mutate(moy = month(date), year = year(date)) %>% 
  group_by(sitename, year, moy) %>% 
  summarise(lue_obs = mean(lue_obs, na.rm=TRUE),
            lue_mod = mean(lue_mod, na.rm=TRUE)) %>% 
  mutate(lue_obs = ifelse(is.nan(lue_obs), NA, lue_obs)) %>% 
  mutate(lue_mod = ifelse(is.nan(lue_mod), NA, lue_mod)) %>% 
  ungroup()

mdf %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  gather(modobs, lue, c(lue_obs, lue_mod)) %>% 
  ggplot(aes(x = modobs, y=lue, fill=modobs)) +
  geom_violin() +
  theme_classic() +
  scale_fill_brewer(palette="Paired", 
                    breaks=c("lue_obs", "lue_mod"), 
                    labels=c("Observed", "Modelled")) +
  labs(fill="", x="Vegetation type", y=expression(paste("LUE (mmol mol" ^{-1}, ")"))) + 
  theme(legend.position="top")

out_mdf <- mdf %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs", type = "heat")

gg_modobs_mdf <- out_mdf$gg +
  labs(x = expression(paste("Modelled LUE (mmol mol" ^{-1}, ")")),
       y = expression(paste("Observed LUE (mmol mol" ^{-1}, ")")))
print(gg_modobs_mdf)

Annual

Sum GPP and APAR and get LUE from sums ==> works much better? But is dominated by large values as for monthly

adf_nice <- ddf %>%
  mutate(year = year(date)) %>%
  group_by(sitename, year) %>%
  summarise( apar = sum(apar)) %>%
  # dplyr::filter(apar > 300) %>%
  right_join(out_eval_FULL$gpp$fluxnet2015$data$adf, by=c("sitename", "year")) %>%
  rename(gpp_obs = obs, gpp_mod = mod) %>%
  mutate(lue_obs = gpp_obs / (mc * apar),
         lue_mod = gpp_mod / (mc * apar) ) %>%
  ungroup()  %>%
  left_join(rsofun::metainfo_Tier1_sites_kgclimate_fluxnet2015 %>% dplyr::select(sitename, classid),
            by="sitename")

adf_nice %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  gather(modobs, lue, c(lue_obs, lue_mod)) %>% 
  ggplot(aes(x = modobs, y=lue, fill=modobs)) +
  geom_violin() +
  theme_classic() +
  scale_fill_brewer(palette="Paired", 
                    breaks=c("lue_obs", "lue_mod"), 
                    labels=c("Observed", "Modelled")) +
  labs(fill="", x="Vegetation type", y=expression(paste("LUE (mmol mol" ^{-1}, ")"))) + 
  theme(legend.position="top")

out_adf <- adf_nice %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs")

gg_modobs_adf <- out_adf$gg +
  labs(x = expression(paste("Modelled LUE (mmol mol" ^{-1}, ")")),
       y = expression(paste("Observed LUE (mmol mol" ^{-1}, ")")))

gg_modobs_adf

Mean across daily values ==> works worse? This is now in the paper.

adf <- ddf %>%
  mutate(year = year(date)) %>% 
  group_by(sitename, year) %>% 
  summarise(lue_obs = mean(lue_obs, na.rm=TRUE),
            lue_mod = mean(lue_mod, na.rm=TRUE)) %>% 
  mutate(lue_obs = ifelse(is.nan(lue_obs), NA, lue_obs)) %>% 
  mutate(lue_mod = ifelse(is.nan(lue_mod), NA, lue_mod)) %>% 
  ungroup()

adf %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  gather(modobs, lue, c(lue_obs, lue_mod)) %>% 
  ggplot(aes(x = modobs, y=lue, fill=modobs)) +
  geom_violin() +
  theme_classic() +
  scale_fill_brewer(palette="Paired", 
                    breaks=c("lue_obs", "lue_mod"), 
                    labels=c("Observed", "Modelled")) +
  labs(fill="", x="Vegetation type", y=expression(paste("LUE (mmol mol" ^{-1}, ")"))) + 
  theme(legend.position="top")

out_adf <- adf %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs")

gg_modobs_adf <- out_adf$gg +
  labs(x = expression(paste("Modelled LUE (mmol mol" ^{-1}, ")")),
       y = expression(paste("Observed LUE (mmol mol" ^{-1}, ")")))
print(gg_modobs_adf)

Some sites have really outlier years. These are the ones where modelled LUE is very low and observed very high.

adf <- adf %>% 
  mutate(res = lue_obs - lue_mod)

adf %>% 
  filter(res > quantile(res, 0.99, na.rm = TRUE))

## # A tibble: 9 x 5
##   sitename  year lue_obs lue_mod    res
##   <chr>    <dbl>   <dbl>   <dbl>  <dbl>
## 1 AR-SLu    2009  0.0380 0.0121  0.0259
## 2 AR-SLu    2010  0.0380 0.0129  0.0251
## 3 AR-SLu    2011  0.0345 0.0127  0.0217
## 4 AU-Lox    2008  0.0296 0.0102  0.0194
## 5 BR-Sa3    2004  0.0360 0.0193  0.0167
## 6 CN-Sw2    2010  0.0407 0.0147  0.0260
## 7 DK-Sor    2014  0.0397 0.0223  0.0174
## 8 SD-Dem    2005  0.0518 0.00779 0.0440
## 9 SN-Dhr    2010  0.0454 0.00963 0.0358

adf_nice <- adf_nice %>% 
  mutate(res = lue_obs - lue_mod)

adf_nice %>% 
  filter(res > quantile(res, 0.99, na.rm = TRUE))

## # A tibble: 6 x 11
##   sitename  year  apar gpp_mod     n date       gpp_obs lue_obs lue_mod
##   <chr>    <dbl> <dbl>   <dbl> <int> <date>       <dbl>   <dbl>   <dbl>
## 1 AR-SLu    2011 5674.    863.   365 2011-01-01   2599.  0.0381 0.0127 
## 2 AT-Neu    2007 5493.   1392.   365 2007-01-01   2409.  0.0365 0.0211 
## 3 AT-Neu    2008 5472.   1406.   365 2008-01-01   2621.  0.0399 0.0214 
## 4 AT-Neu    2011 5677.   1411.   365 2011-01-01   2652.  0.0389 0.0207 
## 5 SN-Dhr    2011 3746.    417.   365 2011-01-01   1584.  0.0352 0.00927
## 6 SN-Dhr    2013 3270.    364.   365 2013-01-01   1151   0.0293 0.00927
## # … with 2 more variables: classid <chr>, res <dbl>

Let’s remove them for good and re-calculate correlation of annual values.

adf <- adf %>% 
  filter(res < quantile(res, 0.99, na.rm = TRUE))

out_adf <- adf %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs")

gg_modobs_adf <- out_adf$gg +
  labs(x = expression(paste("Modelled LUE (mmol mol" ^{-1}, ")")),
       y = expression(paste("Observed LUE (mmol mol" ^{-1}, ")")))
gg_modobs_adf

adf_nice <- adf_nice %>% 
  filter(res < quantile(res, 0.99, na.rm = TRUE))


out_adf <- adf_nice %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs")

gg_modobs_adf <- out_adf$gg +
  labs(x = expression(paste("Modelled LUE (mmol mol" ^{-1}, ")")),
       y = expression(paste("Observed LUE (mmol mol" ^{-1}, ")")))
gg_modobs_adf

By site (spatial) and veg type

Now do this with the version that performed better above (adf_nice).

## mean by site
df <- adf_nice %>% 
  group_by(sitename) %>% 
  summarise(lue_obs = mean(lue_obs, na.rm = TRUE),
            lue_mod = mean(lue_mod, na.rm = TRUE)) %>% 
  left_join(rsofun::metainfo_Tier1_sites_kgclimate_fluxnet2015 %>% dplyr::select(sitename, classid), 
            by="sitename") %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3)

## mean by vegetation type
vegdf <- df %>% 
  ungroup() %>% 
  group_by(classid) %>% 
  summarise(lue_obs = mean(lue_obs, na.rm = TRUE),
            lue_mod = mean(lue_mod, na.rm = TRUE))

##################
# FULL setup
##################
out_df <- df %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs", xlab="Modelled LUE", ylab="Observed LUE") # , filnam = "fig/modobs_lue_bysite.pdf"

rsq_lab <- format( out_df$df_metrics %>% dplyr::filter(.metric=="rsq") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 2 )
rmse_lab <- format( out_df$df_metrics %>% dplyr::filter(.metric=="rmse") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )
bias_lab <- format( out_df$df_metrics %>% dplyr::filter(.metric=="bias") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )
slope_lab <- format( out_df$df_metrics %>% dplyr::filter(.metric=="slope") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )
n_lab <- format( out_df$df_metrics %>% dplyr::filter(.metric=="n") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )

out_vegdf <- vegdf %>% 
  rbeni::analyse_modobs2("lue_mod", "lue_obs", xlab="Modelled LUE", ylab="Observed LUE")

rsq_lab2 <- format( out_vegdf$df_metrics %>% dplyr::filter(.metric=="rsq") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 2 )
rmse_lab2 <- format( out_vegdf$df_metrics %>% dplyr::filter(.metric=="rmse") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )
bias_lab2 <- format( out_vegdf$df_metrics %>% dplyr::filter(.metric=="bias") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )
slope_lab2 <- format( out_vegdf$df_metrics %>% dplyr::filter(.metric=="slope") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )
n_lab2 <- format( out_vegdf$df_metrics %>% dplyr::filter(.metric=="n") %>% dplyr::select(.estimate) %>% unlist() %>% unname(), digits = 3 )

## scatterplot
gg_lue_spatial <- ggplot() +
  geom_point(data=df, aes(x=lue_mod, y=lue_obs, color=classid), size=1.5) +
  geom_point(data=vegdf, aes(x=lue_mod, y=lue_obs, color=classid), size=5) +
  geom_smooth(data=df, aes(x=lue_mod, y=lue_obs), method = "lm", color="red", size=0.5) +
  geom_abline(intercept=0, slope=1, linetype="dotted") +
  theme_classic() +
  scale_color_brewer(palette="Paired") +
  labs(
        subtitle = bquote(
          # atop(
            # "Sites:" ~
            italic(R)^2 == .(rsq_lab) ~ (.(rsq_lab2)) ~~~
            # RMSE == .(rmse_lab) ~~~
            # bias == .(bias_lab) ~~~
            slope == .(slope_lab) ~ (.(slope_lab2)) ~~~
            italic(N) == .(n_lab) ~ (.(n_lab2))
          # , 
          # "Vegetation types:" ~
          # italic(R)^2 == .(rsq_lab2) ~~~
          # RMSE == .(rmse_lab2) ~~~
          # bias == .(bias_lab2) ~~~
          # slope == .(slope_lab2) ~~~
          # italic(N) == .(n_lab2) 
        # )
         ),
        x = expression(paste("Modelled LUE (mmol mol" ^{-1}, ")")), 
        y = expression(paste("Observed LUE (mmol mol" ^{-1}, ")")),
        color="")
print(gg_lue_spatial)

# xxxxxx take this
# ggsave("fig/modobs_lue_bysite_classid.pdf", width=5, height=5)


## boxplot 
gg <- df %>% 
  drop_na(lue_obs, lue_mod) %>% 
  mutate(lue_mod = lue_mod * 1e3, lue_obs = lue_obs * 1e3) %>% 
  gather(modobs, lue, c(lue_obs, lue_mod)) %>% 
  ggplot(aes(x=classid, y=lue, fill=modobs)) +
  geom_boxplot() +
  theme_classic() +
  scale_fill_brewer(palette="Paired", 
                    breaks=c("lue_obs", "lue_mod"), 
                    labels=c("Observed", "Modelled")) +
  labs(fill="", x="Vegetation type", y=expression(paste("LUE (mmol mol" ^{-1}, ")"))) + 
  theme(legend.position="top")
print(gg)

# ggsave("fig/boxplot_lue_bysite_classid.pdf", width=5, height=5)

Publication figure

library(cowplot)

## 
## ********************************************************

## Note: As of version 1.0.0, cowplot does not change the

##   default ggplot2 theme anymore. To recover the previous

##   behavior, execute:
##   theme_set(theme_cowplot())

## ********************************************************

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggthemes':
## 
##     theme_map

## The following object is masked from 'package:lubridate':
## 
##     stamp

plot_grid(gg_modobs_mdf, gg_lue_spatial, labels = "auto", label_size = 12, rel_widths = c(1, 1.1))

ggsave("fig/modobs_lue.pdf", width=10, height=5)