## ROMS Atlantis coupling MDSumner, Jessica Melbourne-Thomas April 2017
## NOTES
#https://confluence.csiro.au/display/Atlantis/Hydro+FAQ
#https://github.com/cecilieha/NoBA/blob/master/create_flux_ncdf_atlantis.R
##corners <- read.delim("https://raw.githubusercontent.com/cecilieha/NoBA/master/corners_neighbours_nordic.txt")
## curvlinear vector rotation: https://www.myroms.org/forum/viewtopic.php?f=3&t=295
## hydro stuff
# https://confluence.csiro.au/display/Atlantis/Forcing+Files
# https://confluence.csiro.au/display/Atlantis/Current+Forcing+File+Structure
# https://confluence.csiro.au/display/Atlantis/HydroConstruct
# https://confluence.csiro.au/display/Atlantis/Matlab+-+Hydro
## dummy hydro
# sohydro <- "~/projects/AtlantisEastAntartica/SOhydrodummy.nc"
# library(ncdump)
# nc <- NetCDF(sohydro)
#' Convenience function to transform map projection .
#'
#' Transform `x` to whatever the projection of `to` is.
#' @param x object to transform
#' @param to object with a map projection
#'
#' @return `x`, transformed
#' @export
#' @importFrom raster projection
#' @importFrom sp spTransform
project_to <- function(x, to) {
spTransform(x, CRS(projection(to)))
}
## devtools::install_github(c("mdsumner/angstroms"))
library(angstroms)
library(rbgm) ## read BGM
library(bgmfiles) ## archive of BGM files
library(raadtools)
library(ncdump)
library(dplyr)
## get the Circumpolar files
cpolar <- raadtools:::cpolarfiles()
## this single dataframe records every time slice across
## all files in a row, so we can loop over time step below
file_db <- bind_rows(lapply(cpolar$fullname, function(x) {
nc <- NetCDF(x)
tlen <- filter(nc$dimension, name == "ocean_time")$len
tibble(fullname = rep(x, tlen), band_level = seq_len(tlen))
}))
roms_file <- cpolar$fullname[2]
## get a BGM and read it
bfile <- bgmfiles::bgmfiles("antarctica_28")
bgm <- bgmfile(bfile)
## get the longitude/latitude arrays
roms_ll <- romscoords(roms_file, transpose = TRUE)
## we need the unsullied boxes to identify points inside them
boxes <- boxSpatial(bgm)
roms_face <- romsmap(project_to(faceSpatial(bgm), "+init=epsg:4326"), roms_ll)
## which box does each point fall in
index_box <- function(box_sp, roms_ll) {
ind <- sp::over(project_to(coords_points(roms_ll), box_sp) , as(box_sp, "SpatialPolygons"))
# tibble(box = ind,
tibble(box = box_sp$label[ind],
cell = seq_len(ncell(roms_ll))) %>%
filter(!is.na(box))
}
## build the index for each box to the ROMS cells it contains
## and each face for the ROMS cells it traverses
library(dplyr)
box_roms_index <- index_box(boxes, roms_ll)
ind_face <- cellFromLine(romsdata(roms_file, "u"), roms_face)
face_roms_index <- tibble(face = roms_face$label[rep(seq_len(nrow(roms_face)), lengths(ind_face))],
cell = unlist(ind_face))
#' return the ramp of positive depths from surface down
#' (so that the order is native to the NetCDF order)
roms_level <- function(Cs_r, h, cell) {
extract(h, cell) * Cs_r
}
## important to readAll here, else extract is very slow in the loop
h <- readAll(raster(roms_file, varname = "h"))
## Cs_r is the S-coord stretching
Cs_r <- rawdata(roms_file, "Cs_r")
## build the level index between Atlantis and ROMS
list_nc_z_index <- vector('list', nrow(box_roms_index))
max_depth <- max(extract(h, unique(box_roms_index$cell)))
atlantis_depths <- -rev(cumsum(c(0, rev(rbgm::build_dz(-max_depth)))))
for (i in seq_len(nrow(box_roms_index))) {
rl <- roms_level(Cs_r, h, box_roms_index$cell[i])
## implicit 0 at the surface, and implicit bottom based on ROMS
list_nc_z_index[[i]] <- length(atlantis_depths) - findInterval(rl, atlantis_depths) + 1
if (i %% 1000 == 0) print(i)
}
#
# ll <- extract(roms_ll, box_roms_index$cell)
# box_roms_index$lon <- ll[, 1]
# box_roms_index$lat <- ll[, 2]
# library(ggplot2)
# ggplot(box_roms_index, aes(lon, lat, colour = box)) + geom_point(pch = ".")
#
## join the box-xy-index to the level index
box_z_index <- bind_rows(lapply(list_nc_z_index,
function(x) tibble(atlantis_level = pmax(1, x), roms_level = seq_along(x))),
.id = "cell_index") %>%
inner_join(mutate(box_roms_index, cell_index = as.character(row_number()))) %>%
select(-cell_index)
# ll <- extract(roms_ll, box_z_index$cell)
# box_z_index$lon <- ll[, 1]
# box_z_index$lat <- ll[, 2]
# library(ggplot2)
# ggplot(box_z_index, aes(lon, lat, colour = roms_level)) + geom_point(pch = ".") +
# facet_wrap(~atlantis_level)
for (i in seq_len(nrow(face_roms_index))) {
rl <- roms_level(Cs_r, h, face_roms_index$cell[i])
## implicit 0 at the surface, and implicit bottom based on ROMS
list_nc_z_index[[i]] <- length(atlantis_depths) -findInterval(rl, atlantis_depths) + 1
if (i %% 1000 == 0) print(i)
}
## join the face-xy-index to the level index
face_z_index <- bind_rows(lapply(list_nc_z_index,
function(x) tibble(atlantis_level = x, roms_level = seq_along(x))),
.id = "cell_index") %>%
inner_join(mutate(face_roms_index, cell_index = as.character(row_number()))) %>%
select(-cell_index)
#z_index$atlantis_depth <- c(0, rev(rbgm::build_dz(box_roms_index$botz[i])))[z_index$atlantis_level]
#z_index$roms_depth <- rl[z_index$roms_level]
## driver function to loop over x raster by levels
## matching $cell for the right group
extract_at_level <- function(x, cell_level) {
ulevel <- unique(cell_level$level)
values <- numeric(nrow(cell_level))
for (ul in seq_along(ulevel)) {
asub <- cell_level$level == ulevel[ul]
values[asub] <- extract(x[[ulevel[ul]]],
cell_level$cell[asub])
}
values
}
#file_db <- file_db[1, ]
box_props <- face_props <- vector("list", nrow(file_db))
i_timeslice <- 1
#for (i_timeslice in seq(nrow(file_db))) {
roms_file <- file_db$fullname[i_timeslice]
level <- file_db$band_level[i_timeslice]
ru <- set_indextent(brick(roms_file, varname = "u", lvar = 4, level = level))
face_z_index$u <- extract_at_level(ru, rename(face_z_index, level = roms_level, cell = cell))
rv <- set_indextent(brick(roms_file, varname = "v", lvar = 4, level = level))
face_z_index$v <- extract_at_level(rv, rename(face_z_index, level = roms_level, cell = cell))
rdata <- set_indextent(brick(roms_file, varname = "temp", lvar = 4, level = level))
box_z_index$temp <- extract_at_level(rdata, rename(box_z_index, level = roms_level, cell = cell))
rdata <- set_indextent(brick(roms_file, varname = "salt", lvar = 4, level = level))
box_z_index$salt <- extract_at_level(rdata, rename(box_z_index, level = roms_level, cell = cell))
box_props[[i_timeslice]] <- box_z_index %>% group_by(atlantis_level, box) %>%
summarize(temp = mean(temp, na.rm = TRUE), salt = mean(salt ,na.rm = TRUE)) %>%
mutate(band_level = level)
face_props[[i_timeslice]] <- face_z_index %>% group_by(atlantis_level, face) %>%
summarize(flux = mean(sqrt(u * u + v * v), na.rm = TRUE) * 24 * 3600) %>%
mutate(band_level = level)
#}
box_props <- bind_rows(box_props)
face_props <- bind_rows(face_props)
library(ggplot2)
ggplot(box_props, aes(x = temp, y = salt, colour = atlantis_level)) + geom_point() + facet_wrap(~box)
ggplot(face_props, aes(x = atlantis_level, y = flux, colour = atlantis_level)) + geom_point()
#
#
# plot(romsmap(boxes[15, ], roms_ll), add = F, border = "red", lwd = 4)
# plot(sqrt(ru[[10]]^2 + rv[[10]]^2), add = TRUE,
# col = colorRampPalette(c("#132B43", high = "#56B1F7"))(10), zlim = c(0, 0.1))
# plot(romsmap(boxes[15, ], roms_ll), add = TRUE, border = "red", lwd = 4)
# ggplot(mutate(z_index %>%
# filter(box == 15),
# x = xFromCell(ru, cell), y = yFromCell(ru, cell)),
# aes(x = x, y = y, fill = sqrt(u * u + v * v))) + geom_tile() + coord_equal() +
# facet_wrap(~roms_level)
#
#
# rm(ru, rv)
#
#
# ggplot(z_index %>% mutate(x = xFromCell(rdata, cell),
# y = yFromCell(rdata, cell)),
# aes(x = x, y = y, fill = temp)) + geom_tile() +
# coord_equal()
#
# ggplot(mutate(z_index %>%
# filter(roms_level %% 3 == 0),
# x = xFromCell(ru, cell), y = yFromCell(ru, cell)),
# aes(x = x, y = y, fill = temp)) + geom_tile() + coord_equal() +
# facet_wrap(~roms_level)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.