In eeholmes/SardineForecast: Indian Oil Sardine Data, Code and References for Papers on Forecasting
# to make the 300dpi for submission
# delete the eval=FALSE
# knit to html
# look in Figure_files/figure_html for the eps files
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)
library(gridExtra)
library(grid)
library(marmap)
Get the data
# Get monthly sst data ERA5
#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")
sstmon <- getdata("hawaii_soest_d124_2bb9_c935", date=c("1979-01-01", "2020-06-01"),
lat=c(7,15), lon=c(70,78), pars="sst",
eserver="http://apdrc.soest.hawaii.edu/erddap")
sstmon$time <- as.Date(sstmon$time)
sstmon$Month <- as.numeric(format(sstmon$time, "%m"))
sstmon$Year <- as.numeric(format(sstmon$time, "%Y"))
sstmon$sst <- sstmon$sst - 273.15
# Get monthly sst data ERA5
#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 <- getEMT(windmon, coast_angle=158)
windmon$EMTperpx <- cos(pi*(180-158)/180)*windmon$EMTperp
windmon$EMTperpy <- sin(pi*(180-158)/180)*windmon$EMTperp
# Get monthly ocean currents data
currentsmon <- getdata("erdTAgeomday", lat=c(7,15), lon=c(70,78), altitude=0)
currentsmon$time <- as.Date(currentsmon$time)
currentsmon$Month <- as.numeric(format(currentsmon$time, "%m"))
currentsmon$Year <- as.numeric(format(currentsmon$time, "%Y"))
currentsmon.mean <- currentsmon %>%
group_by(Month, latitude, longitude) %>%
summarize(u = mean(u_current, na.rm=TRUE),
v = mean(v_current, na.rm=TRUE))
currentsmon.mean$Season = currentsmon.mean$Month
currentsmon.mean$Season[currentsmon.mean$Month %in% 6:9] <- "Summer Jun-Sep"
currentsmon.mean$Season[currentsmon.mean$Month %in% 1:4] <- "Winter Jan-Apr"
currentsmon.mean$Season[currentsmon.mean$Month %in% 10:12] <- "Fall Oct-Dec"
currentsmon.season <- currentsmon.mean %>%
group_by(Season, latitude, longitude) %>%
summarize(u = mean(u, na.rm=TRUE),
v = mean(v, na.rm=TRUE))
Make the mean dataframes
windmon.mean <- windmon %>%
group_by(Month, latitude, longitude) %>%
summarize(v = mean(v, na.rm=TRUE),
u = mean(u, na.rm=TRUE),
EMT = mean(EMTperp, na.rm=TRUE),
We = mean(We, na.rm=TRUE),
EMTx = mean(EMTperpx, na.rm=TRUE),
EMTy = mean(EMTperpy, 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),
EMT = mean(EMT, na.rm=TRUE),
We = mean(We, na.rm=TRUE),
EMTx = mean(EMTx, na.rm=TRUE),
EMTy = mean(EMTy, na.rm=TRUE))
windmon$Season = windmon$Month
windmon$Season[windmon$Month %in% 6:9] <- "Summer Jun-Sep"
windmon$Season[windmon$Month %in% 1:4] <- "Winter Jan-Apr"
windmon.season.2014 <- windmon %>% subset(Year==2014) %>%
group_by(Season, latitude, longitude) %>%
summarize(v = mean(v, na.rm=TRUE),
u = mean(u, na.rm=TRUE),
We = mean(We, na.rm=TRUE),
EMT = mean(EMTperp, na.rm=TRUE))
sstmon.mean <- sstmon %>%
group_by(Month, latitude, longitude) %>%
summarize(sst = mean(sst, na.rm=TRUE))
sstmon.mean$Season = sstmon.mean$Month
sstmon.mean$Season[sstmon.mean$Month %in% 6:9] <- "Summer Jun-Sep"
sstmon.mean$Season[sstmon.mean$Month %in% 1:4] <- "Winter Jan-Apr"
sstmon.season <- sstmon.mean %>%
group_by(Season, latitude, longitude) %>%
summarize(sst = mean(sst, na.rm=TRUE))
Monthly EMT
for(i in 1:12){
rdf <- subset(windmon.mean, Month==i)
rdf$x <- rdf$longitude
rdf$y <- rdf$latitude
rdf <- rdf[, c("x", "y", "We")]
if(i==1){ rst <- rasterFromXYZ(rdf); names(rst) <- month.name[i] }else{
tmp <- rasterFromXYZ(rdf); names(tmp) <- month.name[i]
rst <- stack(rst, tmp)
}
}
proj4string(rst) <- "+proj=longlat +datum=WGS84"
at=seq(-3000,0,50)
pal <- colorRampPalette(c("white","red", "green", "blue"))
spplot(rst, col.regions=rev(pal(length(at))), at=at, zlim=c(-3000,0), scales=list(draw = TRUE),
main="Average Monthly EMT perp to coast (+ in, - out)") + layer(sp.polygons(indiashp2, fill="white"))
Seasonal EMT
for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){
rdf <- subset(windmon.season.2014, Season==i)
rdf$x <- rdf$longitude
rdf$y <- rdf$latitude
rdf <- rdf[, c("x", "y", "EMT")]
if(i=="Summer Jun-Sep"){ rst <- rasterFromXYZ(rdf); names(rst) <- i }else{
tmp <- rasterFromXYZ(rdf); names(tmp) <- i
rst <- stack(rst, tmp)
}
}
at=seq(-3000,0,50)
pal <- colorRampPalette(c("white","red", "green", "blue", "black"))
spplot(rst, col.regions=rev(pal(length(at))), at=at, zlim=c(-3000,0), axes=TRUE,
main="Average Monthly EMT perp to coast (+ in, - out)") + layer(sp.polygons(indiashp2))
Wind with EMT overlaid
plist2 <- list()
windmon.season$EMT[windmon.season$EMT< -4000] <- -4000
at=seq(-4000,50,50)
pal <- colorRampPalette(c("white", "green", "blue", "black", "red"))
for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){
field <- subset(windmon.season, Season==i)[, c("longitude", "latitude", "u", "v")]
field <- rasterFromXYZ(field)
projection(field) <- CRS("+proj=longlat +datum=WGS84")
rdf <- subset(windmon.season, Season==i)
rdf$x <- rdf$longitude
rdf$y <- rdf$latitude
rdf <- rdf[, c("x", "y", "EMT")]
rst <- rasterFromXYZ(rdf)
projection(rst) <- CRS("+proj=longlat +datum=WGS84")
rst <- mask(rst, indiashp2, inverse=TRUE)
field <- mask(field, indiashp2, inverse=TRUE)
field <- mask(field, kerala, inverse=TRUE)
p1 <- vectorplot(field, isField='dXY', region = rst, margin = FALSE, narrows = 75,
length=0.04, main=i, key.arrow = list(size=10, label = 'm/s'), lwd.arrows=1.5,
col.regions=rev(pal(length(at))), at=at,
zlim=c(-3000,50), zlimcol="black", colorkey=FALSE, xlab="") +
layer(sp.polygons(indiashp2, col="#585858")) + 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
}
key <- spplot(rst, scales = list(draw = TRUE),
colorkey = list(space = "bottom"), at=at, zlim=c(-3000,50), col.regions=rev(pal(length(at))))
key <- draw.colorkey(key$legend[[1]]$args$key)
plist2[[3]] <- key
key.label <- grid::textGrob(expression(kg~m^-1~s^-1), x = unit(0.02, "npc"), y = unit(0.6, "npc"), just = "left")
plist2[[4]] <- key.label
lay <- rbind(c(1,1,2,2),
c(NA,3,3,4))
gridExtra::grid.arrange(grobs=plist2[1:4], layout_matrix = lay, widths = c(1,1,1,1),
heights=c(10,2))
Wind with SST overlaid
library(marmap)
bat <- getNOAA.bathy(lon1 = 70, lon2 = 78,
lat1 = 7, lat2 = 15, resolution = 4)
save(kerala, indiashp, indiashp2, bat, file=file.path(here::here(), "inst", "extdata", "get_satellite_data", "Replicate_EMT", "indiashp.RData"))
plist2 <- list()
at=seq(26,31,0.1)
zlims=c(26,31)
pal <- colorRampPalette(c("red", "white", "blue"))
bat.200m <- rasterToContour(as.raster(bat), levels= -200)
for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){
field <- subset(windmon.season, Season==i)[, c("longitude", "latitude", "u", "v")]
field <- rasterFromXYZ(field)
projection(field) <- CRS("+proj=longlat +datum=WGS84")
rdf <- subset(sstmon.season, Season==i)
rdf$x <- rdf$longitude
rdf$y <- rdf$latitude
rdf <- rdf[, c("x", "y", "sst")]
rst <- rasterFromXYZ(rdf)
projection(rst) <- CRS("+proj=longlat +datum=WGS84")
rst <- mask(rst, indiashp2, inverse=TRUE)
field <- mask(field, indiashp2, inverse=TRUE)
field <- mask(field, kerala, inverse=TRUE)
p1 <- vectorplot(field, isField='dXY', region = rst, margin = FALSE, narrows = 75,
length=0.04, main=i, key.arrow = list(size=10, label = 'm/s'), lwd.arrows=1.5,
col.regions=rev(pal(length(at))), at=at,
zlim=zlims, zlimcol="black", colorkey=FALSE, xlab="") +
layer(sp.polygons(indiashp2, col="#585858")) + layer(sp.polygons(kerala, fill="grey", col="#585858")) +
layer(sp.lines(bat.200m, col="grey", alpha=0.9, lty="solid", lwd=2))
p1$par.settings$layout.heights$main.key.padding=-4
p1$main=list(i, x=0.6, just="center")
plist2[[i]] <- p1
}
key <- spplot(rst, scales = list(draw = TRUE),
colorkey = list(space = "bottom"), at=at, zlim=zlims, col.regions=rev(pal(length(at))))
key <- draw.colorkey(key$legend[[1]]$args$key)
plist2[[3]] <- key
key.label <- textGrob("Temperature (°C)", x = unit(0.02, "npc"), y = unit(0.6, "npc"), just = "left")
plist2[[4]] <- key.label
lay <- rbind(c(1,1,2,2),
c(NA,3,3,4))
gridExtra::grid.arrange(grobs=plist2[1:4], layout_matrix = lay, widths = c(1,1,1,1),
heights=c(10,2))
plist2 <- list()
at=seq(26,31,0.1)
zlims=c(26,31)
pal <- colorRampPalette(c("red", "white", "blue"))
for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){
field <- subset(windmon.season, Season==i)[, c("longitude", "latitude", "u", "v")]
field <- rasterFromXYZ(field)
projection(field) <- CRS("+proj=longlat +datum=WGS84")
rdf <- subset(sstmon.season, Season==i)
rdf$x <- rdf$longitude
rdf$y <- rdf$latitude
rdf <- rdf[, c("x", "y", "sst")]
rst <- rasterFromXYZ(rdf)
projection(rst) <- CRS("+proj=longlat +datum=WGS84")
rst <- mask(rst, indiashp2, inverse=TRUE)
field <- mask(field, indiashp2, inverse=TRUE)
field <- mask(field, kerala, inverse=TRUE)
p1 <- vectorplot(field, isField='dXY', region = rst, margin = FALSE, narrows = 75,
length=0.04, main=i, key.arrow = list(size=10, label = 'm/s'), lwd.arrows=1.5,
col.regions=rev(pal(length(at))), at=at,
zlim=zlims, zlimcol="black", colorkey=FALSE, xlab="", ylab="") +
layer(sp.polygons(indiashp2, col="#585858")) + layer(sp.polygons(kerala, fill="grey", col="#585858")) +
layer(sp.lines(bat.200m, col="blue", lty="solid", lwd=2))
p1$par.settings$layout.heights$main.key.padding=-4
p1$main=list(paste("Winds and SST", i), x=0.6, just="center")
plist2[[i]] <- p1
}
key <- spplot(rst, scales = list(draw = TRUE),
colorkey = list(space = "right", height=0.7), at=at, zlim=zlims, col.regions=rev(pal(length(at))))
key <- draw.colorkey(key$legend[[1]]$args$key)
plist2[[3]] <- key
key.label <- textGrob("Temperature (°C)", x = unit(0.02, "npc"), y = unit(0.6, "npc"), just = "left")
plist2[[4]] <- key.label
lay <- rbind(c(1,2,3))
gridExtra::grid.arrange(grobs=plist2[1:3], layout_matrix = lay, widths = c(10,10,1.5),
heights=c(10))
plist3 <- list()
windmon.season$EMT[windmon.season$EMT< -4000] <- -4000
at=seq(-4000,50,50)
pal <- colorRampPalette(c("white", "green", "blue", "black", "red"))
for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){
field <- subset(currentsmon.season, Season==i)[, 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*60, isField='dXY', margin = FALSE, narrows = 500, region=FALSE,
length=0.04, main=paste("Currents", i), key.arrow = list(size=10, label = 'm/min'), lwd.arrows=1.5, xlab="", ylab="") +
layer(sp.polygons(indiashp2, col="#585858")) + layer(sp.polygons(kerala, fill="grey", col="#585858"), under=TRUE) +
layer(sp.lines(bat.200m, col="blue", lty="solid", lwd=2))
p1$par.settings$layout.heights$main.key.padding=-4
p1$main=list(paste("Currents", i), x=0.6, just="center")
plist3[[i]] <- p1
}
gridExtra::grid.arrange(grobs=plist3[c(1,2)], ncol=2)
lay <- rbind(c(1,2,3),c(4,5,NA))
test=gridExtra::arrangeGrob(grobs=c(plist2[1:3], plist3[1:2]), layout_matrix = lay, widths = c(10,10,1.5),
heights=c(10,10), padding=unit(-2,"line"))
grid.draw(test)
grid.text("SST (°C)", x=unit(0.955, "npc"), y=unit(0.56, "npc"), rot=0)
eeholmes/SardineForecast documentation built on July 17, 2021, 2:56 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# to make the 300dpi for submission # delete the eval=FALSE # knit to html # look in Figure_files/figure_html for the eps files 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) library(gridExtra) library(grid) library(marmap)
Get the data
# Get monthly sst data ERA5 #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") sstmon <- getdata("hawaii_soest_d124_2bb9_c935", date=c("1979-01-01", "2020-06-01"), lat=c(7,15), lon=c(70,78), pars="sst", eserver="http://apdrc.soest.hawaii.edu/erddap") sstmon$time <- as.Date(sstmon$time) sstmon$Month <- as.numeric(format(sstmon$time, "%m")) sstmon$Year <- as.numeric(format(sstmon$time, "%Y")) sstmon$sst <- sstmon$sst - 273.15
# Get monthly sst data ERA5 #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 <- getEMT(windmon, coast_angle=158) windmon$EMTperpx <- cos(pi*(180-158)/180)*windmon$EMTperp windmon$EMTperpy <- sin(pi*(180-158)/180)*windmon$EMTperp
# Get monthly ocean currents data currentsmon <- getdata("erdTAgeomday", lat=c(7,15), lon=c(70,78), altitude=0) currentsmon$time <- as.Date(currentsmon$time) currentsmon$Month <- as.numeric(format(currentsmon$time, "%m")) currentsmon$Year <- as.numeric(format(currentsmon$time, "%Y")) currentsmon.mean <- currentsmon %>% group_by(Month, latitude, longitude) %>% summarize(u = mean(u_current, na.rm=TRUE), v = mean(v_current, na.rm=TRUE)) currentsmon.mean$Season = currentsmon.mean$Month currentsmon.mean$Season[currentsmon.mean$Month %in% 6:9] <- "Summer Jun-Sep" currentsmon.mean$Season[currentsmon.mean$Month %in% 1:4] <- "Winter Jan-Apr" currentsmon.mean$Season[currentsmon.mean$Month %in% 10:12] <- "Fall Oct-Dec" currentsmon.season <- currentsmon.mean %>% group_by(Season, latitude, longitude) %>% summarize(u = mean(u, na.rm=TRUE), v = mean(v, na.rm=TRUE))
Make the mean dataframes
windmon.mean <- windmon %>% group_by(Month, latitude, longitude) %>% summarize(v = mean(v, na.rm=TRUE), u = mean(u, na.rm=TRUE), EMT = mean(EMTperp, na.rm=TRUE), We = mean(We, na.rm=TRUE), EMTx = mean(EMTperpx, na.rm=TRUE), EMTy = mean(EMTperpy, 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), EMT = mean(EMT, na.rm=TRUE), We = mean(We, na.rm=TRUE), EMTx = mean(EMTx, na.rm=TRUE), EMTy = mean(EMTy, na.rm=TRUE)) windmon$Season = windmon$Month windmon$Season[windmon$Month %in% 6:9] <- "Summer Jun-Sep" windmon$Season[windmon$Month %in% 1:4] <- "Winter Jan-Apr" windmon.season.2014 <- windmon %>% subset(Year==2014) %>% group_by(Season, latitude, longitude) %>% summarize(v = mean(v, na.rm=TRUE), u = mean(u, na.rm=TRUE), We = mean(We, na.rm=TRUE), EMT = mean(EMTperp, na.rm=TRUE)) sstmon.mean <- sstmon %>% group_by(Month, latitude, longitude) %>% summarize(sst = mean(sst, na.rm=TRUE)) sstmon.mean$Season = sstmon.mean$Month sstmon.mean$Season[sstmon.mean$Month %in% 6:9] <- "Summer Jun-Sep" sstmon.mean$Season[sstmon.mean$Month %in% 1:4] <- "Winter Jan-Apr" sstmon.season <- sstmon.mean %>% group_by(Season, latitude, longitude) %>% summarize(sst = mean(sst, na.rm=TRUE))
Monthly EMT
for(i in 1:12){ rdf <- subset(windmon.mean, Month==i) rdf$x <- rdf$longitude rdf$y <- rdf$latitude rdf <- rdf[, c("x", "y", "We")] if(i==1){ rst <- rasterFromXYZ(rdf); names(rst) <- month.name[i] }else{ tmp <- rasterFromXYZ(rdf); names(tmp) <- month.name[i] rst <- stack(rst, tmp) } } proj4string(rst) <- "+proj=longlat +datum=WGS84" at=seq(-3000,0,50) pal <- colorRampPalette(c("white","red", "green", "blue")) spplot(rst, col.regions=rev(pal(length(at))), at=at, zlim=c(-3000,0), scales=list(draw = TRUE), main="Average Monthly EMT perp to coast (+ in, - out)") + layer(sp.polygons(indiashp2, fill="white"))
Seasonal EMT
for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){ rdf <- subset(windmon.season.2014, Season==i) rdf$x <- rdf$longitude rdf$y <- rdf$latitude rdf <- rdf[, c("x", "y", "EMT")] if(i=="Summer Jun-Sep"){ rst <- rasterFromXYZ(rdf); names(rst) <- i }else{ tmp <- rasterFromXYZ(rdf); names(tmp) <- i rst <- stack(rst, tmp) } } at=seq(-3000,0,50) pal <- colorRampPalette(c("white","red", "green", "blue", "black")) spplot(rst, col.regions=rev(pal(length(at))), at=at, zlim=c(-3000,0), axes=TRUE, main="Average Monthly EMT perp to coast (+ in, - out)") + layer(sp.polygons(indiashp2))
Wind with EMT overlaid
plist2 <- list() windmon.season$EMT[windmon.season$EMT< -4000] <- -4000 at=seq(-4000,50,50) pal <- colorRampPalette(c("white", "green", "blue", "black", "red")) for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){ field <- subset(windmon.season, Season==i)[, c("longitude", "latitude", "u", "v")] field <- rasterFromXYZ(field) projection(field) <- CRS("+proj=longlat +datum=WGS84") rdf <- subset(windmon.season, Season==i) rdf$x <- rdf$longitude rdf$y <- rdf$latitude rdf <- rdf[, c("x", "y", "EMT")] rst <- rasterFromXYZ(rdf) projection(rst) <- CRS("+proj=longlat +datum=WGS84") rst <- mask(rst, indiashp2, inverse=TRUE) field <- mask(field, indiashp2, inverse=TRUE) field <- mask(field, kerala, inverse=TRUE) p1 <- vectorplot(field, isField='dXY', region = rst, margin = FALSE, narrows = 75, length=0.04, main=i, key.arrow = list(size=10, label = 'm/s'), lwd.arrows=1.5, col.regions=rev(pal(length(at))), at=at, zlim=c(-3000,50), zlimcol="black", colorkey=FALSE, xlab="") + layer(sp.polygons(indiashp2, col="#585858")) + 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 } key <- spplot(rst, scales = list(draw = TRUE), colorkey = list(space = "bottom"), at=at, zlim=c(-3000,50), col.regions=rev(pal(length(at)))) key <- draw.colorkey(key$legend[[1]]$args$key) plist2[[3]] <- key key.label <- grid::textGrob(expression(kg~m^-1~s^-1), x = unit(0.02, "npc"), y = unit(0.6, "npc"), just = "left") plist2[[4]] <- key.label
lay <- rbind(c(1,1,2,2), c(NA,3,3,4)) gridExtra::grid.arrange(grobs=plist2[1:4], layout_matrix = lay, widths = c(1,1,1,1), heights=c(10,2))
Wind with SST overlaid
library(marmap) bat <- getNOAA.bathy(lon1 = 70, lon2 = 78, lat1 = 7, lat2 = 15, resolution = 4) save(kerala, indiashp, indiashp2, bat, file=file.path(here::here(), "inst", "extdata", "get_satellite_data", "Replicate_EMT", "indiashp.RData"))
plist2 <- list() at=seq(26,31,0.1) zlims=c(26,31) pal <- colorRampPalette(c("red", "white", "blue")) bat.200m <- rasterToContour(as.raster(bat), levels= -200) for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){ field <- subset(windmon.season, Season==i)[, c("longitude", "latitude", "u", "v")] field <- rasterFromXYZ(field) projection(field) <- CRS("+proj=longlat +datum=WGS84") rdf <- subset(sstmon.season, Season==i) rdf$x <- rdf$longitude rdf$y <- rdf$latitude rdf <- rdf[, c("x", "y", "sst")] rst <- rasterFromXYZ(rdf) projection(rst) <- CRS("+proj=longlat +datum=WGS84") rst <- mask(rst, indiashp2, inverse=TRUE) field <- mask(field, indiashp2, inverse=TRUE) field <- mask(field, kerala, inverse=TRUE) p1 <- vectorplot(field, isField='dXY', region = rst, margin = FALSE, narrows = 75, length=0.04, main=i, key.arrow = list(size=10, label = 'm/s'), lwd.arrows=1.5, col.regions=rev(pal(length(at))), at=at, zlim=zlims, zlimcol="black", colorkey=FALSE, xlab="") + layer(sp.polygons(indiashp2, col="#585858")) + layer(sp.polygons(kerala, fill="grey", col="#585858")) + layer(sp.lines(bat.200m, col="grey", alpha=0.9, lty="solid", lwd=2)) p1$par.settings$layout.heights$main.key.padding=-4 p1$main=list(i, x=0.6, just="center") plist2[[i]] <- p1 } key <- spplot(rst, scales = list(draw = TRUE), colorkey = list(space = "bottom"), at=at, zlim=zlims, col.regions=rev(pal(length(at)))) key <- draw.colorkey(key$legend[[1]]$args$key) plist2[[3]] <- key key.label <- textGrob("Temperature (°C)", x = unit(0.02, "npc"), y = unit(0.6, "npc"), just = "left") plist2[[4]] <- key.label
lay <- rbind(c(1,1,2,2), c(NA,3,3,4)) gridExtra::grid.arrange(grobs=plist2[1:4], layout_matrix = lay, widths = c(1,1,1,1), heights=c(10,2))
plist2 <- list() at=seq(26,31,0.1) zlims=c(26,31) pal <- colorRampPalette(c("red", "white", "blue")) for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){ field <- subset(windmon.season, Season==i)[, c("longitude", "latitude", "u", "v")] field <- rasterFromXYZ(field) projection(field) <- CRS("+proj=longlat +datum=WGS84") rdf <- subset(sstmon.season, Season==i) rdf$x <- rdf$longitude rdf$y <- rdf$latitude rdf <- rdf[, c("x", "y", "sst")] rst <- rasterFromXYZ(rdf) projection(rst) <- CRS("+proj=longlat +datum=WGS84") rst <- mask(rst, indiashp2, inverse=TRUE) field <- mask(field, indiashp2, inverse=TRUE) field <- mask(field, kerala, inverse=TRUE) p1 <- vectorplot(field, isField='dXY', region = rst, margin = FALSE, narrows = 75, length=0.04, main=i, key.arrow = list(size=10, label = 'm/s'), lwd.arrows=1.5, col.regions=rev(pal(length(at))), at=at, zlim=zlims, zlimcol="black", colorkey=FALSE, xlab="", ylab="") + layer(sp.polygons(indiashp2, col="#585858")) + layer(sp.polygons(kerala, fill="grey", col="#585858")) + layer(sp.lines(bat.200m, col="blue", lty="solid", lwd=2)) p1$par.settings$layout.heights$main.key.padding=-4 p1$main=list(paste("Winds and SST", i), x=0.6, just="center") plist2[[i]] <- p1 } key <- spplot(rst, scales = list(draw = TRUE), colorkey = list(space = "right", height=0.7), at=at, zlim=zlims, col.regions=rev(pal(length(at)))) key <- draw.colorkey(key$legend[[1]]$args$key) plist2[[3]] <- key key.label <- textGrob("Temperature (°C)", x = unit(0.02, "npc"), y = unit(0.6, "npc"), just = "left") plist2[[4]] <- key.label
lay <- rbind(c(1,2,3)) gridExtra::grid.arrange(grobs=plist2[1:3], layout_matrix = lay, widths = c(10,10,1.5), heights=c(10))
plist3 <- list() windmon.season$EMT[windmon.season$EMT< -4000] <- -4000 at=seq(-4000,50,50) pal <- colorRampPalette(c("white", "green", "blue", "black", "red")) for(i in c("Summer Jun-Sep", "Winter Jan-Apr")){ field <- subset(currentsmon.season, Season==i)[, 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*60, isField='dXY', margin = FALSE, narrows = 500, region=FALSE, length=0.04, main=paste("Currents", i), key.arrow = list(size=10, label = 'm/min'), lwd.arrows=1.5, xlab="", ylab="") + layer(sp.polygons(indiashp2, col="#585858")) + layer(sp.polygons(kerala, fill="grey", col="#585858"), under=TRUE) + layer(sp.lines(bat.200m, col="blue", lty="solid", lwd=2)) p1$par.settings$layout.heights$main.key.padding=-4 p1$main=list(paste("Currents", i), x=0.6, just="center") plist3[[i]] <- p1 }
gridExtra::grid.arrange(grobs=plist3[c(1,2)], ncol=2)
lay <- rbind(c(1,2,3),c(4,5,NA)) test=gridExtra::arrangeGrob(grobs=c(plist2[1:3], plist3[1:2]), layout_matrix = lay, widths = c(10,10,1.5), heights=c(10,10), padding=unit(-2,"line")) grid.draw(test) grid.text("SST (°C)", x=unit(0.955, "npc"), y=unit(0.56, "npc"), rot=0)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.