In energyRt/merra2ools: R-tools to process MERRA-2 data (NASA) with focus on renewable energy

options(rmarkdown.html_vignette.check_title = FALSE)
library(knitr)
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.align = "center") 

library(merra2ools)
library(merra2sample)
library(tidyverse)
library(data.table)

Wind speed extrapolation

Hellmann constant

dat <- merra2_mar %>%
  select(-SWGDN,    -ALBEDO) %>% # drop solar radiation
  mutate(h = fH(W10M = W10M, W50M = W50M)) # Hellmann constant

summary(dat$h) # check basic statistics

kableExtra::kbl(
  head(dat), 
  caption = "MERRA-2 subset with estimated Hellmann exponent (h)") %>% 
  kableExtra::kable_styling(position = "center")

Figure code

(h_mean = mean(dat$h))
(h_median = median(dat$h))

# .GeomVline_draw_key <- GeomVline$draw_key
GeomVline$draw_key <- GeomHline$draw_key # adjusting legend

fig <- ggplot(dat) +
  geom_histogram(aes(h, after_stat(density), 
                     fill = "Hellmann", colour = "NA"),
                 alpha = .75, show.legend = F, bins = 50) +
  labs(x = "Hellmann exponent (h)") +
  geom_vline(aes(xintercept = x, linetype = name),
             data = tibble(name = c("mean", "median"),
                           x = c(h_mean, h_median)),
             colour = "blue") +
  scale_linetype(name = "Statistics") +
  theme_bw() + 
  theme(legend.position = c(0.9, .8),
        legend.box.background = element_rect(
          colour = "black", fill = "grey"
          ))

fig

Extrapolation

# extrapolation
for (i in seq(10, 200, 10)) {
  if (i == 10) cat("Height, m: ")
  cat(i)
  dat[[paste0("w", i)]] <- fWSE(i, dat$W10M, dat$h)
  cat(" ")
}

Figure below shows a sample of extrapolated speed of wind from 10 and 50 to up to 200 meters height.

Figure code

dd <- dat %>% ungroup() %>%
  select(UTC, locid, starts_with("w", ignore.case = F)) %>%
  pivot_longer(cols = starts_with("w"), names_prefix = "w") %>%
  mutate(name = as.integer(name))

set.seed(0)
fig <- ggplot(filter(dd, locid %in% sample(locid, 100), UTC == dd$UTC[1])) +
  geom_line(aes(value, name, group = locid, colour = locid), 
            alpha = .75, show.legend = F, size = .75) +
  scale_colour_distiller(palette = "Spectral") +
  geom_hline(yintercept = 10, linetype = "dashed", colour = "blue", alpha = 1) + 
  geom_hline(yintercept = 50, linetype = "dashed", colour = "blue", alpha = 1) +
  # ylim(0, NA) +
  labs(x = "Wind speed, m/s", y = "Height, m") +
  scale_x_continuous(expand = c(0,0)) + 
  scale_y_continuous(expand = c(0,0), limits = c(0, NA), 
                     breaks = c(0, 10, 50, 100, 150, 200),
                     minor_breaks = c(0, 10, seq(25, 200, by = 25)),
                     labels = c(0, 10, 50, 100, 150, 200)) +
  geom_text(data = tibble(name = c("wind speed at 10m from MERRA-2 (W10M)", 
                                   "wind speed at 10m from MERRA-2 (W50M)"), 
                          x = c(19.2, 19.2), y = c(11, 51)),
                          size = 3, colour = "black", alpha = .85,
             aes(x, y, label = name), vjust = 0, hjust = 1) +
             # nudge_x = 0.25, nudge_y = 7) +
  theme_bw()

fig

Wind speed at 50m

locid <- merra2ools::locid %>% select(locid, lon, lat) 
dat2 <- select(dat, UTC, locid, w10, w50, w100, w200) %>%
  full_join(locid)

Figure code

gif_merra(dat2, name = "w50", 1, fps = 12, filename = "merra_wind_50m_12fps.gif", dirname = "tmp")

gif1 <- "images/merra_wind_50m_12fps.gif"
if (file.exists(gif1)) include_graphics(gif1)

Wind power curve

Figure code

# Simplified wind power curve
wpc <- tibble( 
  mps = c(0:25, NA, 25.01, 26:30), # wind speed m/s
  af = fWPC(mps) # availability factor [0,1]
)

fig <- ggplot(wpc) + 
  geom_line(aes(x = mps, y = af, colour = "red"), size = 2, show.legend = F) +
  labs(x = "Wind speed, m/s", 
       y = "Wind turbine availability (capacity) factor") +
  theme_bw()

fig

Estimated hourly availabillity factor

dat2 <- mutate(
  dat2,
  af10 = fWPC(w10),
  af50 = fWPC(w50),
  af100 = fWPC(w100),
  af200 = fWPC(w200)
)

summary(dat2$w50); summary(dat2$af50)
summary(dat2$af50[dat2$w50 > 25])

plot_merra(dat2, "w50", 1, map.border = NA)
plot_merra(dat2, "af50", 1, limits = c(0, 1), 
          palette = "Blues", direction = 1, legend.name = "AF")
plot_merra(dat2, "af10", 1, limits = c(0, 1),
          palette = "Blues", direction = 1, legend.name = "AF")
plot_merra(dat2, "af100", 1, limits = c(0, 1),
          palette = "Blues", direction = 1, legend.name = "AF")
plot_merra(dat2, "af200", 1, limits = c(0, 1), 
          palette = "Blues", direction = 1, legend.name = "AF")

gif_merra(dat2, "af50", 1, limits = c(0, 1), fps = 12, 
         palette = "PuBu", direction = 1,
         legend.name = "AF",
         filename = "merra_wind_af50_12fps.gif")

gif2 <- "images/merra_wind_af50_12fps.gif"
if (file.exists(gif2)) include_graphics(gif2)

# load("data/merra_mean_40y.RData")
d <- merra_mean_40y
summary(d$W50M)

d$w50 <- cut(d$W50M, breaks = wind_levels, include.lowest = T, 
               labels = names(palette.windatlas), right = F)
summary(d$w50)

world <- rnaturalearth::ne_countries(scale = "small", returnclass = "sf")

ggplot(d) +
  geom_tile(aes(lon, lat, fill = w50)) +
  # geom_tile(aes(lon, lat, fill = w50)) +
  # scale_fill_brewer(palette = "Spectral", direction = -1)
  scale_fill_manual(values = palette.windatlas,
                    guide = guide_legend(reverse = TRUE),
                    name = "m/s") + 
  theme_void() + 
  labs(title = "Wind speed at 50 meters height, long-term average (1980-2019)") +
  theme(legend.position = "bottom", legend.spacing.x = unit(1, "mm"),
        legend.margin = margin(0,0,0,0),
        legend.box.margin = margin(0,0,0,0),
        plot.margin = unit(c(0,0,0,0), "mm"),
        plot.title = element_text(hjust = 0.5)) +
  guides(fill = guide_legend(nrow = 1, label.position = "bottom", 
                             keywidth = unit(3, "mm"), keyheight = 5,
                             default.unit = "mm")) +
  scale_x_continuous(expand = c(0, 0), limits = c(-180, 180)) +
  scale_y_continuous(expand = c(0, 0), limits = c(-90, 90)) +
  geom_sf(data = world, color = alpha("black", 1), fill = NA, 
          inherit.aes = F, size = unit(.25, "mm"))

ggsave("images/wind_50m_40y_avr_24steps.png", width = 10.5, height = 6.)

# summary(d$AF100M)

energyRt/merra2ools documentation built on May 2, 2024, 4:53 a.m.