In EarthSystemDiagnostics/orbitalforcing: Orbital Forcing from Insolation

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

orbitalforcing

The goal of orbitalforcing is to ...

library(orbitalforcing)
library(tidyverse)

Annual insolation

tb <- crossing(kyear = 0:1e02,
             lat = seq(0, 90, 10)) %>% 
  group_by(lat) %>% 
  mutate(annual.insolation = AnnualInsolation(kyear, lat))

tb %>% 
  ggplot(aes(kyear, annual.insolation, colour = factor(lat))) +
  geom_line() +
  expand_limits(y = 0)

Daily insolation

tb <- crossing(kyear = 0:5e02,
               lat = c(-45, 0, 45),
               day = c(1, 182)) %>% 
  group_by(lat) %>% 
  mutate(daily.insolation = DailyInsolation(kyear, lat, day)$Fsw,
         season = ifelse(daily.insolation > 200, "Summer", "Winter"),
         season = ifelse(lat == 0, "Equatorial", season))

tb %>% 
  #filter(lat == 50) %>% 
  ggplot(aes(kyear, daily.insolation, colour = factor(lat), group = paste(day, season, lat), linetype = factor(day))) +
  geom_line() +
  facet_grid(season~.)

Max - Min Insolation

tb.range <- tb %>% 
  select(-season) %>% 
  spread(day, daily.insolation) %>% 
  mutate(range.insolation = abs(`1` - `182`))

tb.range %>% 
  ggplot(aes(x = kyear, y = range.insolation,
             colour = factor(lat), group = lat)) +
  geom_line()

PlotWorld <- function(){
  world <- map_data("world")
  worldmap <- ggplot(world, aes(x = long, y = lat, group = group)) +
    geom_polygon(fill = "Grey") +
    coord_quickmap() +
    scale_x_continuous(limits = c(-180, 180), breaks = seq(-180, 180, 60))+
    scale_y_continuous(limits = c(-90, 90), breaks = c(-45, 0, 45)) +
    labs(x = "Longitude", y = "Latitude") +
    theme_bw()
  worldmap
}

PlotWorld() +
  geom_tile(data = glob.grid, aes(x = lon, group = max.amp, fill = max.amp)) +
  scale_fill_viridis_c(option = "inferno") +
  expand_limits(fill = 0)


PlotWorld() +
  geom_tile(data = glob.grid, aes(x = lon, group = min.amp, fill = min.amp))+
  scale_fill_viridis_c(option = "inferno")+
  expand_limits(fill = 0)

PlotWorld() +
  geom_tile(data = glob.grid, aes(x = lon, group = amp.diff, fill = amp.diff))+
  scale_fill_viridis_c(option = "inferno")+
  expand_limits(fill = 0)