In eeholmes/SardineForecast: Indian Oil Sardine Data, Code and References for Papers on Forecasting

knitr::opts_chunk$set(echo = FALSE, warning=FALSE, results='hide', message=FALSE,
                      dev="postscript", dpi=300)

To create an eps file, make sure to knit to html. The eps created when using knit to pdf is very grainy. The file will be in the figures folder. EPS doesn't seem to respect transparency, fyi.

# the getdata() function
source(file.path(here::here(), "inst", "extdata", "get_satellite_data", "Replicate_EMT", "get_EMT_functions.R"))
load(file.path(here::here(), "inst", "extdata", "get_satellite_data", "Replicate_EMT", "indiashp.RData"))
library(ggplot2)
library(tidyr)
library(dplyr)
library(raster)
library(rasterVis)

# Get monthly wind data
#getdata("hawaii_soest_66d3_10d8_0f3c", lat=c(7,15), lon=c(70,78), 
#       altitude=1000, eserver="http://apdrc.soest.hawaii.edu/erddap", alt.name="LEV")
# dfil <- file.path(here::here(), "inst", "extdata", "get_satellite_data", "Replicate_EMT", "hawaii_soest_66d3_10d8_0f3c-7-15-70-78-1979-01-01-2020-06-01.csv")
windmon <- getdata("hawaii_soest_66d3_10d8_0f3c", date=c("1979-01-01", "2020-06-01"), 
        lat=c(7,15), lon=c(70,78), 
        altitude=1000, eserver="http://apdrc.soest.hawaii.edu/erddap", alt.name="LEV")
windmon$time <- as.Date(windmon$time)
windmon$Month <- as.numeric(format(windmon$time, "%m"))
windmon$Year <- as.numeric(format(windmon$time, "%Y"))

windmon.mean <- windmon %>% 
  group_by(Month, latitude, longitude) %>% 
  summarize(v = mean(v, na.rm=TRUE),
            u = mean(u, na.rm=TRUE))
windmon.mean$Season = windmon.mean$Month
windmon.mean$Season[windmon.mean$Month %in% 6:9] <- "Summer Jun-Sep"
windmon.mean$Season[windmon.mean$Month %in% 1:4] <- "Winter Jan-Apr"
windmon.season <- windmon.mean %>% 
  group_by(Season, latitude, longitude) %>% 
  summarize(v = mean(v, na.rm=TRUE),
            u = mean(u, na.rm=TRUE))

plist=list()
for(i in 1:12){
  field <- subset(windmon.mean, Month==i)[, c("longitude", "latitude", "u", "v")]
  field <- rasterFromXYZ(field)
  projection(field) <- CRS("+proj=longlat +datum=WGS84")
  field <- mask(field, indiashp2, inverse=TRUE)
  p1 <- vectorplot(field, isField='dXY', region = FALSE, margin = FALSE, narrows = 100, length=0.04, main=month.name[i], key.arrow = list(size=10, label = 'm/s')) +
    layer(sp.polygons(indiashp2)) + layer(sp.polygons(kerala, fill="grey"))
  plist[[month.name[i]]] <- p1
}
plist2=list()
for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){
  field <- subset(windmon.season, Season==i & longitude>74 & latitude<13 & latitude>7.5)[, c("longitude", "latitude", "u", "v")]
  field <- rasterFromXYZ(field)
  projection(field) <- CRS("+proj=longlat +datum=WGS84")
  field <- mask(field, indiashp2, inverse=TRUE)
  field <- mask(field, kerala, inverse=TRUE)
  p1 <- vectorplot(field, isField='dXY', region = FALSE, margin = FALSE, narrows = 125, length=0.04, main=i, key.arrow = list(size=10, label = 'm/s'), lwd.arrows=1.5) +
    layer(sp.polygons(indiashp2, col="#585858"), under=TRUE) + layer(sp.polygons(kerala, fill="grey", col="#585858"), under=TRUE)
  p1$par.settings$layout.heights$main.key.padding=-4
  p1$main=list(i, x=0.6, just="center")
  plist2[[i]] <- p1
}

gridExtra::grid.arrange(grobs=plist2[1:2], ncol=2)

eeholmes/SardineForecast documentation built on July 17, 2021, 2:56 a.m.