Gets global solar irradiance data from the NASA POWER agroclimatology data set, then plots it. This gist compares daily light integral near London, England, in the first 149 days of 2018 to the first 149 days of 2019.

nasa_power_wisley.R

# to look at The Wisley location data 2018 and 2019 based on 
# this Q from Tom Coulson https://twitter.com/TomCoulson85/status/1135905779242995712

library(nasapower)
library(ggplot2)
library(cowplot)
library(RColorBrewer)
library(zoo)

# get the 2018 data
d18 <- get_power(
  community = "AG",
  lonlat = c(-0.48, 51.3),
  pars = c('ALLSKY_SFC_SW_DWN'),
  dates = c("2018-01-01", "2018-05-31"),
  temporal_average = "DAILY"
)

# get the 2019 data
d19 <- get_power(
  community = "AG",
  lonlat = c(-0.48, 51.3),
  pars = c('ALLSKY_SFC_SW_DWN'),
  dates = c("2019-01-01", "2019-05-31"),
  temporal_average = "DAILY"
)

# 6 February 2018 has missing data. I use the na.approx function to approximate it
d18$ALLSKY_SFC_SW_DWN <- ifelse(d18$ALLSKY_SFC_SW_DWN < 0, NA, d18$ALLSKY_SFC_SW_DWN)
d18$ALLSKY_SFC_SW_DWN <- na.approx(d18$ALLSKY_SFC_SW_DWN)

# put the 2018 and 2019 data together for easy plotting
d <- rbind.data.frame(d18, d19)

# make a vector with all dates on the same year, so the x-axis labels
# can easily show as month, rather than as day of year
d$fakeDate <- as.Date('2018-12-31') + d$DOY

# data for the last 2 days of May 2019 aren't in the dataset yet, so omit data from May 30 and May 31
d <- subset(d, DOY < 150)

# show the MJ/m2 global solar irradiance in photosynthetic units, using conversion of Meek et al 1984
d$dli <- d$ALLSKY_SFC_SW_DWN * 2.04

# make a plot to show the DLI from Jan 1 to May 29 in 2018 and 2019
p <- ggplot(data = d, aes(x = fakeDate, y = dli))
yo <- p + theme_cowplot(font_family = 'Fira Sans') +
  background_grid(major = 'xy') +
  geom_smooth(se = FALSE, aes(colour = as.factor(YEAR))) +
  geom_point(shape = 1, aes(colour = as.factor(YEAR))) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title = element_blank(),
        plot.title = element_text(hjust = 0),
        plot.caption = element_text(size = 9),
        legend.position = c(0.2, 0.8)) +
  labs(x = 'Date',
       y = expression(paste('Daily light integral (mol ', m^{-2}, d^{-1}, ')')),
       caption = 'calculated from NASA POWER satellite data, https://power.larc.nasa.gov/\nfor 1/2 degree by 1/2 degree resolution region containing The Wisley, Woking, UK',
       title = 'Daily light integral (DLI) in 2018 & 2019')

yo

save_plot('~/Desktop/dli_by_day.png', yo, base_aspect_ratio = 1.78)

# calulate the cumulative sum of DLI for each year
d$runsum <- ave(d$dli, d$YEAR, FUN = cumsum)

# plot the cumulative sum
p <- ggplot(data = d, aes(x = fakeDate, y = runsum))
yo2 <- p + theme_cowplot(font_family = 'Fira Sans') +
  background_grid(major = 'xy') +
  geom_step(aes(colour = as.factor(YEAR))) +
  scale_colour_brewer(palette = 'Set1') +
  theme(legend.title = element_blank(),
        plot.title = element_text(hjust = 0),
        plot.caption = element_text(size = 9),
        legend.position = c(0.8, 0.2)) +
  labs(x = 'Date',
       y = expression(paste('Cumulative sum of DLI (mol ', m^{-2}, ')')),
       caption = 'calculated from NASA POWER satellite data, https://power.larc.nasa.gov/\nfor 1/2 degree by 1/2 degree resolution region containing The Wisley, Woking, UK',
       title = 'Daily light integral (DLI) through 29 May in 2018 & 2019')

yo2

save_plot('~/Desktop/dli_cumsum.png', yo2, base_aspect_ratio = 1.78)