#' @title Load Coastline Data
#' @description Low res to high res: coastlineWorld, coastlineWorldCoarse, coastlinewWorldMedium, coastlineWorldFine
#' @author Laura Whitmore
#' @author Thomas Bryce Kelly
#' @param cost a string referring to a coastline data object such as those provided by the OCE package.
#' @import oce
#' @import ocedata
get.coast = function(coast = 'coastlineWorld') {
do.call('data', list(coast))
}
#' @title Make Map Projection
#' @author Thomas Bryce Kelly
#' @param projection A string or number corresponding to a base projection (e.g. 1 = 'merc')
#' @param lat the primary latitude for the projection, may or may not be applicable based on projection
#' @param lon same as lat for longitude
#' @param h The viewing height in meters (used only for some projections)
#' @param dlat The distance to the secondary latitude (in degrees). Only applicable to some projections
#' @export
make.proj = function(projection = NULL, lat = NULL, lon = NULL, h = NULL, dlat = 10) {
if (is.null(projection)) { projection = 'merc' }
projections = c('merc', 'aea', 'eck3', 'eck4', 'eqc', 'geos', 'lonlat', 'mill', 'natearth', 'nsper', 'stere', 'ortho')
projection.list = c('1) merc', '2) aea', '3) eck3', '4) eck4', '5) eqc', '6) geos', '7) lonlat', '8) mill', '9) natearth', '10) nsper', '11) stere', '12) ortho')
if (is.numeric(projection)) {
projection = projections[projection]
}
if (!projection %in% projections) {
message('Unknown projection type, recommend to use:', paste(projection.list, collapse = ', '))
}
if (projection == 'geos' & is.null(h)) { h = 1e8 }
if (projection == 'nsper' & is.null(h)) { h = 1e8 }
# Default
if (is.null(lat)) { lat = 0 }
if (is.null(lon)) { lon = 0 }
if (is.null(h)) { h = 1e8 }
paste0('+proj=', projection, ' +lon_0=', lon, ' +lat_0=', lat, ' +lat_1=', lat, ' +lat_2=', lat + dlat, ' +h=', h)
}
#' @title Bilinear Interpolation of a Grid
#' @author Thomas Bryce Kelly
#' @export
interp.bilinear = function(x, y, gx, gy, z) {
z.out = rep(0, length(x))
for (i in 1:length(x)) {
x1 = max(0, which(gx <= x[i]))
x2 = x1 + 1
y1 = max(0, which(gy <= y[i]))
y2 = y1 + 1
wx1 = (gx[x2] - x[i]) / (gx[2] - gx[1])
wy1 = (gy[y2] - y[i]) / (gy[2] - gy[1])
if (x1 == 0) {
x1 = 1
wx1 = 1
}
if (y1 == 0) {
y1 = 1
wy1 = 1
}
if(x1 == length(gx)) {
x2 = x1
wx1 = 1
}
if(y1 == length(gy)) {
y2 = y1
wy1 = 1
}
z.out[i] = wy1 * (wx1 * z[x1, y1] + (1 - wx1) * z[x2,y1]) + (1 - wy1) * (wx1 * z[x1, y2] + (1 - wx1) * z[x2,y2])
}
z.out
}
#' @title Interpolate from fractional index onto grid
#' @author Thomas Bryce Kelly
#' @export
grid.interp = function(grid, i, j) {
val = rep(NA, length(i))
x = round(i)
y = round(j)
dx = i - x
dy = j - y
for (k in 1:length(i)) {
val[k] = (1 - abs(dy[k])) * ((1 - abs(dx[k])) * grid[cbind(x[k], y[k])] + abs(dx[k]) * grid[cbind(x[k] + sign(dx[k]), y[k])]) +
abs(dy[k]) * ((1 - abs(dx[k])) * grid[cbind(x[k], y[k] + sign(dy[k]))] + abs(dx[k]) * grid[cbind(x[k] + sign(dx[k]), y[k] + sign(dy[k]))])
}
## Return
val
}
#' @title Calculate extended grid
#' @author Thomas Bryce Kelly
#' @export
calc.vertex = function(x, y) {
## Diffs
dx.dx = t(diff(t(x))) / 2
dx.dy = diff(x) / 2
dy.dx = t(diff(t(y))) / 2
dy.dy = diff(y) / 2
## Vertex
vertex.x = matrix(NA, nrow = dim(x)[1]+1, ncol = dim(x)[2]+1)
vertex.y = matrix(NA, nrow = dim(y)[1]+1, ncol = dim(y)[2]+1)
## Field
for (i in 2:(dim(vertex.x)[1] - 1)) {
for (j in 2:(dim(vertex.x)[2] - 1)) {
ii = max(i-1, 1)
jj = max(j-1, 1)
vertex.x[i,j] = x[ii,jj] + dx.dx[ii,jj] + dx.dy[ii,jj]
vertex.y[i,j] = y[ii,jj] + dy.dx[ii,jj] + dy.dy[ii,jj]
}
}
## Fill in perimeter
# i = 1
for (j in 2:(dim(vertex.x)[2] - 1)) {
jj = max(j-1, 1)
vertex.x[1,j] = x[1,jj] + dx.dx[1,jj] - dx.dy[1,jj]
vertex.y[1,j] = y[1,jj] + dy.dx[1,jj] - dy.dy[1,jj]
}
# j = 1
for (i in 2:(dim(vertex.x)[1] - 1)) {
ii = max(i-1, 1)
vertex.x[i,1] = x[ii,1] - dx.dx[ii,1] + dx.dy[ii,1]
vertex.y[i,1] = y[ii,1] - dy.dx[ii,1] + dy.dy[ii,1]
}
# j = dim(vertex.x)[2]
for (i in 1:(dim(vertex.x)[1] - 1)) {
ii = max(i-1, 1)
j = dim(vertex.x)[2]
vertex.x[i,j] = x[ii,j-1] + dx.dx[ii,j-2] + dx.dy[ii,j-1]
vertex.y[i,j] = y[ii,j-1] + dy.dx[ii,j-2] + dy.dy[ii,j-1]
}
# i = dim(vertex.x)[2]
for (j in 1:(dim(vertex.x)[2] - 1)) {
jj = max(j-1, 1)
i = dim(vertex.x)[1]
vertex.x[i,j] = x[i-1,jj] + dx.dx[i-1,jj] + dx.dy[i-2,jj]
vertex.y[i,j] = y[i-1,jj] + dy.dx[i-1,jj] + dy.dy[i-2,jj]
}
## Fill in corners
## both = 1
vertex.x[1,1] = x[1,1] - dx.dx[1,1] - dx.dy[1,1]
vertex.y[1,1] = y[1,1] - dy.dx[1,1] - dy.dy[1,1]
i = dim(vertex.x)[1]
vertex.x[i,1] = x[i-1,1] - dx.dx[i-1,1] + dx.dy[i-2,1]
vertex.y[i,1] = y[i-1,1] - dy.dx[i-1,1] + dy.dy[i-2,1]
i = dim(vertex.x)[1]
j = dim(vertex.x)[2]
vertex.x[i,j] = x[i-1,j-1] + dx.dx[i-1,j-2] + dx.dy[i-2,j-1]
vertex.y[i,j] = y[i-1,j-1] + dy.dx[i-1,j-2] + dy.dy[i-2,j-1]
j = dim(vertex.x)[2]
vertex.x[1,j] = x[1,j-1] + dx.dx[1,j-2] - dx.dy[1,j-1]
vertex.y[1,j] = y[1,j-1] + dy.dx[1,j-2] - dy.dy[1,j-1]
list(x = vertex.x, y = vertex.y)
}
#' @title Calcuye extended grid
#' @author Thomas Bryce Kelly
#' @export
grid.refinement = function(x = NULL, y = NULL, z) {
if (is.null(dim(z))) {stop('grid.refinement: z must be an array object of two dimensions.')}
dim = dim(z)
if (is.null(x) & is.null(y)) {
x = c(1:dim[1])
y = c(1:dim[2])
}
if (is.null(dim(x)) & is.null(dim(y))) {
x = array(x, dim = dim)
y = t(array(y, dim = rev(dim)))
}
## Vertex
vertex.x = array(0, dim = c(2*dim(x)[1]-1, 2*dim(x)[2]-1))
vertex.y = vertex.x
vertex.z = vertex.x
## fill in known values
for (i in 1:dim(x)[1]) {
for (j in 1:dim(x)[2]) {
vertex.x[2*i-1, 2*j-1] = x[i,j]
vertex.y[2*i-1, 2*j-1] = y[i,j]
vertex.z[2*i-1, 2*j-1] = z[i,j]
}
}
## Interpolate x
for (i in 1:(dim(x)[1]-1)) {
for (j in 1:dim(x)[2]) {
vertex.x[2*i, 2*j-1] = 0.5 * (x[i,j] + x[i+1,j])
vertex.y[2*i, 2*j-1] = 0.5 * (y[i,j] + y[i+1,j])
vertex.z[2*i, 2*j-1] = 0.5 * (z[i,j] + z[i+1,j])
}
}
## Interpolate y
for (i in 1:dim(x)[1]) {
for (j in 1:(dim(x)[2]-1)) {
vertex.x[2*i-1, 2*j] = 0.5 * (x[i,j] + x[i,j+1])
vertex.y[2*i-1, 2*j] = 0.5 * (y[i,j] + y[i,j+1])
vertex.z[2*i-1, 2*j] = 0.5 * (z[i,j] + z[i,j+1])
}
}
## corners
for (i in 1:(dim(x)[1]-1)) {
for (j in 1:(dim(x)[2]-1)) {
vertex.x[2*i, 2*j] = 0.25 * (x[i,j] + x[i,j+1] + x[i+1,j] + x[i+1,j+1])
vertex.y[2*i, 2*j] = 0.25 * (y[i,j] + y[i,j+1] + y[i+1,j] + y[i+1,j+1])
vertex.z[2*i, 2*j] = 0.25 * (z[i,j] + z[i,j+1] + z[i+1,j] + z[i+1,j+1])
}
}
list(x = vertex.x, y = vertex.y, z = vertex.z)
}
#' @title Calculate subsampled grid
#' @author Thomas Bryce Kelly
#' @export
grid.subsample = function(x = NULL, y = NULL, z, approx = F) {
if (is.null(dim(z))) {stop('grid.refinement: z must be an array object of two dimensions.')}
dim = dim(z)
if (is.null(x) & is.null(y)) {
x = c(1:dim[1])
y = c(1:dim[2])
}
if (is.null(dim(x)) & is.null(dim(y))) {
x = array(x, dim = dim)
y = t(array(y, dim = rev(dim)))
}
## Vertex
vertex.x = array(0, dim = floor(c(dim(x)[1], dim(x)[2]) / 2))
vertex.y = vertex.x
vertex.z = vertex.x
if (approx) {
for (i in 1:dim(vertex.x)[1]) {
for (j in 1:dim(vertex.x)[2]) {
vertex.x[i, j] = x[2*i, 2*j]
vertex.y[i, j] = y[2*i, 2*j]
vertex.z[i, j] = z[2*i, 2*j]
}
}
} else {
for (i in 1:dim(vertex.x)[1]) {
for (j in 1:dim(vertex.x)[2]) {
vertex.x[i, j] = 0.25 * (x[2*i, 2*j] + x[2*i-1, 2*j] + x[2*i, 2*j-1] + x[2*i-1, 2*j-1])
vertex.y[i, j] = 0.25 * (y[2*i, 2*j] + y[2*i-1, 2*j] + y[2*i, 2*j-1] + y[2*i-1, 2*j-1])
vertex.z[i, j] = mean(c(z[2*i, 2*j], z[2*i-1, 2*j], z[2*i, 2*j-1], z[2*i-1, 2*j-1]), na.rm = T)
}
}
}
list(x = vertex.x, y = vertex.y, z = vertex.z)
}
#' @export
calc.section.dist = function(lon, lat) {
if (length(lon) == 1) { lon = rep(lon, length(lat))}
if (length(lat) == 1) { lat = rep(lat, length(lon))}
if (length(lon) != length(lat)) { stop('Length of lon/lat are not the same!')}
d = rep(NA, length(lon))
l = which(diff(lon) != 0 | diff(lat) != 0)
d[1:(l[1]-1)] = 0 ## same station
if (lenght(l) > 1) {
for (i in 2:length(l)) {
d[l[i-1]:(l[i]-1)] = calc.dist(lon[c(l[i],l[i]+1)], lat[c(l[i],l[i]+1)])
}
}
## return
d
}
#' @export
calc.dist = function(lon, lat) {
sapply(1:(length(lon)-1), function(x) {abs(lon[x+1] - lon[x])})
}
#' @title Retreive depth value from bathymetric grid.
#' @export
get.depth = function(lon, lat, bathy) {
depths = rep(NA, length(lon))
for (i in 1:lnegth(lon)) {
k1 = which.min(abs(lon[i] - bathy$Lon))
k2 = which.min(abs(lat[i] - bathy$Lat))
depths[i] = bathy$Z[k1,k2]
}
## Return
depths
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.