R/map.plot.R
In tbrycekelly/TheSource: TheSource: It Will Hold Your Beer
Defines functions
redraw.map
add.map.bathy
add.map.eez
add.map.layer
add.map.contour
add.map.quiver
add.map.scale
add.map.grid
add.map.axis
add.map.axis.label
add.map.polygon
add.map.text
add.map.arrows
add.map.line
add.map.points
get.map.line
add.map.coastline
Documented in
add.map.arrows
add.map.axis
add.map.axis.label
add.map.bathy
add.map.coastline
add.map.contour
add.map.grid
add.map.layer
add.map.line
add.map.points
add.map.polygon
add.map.quiver
add.map.scale
add.map.text
get.map.line
redraw.map
#' @title Make Map
#' @description Default settings are for Arctic, lowest resolution coastline.
#' Exaples of projection arguments include: `+proj=stere +lat_0=90`, `+proj=merc`, `+proj=aea +lat_1=30 +lat_2=45 +lon_0=-120`.
#' @param coast Should be the name of a coastline data object. A value of NULL sets the default cosatline to 'coastline2'.
#' @param lon.min The minimum longitude displayed on the map.
#' @param lon.max The maximum longitude on the map
#' @param lat.min The minimum latitude shown on the map
#' @param lat.max The maximum latitude shown on the map
#' @param p a projection string, such as those generated by 'make.proj()'
#' @param land.col A color string or value for the land polygon
#' @param grid Boolean, draw lat/lon grid?
#' @param dlon The spacing for the longitude grid (in degrees)
#' @param dlat The spacing for the latitude grid (in degrees)
#' @import oce
#' @import ocedata
#' @author Laura Whitmore
#' @author Thomas Bryce Kelly
#' @export
make.map = function (coast = NULL,
lon.min = -180,
lon.max = 180,
lat.min = -60,
lat.max = 80,
p = NULL,
land.col = 'lightgray',
draw.grid = T,
draw.axis = T,
dlon = 60,
dlat = 20,
verbose = T) {
## Apply Defaults
if (all(is.null(c(lon.min, lon.max, lat.min, lat.max)))) {
lon.min = -180
lon.max = 180
lat.min = -60
lat.max = 80
} else if (any(is.null(c(lon.min, lon.max, lat.min, lat.max)))) {
stop('All lon/lat bounds must be provided!')
}
if (is.null(coast)) { coast = 'coastline2' }
if (is.null(p)) {
p = make.proj(projection = 1,
lat = mean(c(lat.max, lat.min)),
lon = mean(c(lon.min, lon.max)))
}
## Determine lat/lon axis details
lons = seq(lon.min - 5*dlon, lon.max + 5*dlon, by = dlon)
lats = seq(lat.min - 5*dlat, lat.max + 5*dlon, by = dlat)
if (lon.min < 0) {
lons = lons[lons >= -180 & lons <= 180]
} else {
lons = lons[lons >= 0 & lons <= 360]
}
lats = lats[lats >= -90 & lats <= 90]
lims = expand.grid(lon = lons, lats = lats)
lon.axis = rgdal::project(cbind(lons, rep(lat.min, length(lons))), proj = p)[,1]
lat.axis = rgdal::project(cbind(rep(lon.min, length(lats)), lats), proj = p)[,2]
lims = rgdal::project(cbind(lims$lon, lims$lat), proj = p)
## Start plotting
# Calculate projections of plotting area for plot limits
field = expand.grid(lon = seq(lon.min, lon.max, length.out = 10),
lat = seq(lat.min, lat.max, length.out = 10))
field = rgdal::project(cbind(field$lon, field$lat), proj = p)
## make plot
plot(NULL,
NULL,
xlim = range(field[,1], na.rm = T),
ylim = range(field[,2], na.rm = T),
xaxt = 'n', yaxt = 'n',
xlab = '', ylab = '',
xaxs = 'i', yaxs = 'i')
coast = add.map.coastline(coast, p = p, land.col = land.col)
## Setup map object
map = list(coastline = coast,
lon = (lon.min + lon.max)/2,
lat = (lat.min + lat.max)/2,
scale = NA,
p = p,
land.col = land.col,
grid = list(
draw.grid = draw.grid,
draw.axis = draw.axis,
dlon = dlon,
dlat = dlat
),
meta = list(
Source.version = packageVersion('TheSource'),
R.version = R.version.string,
Sys.time = Sys.time()
))
if (draw.grid) { add.map.grid(map) }
if (draw.axis) { add.map.axis(map, sides = c(1,2)) }
box()
map
}
#' @title Add Coastline to Map
#' @author Thomas Bryce Kelly
#' @description Adds a coastline list to a map given a projection and land color.
#' @export
add.map.coastline = function(coast, p, land.col, field.inv = NULL) {
## helper function
break.polygon = function(projected.coast) {
for (i in 1:length(projected.coast)) {
k = c(which(diff(projected.coast[[i]][,1])^2 + diff(projected.coast[[i]][,2])^2 > diff(par('usr')[1:2])^2 / 10), length(projected.coast[[i]][,1]))
if (length(k) > 1) {
for (j in 2:length(k)) {
## split off parts of polygons with more than 5 points
if (k[j] - k[j-1] > 5) {
projected.coast[[length(projected.coast) + 1]] = projected.coast[[i]][(k[j-1]+1):k[j],]
projected.coast[[i]][(k[j-1]+1):k[j],] = NA ## Remove newly split areas
}
}
## Trim original down to just first likely contiguous area
projected.coast[[i]] = array(projected.coast[[i]][1:k[1],], dim = c(k[1], 2))
}
}
## Return
projected.coast
}
# Load data if necessary
if (typeof(coast) != "character") {
coastline = coast
} else {
do.call('data', list(coast))
coastline = eval(parse(text = coast))
}
## First pass at trimming domain
if (!is.null(field.inv)) {
keep = rep(T, length(coastline$data))
dx = median(diff(field.inv[,1])) * 2
dy = median(diff(field.inv[,2])) * 2
for (i in 1:length(coastline$data)) {
if (all(coastline$data[[i]][,1] + dx < min(field.inv[,1], na.rm = T)) | all(coastline$data[[i]][,1] - dx > max(field.inv[,1], na.rm = T)) |
all(coastline$data[[i]][,2] + dy < min(field.inv[,2], na.rm = T)) | all(coastline$data[[i]][,2] - dy > max(field.inv[,2]), na.rm = T)) {
keep[i] = F
}
}
if (sum(keep) == 0) {
return(list(data = NA, meta = coastline$meta, projected = NA))
}
coastline$data = coastline$data[keep]
coastline$meta$trim1 = sum(keep)
}
## Project coastline (takes a while!)
projected.coast = lapply(coastline$data, function(x) {
rgdal::project(as.matrix(x)+1e-16, proj = p)
})
## Second round of trimming
keep = rep(T, length(projected.coast))
usr = par('usr')
for (i in 1:length(projected.coast)) {
if (all(projected.coast[[i]][,1] < usr[1]) | all(projected.coast[[i]][,1] > usr[2]) |
all(projected.coast[[i]][,2] < usr[3]) | all(projected.coast[[i]][,2] > usr[4])) {
keep[i] = F
}
}
if (sum(keep) == 0) {
return(list(data = NA, meta = coastline$meta, projected = NA))
}
coastline$data = coastline$data[keep]
projected.coast = projected.coast[keep]
projected.coast = break.polygon(projected.coast)
coastline$meta$trim2 = sum(keep)
## Add coastlines (takes a while!)
for (i in 1:length(projected.coast)) {
polygon(projected.coast[[i]][,1], projected.coast[[i]][,2], col = land.col)
}
list(data = coastline$data, meta = coastline$meta, projected = projected.coast)
}
#' @title Make Map
#' @description Default settings are for Arctic, lowest resolution coastline.
#' Exaples of projection arguments include: `+proj=stere +lat_0=90`, `+proj=merc`, `+proj=aea +lat_1=30 +lat_2=45 +lon_0=-120`.
#' @param coast Should be the name of a coastline data object. A value of NULL sets the default cosatline to 'coastline2'.
#' @param lon.min The minimum longitude displayed on the map.
#' @param lon.max The maximum longitude on the map
#' @param lat.min The minimum latitude shown on the map
#' @param lat.max The maximum latitude shown on the map
#' @param p a projection string, such as those generated by 'make.proj()'
#' @param land.col A color string or value for the land polygon
#' @param grid Boolean, draw lat/lon grid?
#' @param dlon The spacing for the longitude grid (in degrees)
#' @param dlat The spacing for the latitude grid (in degrees)
#' @import oce
#' @import ocedata
#' @author Laura Whitmore
#' @author Thomas Bryce Kelly
#' @export
make.map2 = function (coast = NULL,
lon = 0,
lat = 0,
scale = 1000,
p = NULL,
land.col = 'lightgray',
draw.grid = T,
draw.axis = T,
dlon = 60,
dlat = 20,
verbose = T) {
## Apply Defaults
if (is.null(coast)) { coast = 'coastline2' }
if (is.null(p)) {
p = make.proj(projection = 1,
lat = lat,
lon = lon)
}
center.point = rgdal::project(cbind(lon, lat), proj = p)
x.y.ratio = par()$pin[1] / par()$pin[2]
## Start plotting
# Calculate projections of plotting area for plot limits
field = expand.grid(lon = seq(center.point[,1] - scale * x.y.ratio * 1e3, center.point[,1] + scale * x.y.ratio * 1e3, length.out = 10),
lat = seq(center.point[,2] - scale * 1e3, center.point[,2] + scale * 1e3, length.out = 10))
field.inv = rgdal::project(cbind(field$lon, field$lat), proj = p, inv = T)
## make plot
plot(NULL,
NULL,
xlim = range(field[,1], na.rm = T),
ylim = range(field[,2], na.rm = T),
xaxt = 'n', yaxt = 'n',
xlab = '', ylab = '',
xaxs = 'i', yaxs = 'i')
coast = add.map.coastline(coast, p = p, land.col = land.col, field.inv = field.inv)
## Setup map object
map = list(coastline = coast,
lon = lon,
lat = lat,
scale = scale,
p = p,
land.col = land.col,
grid = list(
draw.grid = draw.grid,
draw.axis = draw.axis,
dlon = dlon,
dlat = dlat
),
meta = list(
Source.version = packageVersion('TheSource'),
R.version = R.version.string,
Sys.time = Sys.time()
))
if (draw.grid) { add.map.grid(map) }
if (draw.axis) { add.map.axis(map, sides = c(1,2)) }
box()
map
}
#' @title Get Map Line
#' @author Thomas Bryce Kelly
#' @import rgdal
#' @param lon a set of longitudes to draw a line between
#' @param lat a set of latitudes to draw a lien between
#' @param col the color of the line to be drawn
#' @param lty the type of line to be drawn
#' @param lty the type of line to be drawn
#' @param lwd the width of the line to be drawn
#' @export
get.map.line = function(map,
lon,
lat,
greatCircle = F,
N = 1e3) {
if (length(lon) != length(lat)) {stop('get.map.line: length of lon/lat are not the same.')}
if (greatCircle) {
lons = approx(1:length(lon), lon, xout = seq(1, length(lon), length.out = max(length(lon), 1e4)))$y
lats = approx(1:length(lat), lat, xout = seq(1, length(lat), length.out = max(length(lon), 1e4)))$y
} else {
lons = lon
lats = lat
}
## project
rgdal::project(cbind(lons, lats), proj = map$p)
}
#' @title Add Map Points
#' @author Laura Whitmore
#' @description Add station points to a map
#' @param lon the longitudes of the points to be drawn
#' @param lat the latitude of the points to be drawn
#' @param stn.lon (depreciated) the longitudes of the points to be drawn
#' @param stn.lat (depreciated) the latitude of the points to be drawn
#' @param col the color of the points
#' @param cex the size of point to be drawn
#' @param pch the point character to be used
#' @param ... optional arguments passed to points().
#' @import rgdal
#' @export
add.map.points = function(map,
lon,
lat,
col = 'black',
cex = 1,
pch = 16, ...){
if (length(lon) != length(lat)) { stop('add.map.points: lengths of lon/lat are not the same.') }
xy = rgdal::project(cbind(lon, lat), proj = map$p)
points(xy[,1], xy[,2], col = col, cex = cex, pch = pch, ...)
}
#' @title Add Map Line
#' @author Thomas Bryce Kelly
#' @param lon a set of longitudes to draw a line between
#' @param lat a set of latitudes to draw a lien between
#' @param col the color of the line to be drawn
#' @param lty the type of line to be drawn
#' @param lty the type of line to be drawn
#' @param lwd the width of the line to be drawn
#' @export
add.map.line = function(map,
lon,
lat,
col = 'black',
lty = 1,
lwd = 1,
greatCircle = T,
N = 1e3) {
line = get.map.line(map = map, lon = lon, lat = lat, greatCircle = greatCircle, N = N)
## Plot
lines(line[,1], line[,2], col = col, lty = lty, lwd = lwd)
}
#' @title Add Map Arrows
#' @author Thomas Bryce Kelly
#' @description Add station points to a map
#' @param ... optional arguments passed to shape::Arrows().
#' @import rgdal shape
#' @export
add.map.arrows = function(map,
lon,
lat,
lon2,
lat2,
col = 'black',
cex = 1,
lwd = 1,
...){
if (length(lon) != length(lat) | length(lon) != length(lat2) | length(lon2) != length(lat2)) { stop('add.map.arrows: lengths of lon/lat are not the same.') }
xy = rgdal::project(cbind(lon, lat), proj = map$p)
xy2 = rgdal::project(cbind(lon2, lat2), proj = map$p)
shape::Arrows(x0 = xy[,1],
y0 = xy[,2],
x1 = xy2[,1],
y1 = xy2[,2],
col = col,
lwd = lwd,
arr.adj = 1,
arr.length = 0.4 * cex,
arr.width = 0.4 * cex,
...)
}
#' @title Add Map Text
#' @author Thomas Bryce Kelly
#' @author Laura Whitmore
#' @import rgdal
#' @param lon longitude position of the text or a vector of positions
#' @param lat latitude position of the text or a vector of positions
#' @param text a string or vector of strings for the text to be written
#' @param col text color
#' @param cex size of the text to be written
#' @param adj one or two values in [0, 1] which specify the x (and optionally y) adjustment of the labels.
#' @param pos a position specifier for the text. If specified this overrides any adj value given. Values of 1, 2, 3 and 4, respectively indicate positions below, to the left of, above and to the right of the specified coordinates.
#' @export
add.map.text = function(map,
lon,
lat,
text,
col = 'black',
cex = 1,
adj = NULL,
pos = NULL,
...){
xy = rgdal::project(cbind(lon, lat), proj = map$p)
text(xy[,1], xy[,2], text, col = col, cex = cex, adj = adj, pos = pos, ...)
}
#' @title Add Map Polygon
#' @author Laura Whitmore
#' @import rgdal
#' @export
add.map.polygon = function(map,
lon,
lat,
col = "#00000020",
lty = 1,
lwd = 1,
border = NULL,
density = NULL,
angle = 45,
fillOddEven = FALSE) {
## Project and plot normally
xy = rgdal::project(cbind(lon, lat), proj = map$p)
polygon(xy[,1], xy[,2], col = col, lty = lty, lwd = lwd, border = border, density = density, angle = angle, fillOddEven = fillOddEven)
}
#' @title Add Axis Label to Map
#' @author Thomas Bryce Kelly
#' @export
add.map.axis.label = function(map, temp, label, sides) {
usr = par('usr')
dx = 2e-2 * (usr[2] - usr[1])
dy = 2e-2 * (usr[4] - usr[3])
## bottom
if (1 %in% sides & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] > usr[3]) & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] <= usr[3])) {
if (min(abs(temp[,2] - usr[3])) < dx) {
## Check angle
k = which.min(abs(temp[-nrow(temp),2] - usr[3]))
dd = abs(diff(temp)[k,])
if (dd[1]/dd[2] < dy / dx) {
axis(1, at = temp[which.min((temp[,2] - usr[3])^2),1], labels = label)
}
}
}
## top
if (3 %in% sides & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] < usr[4]) & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] >= usr[4])) {
if (min(abs(temp[,2] - usr[4])) < dx) {
## Check angle
k = which.min(abs(temp[-nrow(temp),2] - usr[4]))
dd = abs(diff(temp)[k,])
if (dd[1]/dd[2] < dy / dx) {
axis(3, at = temp[which.min((temp[,2] - usr[4])^2),1], labels = label)
}
}
}
## Left
if (2 %in% sides & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] > usr[1]) & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] <= usr[1])) {
if (min(abs(temp[,1] - usr[1])) < dy) {
## Check angle
k = which.min(abs(temp[-nrow(temp),1] - usr[1]))
dd = abs(diff(temp)[k,])
if (dd[2]/dd[1] < dx / dy) {
axis(2, at = temp[which.min((temp[,1] - usr[1])^2),2], labels = label, las = 2)
}
}
}
## Right
if (4 %in% sides & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] < usr[2]) & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] >= usr[2])) {
if (min(abs(temp[,1] - usr[2])) < dy) {
## Check angle
k = which.min((temp[-nrow(temp),1] - usr[2])^2)
dd = abs(diff(temp)[k,])
if (dd[2]/dd[1] < dx / dy) {
axis(4, at = temp[which.min((temp[,1] - usr[2])^2),2], labels = label, las = 2)
}
}
}
}
#' @title Add Axis to Map
#' @author Thomas Bryce Kelly
#' @export
add.map.axis = function(map, lons = NULL, lats = NULL, sides = c(1:4), N = 100) {
if (is.null(lons)) { lons = seq(-180, 180, by = map$grid$dlon) }
if (is.null(lats)) { lats = seq(-90, 90, by = map$grid$dlat) }
## Add lon axes
for (ll in lons) {
temp = get.map.line(map, lon = rep(ll, 2), lat = c(-90,90), N = N)
#temp = temp[is.finite(temp[,1]) & is.finite(temp[,2]),]
s = 'E'
if (ll < 0 | ll %% 360 >= 180) { s = 'W'; ll = -ll }
add.map.axis.label(map, temp, label = paste0(ll, s), sides = sides)
}
## Add lats
for (ll in lats) {
temp = get.map.line(map, lon = c(-180, 180), lat = rep(ll, 2), N = N)
#temp = temp[is.finite(temp[,1]) & is.finite(temp[,2]),]
s = 'N'
if (ll < 0) { s = 'S'; ll = -ll }
add.map.axis.label(map, temp, label = paste0(ll, s), sides = sides)
}
}
#' @title Add Grid to Map
#' @author Thomas Bryce Kelly
#' @export
add.map.grid = function(map, lons = NULL, lats = NULL, col = 'darkgrey', lty = 1, lwd = 1) {
if (is.null(lons)) { lons = seq(-180, 180, by = map$grid$dlon) }
if (is.null(lats)) { lats = seq(-90, 90, by = map$grid$dlat) }
## Add lon grid
for (ll in lons) {
add.map.line(map, lon = rep(ll, 2), lat = c(-90,90), lwd = lwd, lty = lty, col = col)
}
## Add lats
for (ll in lats) {
add.map.line(map, lon = c(-180, 180), lat = rep(ll, 2), lwd = lwd, lty = lty, col = col)
}
}
#########################################3
############ Compound Functions ##########
##########################################
#' @title Add Map Scalebar
#' @author Thomas Bryce Kelly
#' @export
add.map.scale = function(pos = 1,
scale = NULL,
verbose = T,
col = 'black',
cex = 1) {
usr = par('usr')
if (is.null(scale)) {
scale = (usr[2] - usr[1]) * 0.001 #km
scale = 10^(round(log10(scale)-1))
}
sign = 1
if (pos == 1) {
x.origin = 0.95 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.05 * (usr[4] - usr[3]) + usr[3]
} else if (pos == 2) {
x.origin = 0.05 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.05 * (usr[4] - usr[3]) + usr[3]
sign = -1
} else if (pos == 3) {
x.origin = 0.05 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.925 * (usr[4] - usr[3]) + usr[3]
sign = -1
} else if (pos == 4) {
x.origin = 0.95 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.925 * (usr[4] - usr[3]) + usr[3]
} else {
stop ('pos should be 1, 2, 3, or 4!')
}
lines(x = c(x.origin, x.origin - sign*scale*1e3), y = rep(y.origin,2), lwd = 3*cex, col = col)
points(x = c(x.origin, x.origin - sign*scale*1e3), y = rep(y.origin,2), pch = '|', cex = cex, col = col)
text(x.origin - sign*scale*1e3*0.5, y = y.origin, pos = 3, paste0(scale, ' km'), col = col, cex = cex)
}
#' @title Add Quiver Lines
#' @export
#' @author Thomas Bryce Kelly
add.map.quiver = function(map,
lon,
lat,
u,
v,
zscale = 1,
cex = 1,
col = 'black',
lwd = 1,
show.scale = T,
verbose = T,
...) {
if (length(lon) != length(lat) | length(lon) != length(u) | length(lon) != length(v)) {stop('add.map.quiver: length of lat, lon, u and v must be the same.')}
lon = as.numeric(lon)
lat = as.numeric(lat)
u = as.numeric(u)
v = as.numeric(v)
## Filter our NAs
k = !is.na(u) & !is.na(v)
add.map.arrows(map, lon[k], lat[k], lon[k] + u[k] * zscale, lat[k] + v[k] * zscale, col = col, cex = cex / 2, lwd = lwd, ...)
if (show.scale) {
usr = par('usr')
dx = median(diff(unique(lon)))
dy = median(diff(unique(lon)))
x.origin = usr[2] - (usr[2] - usr[1])/10
y.origin = usr[3] + (usr[4] - usr[3])/10
lines(x = c(x.origin, x.origin - sign * scale * 1e3), y = rep(y.origin,2), lwd = 3 * cex, col = col)
points(x = c(x.origin, x.origin - sign * scale * 1e3), y = rep(y.origin,2), pch = '|', cex = cex, col = col)
text(x.origin - sign * scale * 1e3 * 0.5, y = y.origin, pos = 3, paste0(scale, ' km'), col = col, cex = cex)
}
}
#' @title Add Map Contours
#' @export
#' @author Thomas Bryce Kelly
add.map.contour = function(map,
lon,
lat,
z,
trim = T,
col = 'black',
levels = NULL,
zlim = NULL,
n = 3,
labels = TRUE,
lty = 1,
lwd = 1,
verbose = T) {
if (verbose) { message('Add Map Contours:')}
lon = lon %% 360
map$lon.min = map$lon.min %% 360
map$lon.max = map$lon.max %% 360
lon = as.array(lon)
lat = as.array(lat)
z = as.array(z)
if (is.null(zlim)) {
zlim = range(z, na.rm = T)
}
## Set default values
if (is.null(levels)) {
levels = pretty.default(zlim, n = n + 2)[-c(1, n+2)]
if (verbose) { message(' Levels: ', paste(levels, collapse = ', '))}
}
n = length(levels)
if (length(col) != n) { col = rep(col, n) }
if (length(lty) != n) { lty = rep(lty, n) }
if (length(lwd) != n) { lwd = rep(lwd, n) }
## make sure everything is a full grid
if(is.na(dim(lon)[2]) & is.na(dim(lat)[2])) {
if (length(z) == length(lon) * length(lat)) {
z = array(z, dim = c(length(lon), length(lat)))
lon = matrix(lon, nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(lat, nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
} else {
z = array(z, dim = c(length(unique(lon)), length(unique(lat))))
lon = matrix(unique(lon), nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(unique(lat), nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
}
}
## Trim
if (trim) {
nz = length(z)
if (verbose) { message(' Starting domain trimming... ', appendLF = F)}
corners = expand.grid(lon = c(par('usr')[1], par('usr')[2]),
lat = c(par('usr')[3], par('usr')[4]))
corners = rgdal::project(cbind(corners$lon, corners$lat), proj = map$p, inv = T)
field = expand.grid(lon = seq(par('usr')[1], par('usr')[2], length.out = 360),
lat = seq(par('usr')[3], par('usr')[4], length.out = 10))
field = rgdal::project(cbind(field$lon, field$lat), proj = map$p, inv = T)
field[,1] = field[,1] %% 360
field.lon = range(field[,1], na.rm = T)
field.lat = range(field[,2], na.rm = T)
## Trim longitude
if (corners[1,1] > corners[2,1]) { ## antimeridian
if (verbose) { message(' antimeridian... ', appendLF = F)}
k = apply(lon, 1, function(x) {any(x >= field.lon[1] & x <= field.lon[2])})
if (length(k) > 0) {
z = z[k,]
lon = lon[k,]
lat = lat[k,]
}
} else {
if (verbose) { message(' longitude... ', appendLF = F)}
k = apply(lon, 1, function(x) {any(x <= field.lon[2] & x >= field.lon[1])})
if (length(k) > 0) {
z = z[k,]
lon = lon[k,]
lat = lat[k,]
}
}
## Trim latitude
if (verbose) { message(' latitude... ', appendLF = F) }
k = apply(lat, 2, function(x) {any(x >= field.lat[1] & x <= field.lat[2])})
if (length(k) > 0) {
z = z[,k]
lon = lon[,k]
lat = lat[,k]
}
if (verbose) { message(' complete, n = ', length(z), ' (', 100*round(1 - n/nz, digits = 3), ' %)')}
}
## Calculate contors and plot
contour = contourLines(x = c(1:dim(z)[1]),
y = c(1:dim(z)[2]),
z = z,
levels = levels)
if (verbose) { message(' Starting plotting... ', appendLF = F) }
for (i in 1:length(contour)) {
n = which(contour[[i]]$level == levels)
if (length(contour[[i]]$x) > 5) {
xx = grid.interp(lon, contour[[i]]$x, contour[[i]]$y)
yy = grid.interp(lat, contour[[i]]$x, contour[[i]]$y)
add.map.line(map,
xx,
yy,
col = col[n],
lty = lty[n],
lwd = lwd[n],
greatCircle = F)
}
}
redraw.map(map)
if (verbose) { message(' done.') }
}
#' @title Add Map Layer
#' @description Add a image layer to the map!
#' @author Thomas Bryce Kelly
#' @param lon longitude of the data
#' @param lat latitude of the data
#' @param z the data
#' @param zlim the limits of the z-axis. A value of NULL indicates that zlim should equal the range of z values.
#' @param pal the name of a palette generating function. Should be a string.
#' @param col.na the color with which any NAs should be drawn with. A value of NA skips this step.
#' @param n the number of distinct colors to request from the palette function.
#' @param refinement the level of bilinear refinement to apply to the image grid
#' @export
add.map.layer = function(map,
lon,
lat,
z,
zlim = NULL,
pal = 'greyscale',
col.na = NA,
n = 255,
refine = 0,
trim = T,
col.low = '',
col.high = '',
rev = F,
indicate = T,
verbose = T) {
## Misc corrections
lon = lon %% 360
map$lon.min = map$lon.min %% 360
map$lon.max = map$lon.max %% 360
lon = as.array(lon)
lat = as.array(lat)
z = as.array(z)
if (length(dim(lon)) < 2) {
if (verbose) { message(' Establishing Grid...', appendLF = F) }
if (length(z) == 1.0 * length(lon) * length(lat)) {
z = array(z, dim = c(length(lon), length(lat)))
lon = matrix(lon, nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(lat, nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
} else { ## TODO add check here that the dimentions are actually correct!
z = array(z, dim = c(length(unique(lon)), length(unique(lat))))
lon = matrix(unique(lon), nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(unique(lat), nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
}
if (verbose) { message(' Done. ') }
}
nz = length(z)
## Trim
if (trim & length(z) > 100) {
if (verbose) { message(' Starting domain trimming... ', appendLF = F)}
usr = par('usr')
corners = expand.grid(lon = c(usr[1], usr[2]),
lat = c(usr[3], usr[4]))
corners = rgdal::project(cbind(corners$lon, corners$lat), proj = map$p, inv = T)
if (corners[1,1] > corners[2,1]) { ## antimeridian
if (verbose) { message(' antimeridian... ', appendLF = F)}
antimeridian = T
} else {
antimeridian = F
lon = lon %% 360
lon[lon > 180] = lon[lon > 180] - 360
}
field = expand.grid(lon = seq(usr[1], usr[2], length.out = 100),
lat = seq(usr[3], usr[4], length.out = 100))
field = rgdal::project(cbind(field$lon, field$lat), proj = map$p, inv = T)
field[,1] = field[,1] %% 360
if (!antimeridian & any(field[,1] > 180)) {
l = which(!is.na(field[,1]) & field[,1] > 180)
field[l,1] = field[l,1] - 360
}
field.lon = range(field[,1], na.rm = T) %% 360
field.lat = range(field[,2], na.rm = T)
## Trim longitude
if (field.lat[1] > -80 & field.lat[2] < 80) { ## only if a pole isn't visible!
if (verbose) { message(' longitude... ', appendLF = F)}
if (corners[1,1] > corners[2,1]) { ## antimeridian
k = apply(lon, 1, function(x) {any(x >= field.lon[1] & x <= field.lon[2])})
} else {
k = apply(lon, 1, function(x) {any(x <= field.lon[2] & x >= field.lon[1])})
}
if (sum(k) > 2) {
z = z[k,]
lon = lon[k,]
lat = lat[k,]
}
}
## Trim latitude
if (verbose) { message(' latitude... ', appendLF = F) }
k = apply(lat, 2, function(x) {any(x >= field.lat[1] & x <= field.lat[2])})
if (sum(k) > 2) {
z = z[,k]
lon = lon[,k]
lat = lat[,k]
}
if (verbose) { message(' complete, n = ', length(z), ' (', 100 - round(100 * length(z) / nz), '% trimmed)')}
}
if (refine > 0) {
if (!antimeridan) {
lon[lon > 180] = lon[lon>180] - 360
}
for (i in 1:refine) {
temp = grid.refinement(lon, lat, z)
lon = temp$x
lat = temp$y
z = temp$z
}
if (verbose) {message(' Refined grid by ', 2^refine,'x.')}
} else if (refine < 0) {
for (i in 1:abs(refine)) {
temp = grid.subsample(lon, lat, z)
lon = temp$x
lat = temp$y
z = temp$z
}
if (verbose) {message(' Refined grid by 1:', 2^-refine,'x.')}
}
if (is.null(zlim)) { zlim = range(pretty(as.numeric(z), na.rm = TRUE)) }
## Order: Not needed?
#l = order(rgdal::project(cbind(lon[,round(dim(lon)[2]/2)], lat[,round(dim(lon)[2]/2)]), p = map$p)[,1])
#l = order(lon[,round(dim(lon)[2]/2)])
#lon = lon[l,]
#lat = lat[l,]
#z = z[l,]
## Color scale
col = make.pal(z, pal = pal, rev = rev, n = n, min = zlim[1], max = zlim[2])
col = array(col, dim = dim(lon))
## Apply out of range colors if provided:
if (is.na(col.low)) { col[z < zlim[1]] = NA }
if (is.na(col.high)) { col[z > zlim[2]] = NA }
## Project and Plot
if (verbose) {message(' Projecting grid...')}
#vertex = calc.vertex(lon, lat)
#xy = rgdal::project(cbind(as.numeric(vertex$x), as.numeric(vertex$y)), p = map$p)
#xy = list(x = array(xy[,1], dim = dim(vertex$x)),
# y = array(xy[,2], dim = dim(vertex$x)))
xy = rgdal::project(cbind(as.numeric(lon), as.numeric(lat)), p = map$p)
xy = list(x = array(xy[,1], dim = dim(lon)),
y = array(xy[,2], dim = dim(lat)))
xy = calc.vertex(xy$x, xy$y)
if (verbose) { message(' Starting plotting... ', appendLF = F) }
for (i in 1:dim(col)[1]) {
for (j in 1:dim(col)[2]) {
if (!is.na(col[i,j])) {
polygon(x = c(xy$x[i,j], xy$x[i,j+1], xy$x[i+1,j+1], xy$x[i+1,j]),
y = c(xy$y[i,j], xy$y[i,j+1], xy$y[i+1,j+1], xy$y[i+1,j]),
col = col[i,j], border = NA)
}
}
}
if (verbose) { message(' complete.') }
## Extras
if (verbose) { message(' Final stats: \tN.low: ', sum(z < zlim[1]), '\tN.high: ', sum(z > zlim[2]), '\tN: ', length(z)) }
if (indicate) {
st = paste0('Data range: (', round(min(z, na.rm = TRUE), 3), ', ', round(max(z, na.rm = TRUE), 3),
') Z range: (', round(zlim[1], 3), ', ', round(zlim[2], 3), ')'
)
mtext(st, line = 0.25, adj = 1, cex = 0.7)
}
}
#' @export
add.map.eez = function(map, eez, names = NULL, col = 'black', lwd = 1, lty = 1, verbose = T) {
if (is.null(names)) {
names = names(eez)
if (verbose) {
message('No names provided, using all available: ', paste0(names, collapse = ', '))
}
}
for (n in names) {
l = which(names(eez) == n)
if (length(l) > 0) {
for (i in 1:length(eez[[n]])) {
add.map.line(map, eez[[n]][[i]]$lon, eez[[n]][[i]]$lat, col = col, lwd = lwd, lty = lty)
}
}
}
}
#' @title Add Map Bathymetry (shaded)
#' @author Thomas Bryce Kelly
#' @param map a map object returned by the function make.map()
#' @param bathy a bathymetric dataframe returned by get.map.bathy
#' @param pal the color palette used for coloring the shading. Should consist of a function name as a string
#' @param n the number of distinct colors to use
#' @param zlim the z-axis range that the colormap is projected upon
#' @param rev a boolean flag to flip the color axis
#' @param filled a flag to turn on color smoothing and filling from the OCE package. Works for some datasets but may cause issues.
#' @param col.low the color used when painting data that is below the zlim range. A value of NA causes the data not to be drawn, and a value of '' uses the lowest color value for out of range data.
#' @param col.high same as \emph{col.low} but for data that is out of range above the zlim.
#' @param refinement the level of bilinear refinement to apply to the image grid
#' @export
add.map.bathy = function(map,
bathy,
refine = 0,
trim = T,
zlim = NULL,
pal = 'greyscale',
n = 255,
rev = FALSE,
col.low = '',
col.high = '',
verbose = T) {
add.map.layer(map,
lon = bathy$Lon,
lat = bathy$Lat,
z = bathy$Z,
refine = refine,
pal = pal,
rev = rev,
zlim = zlim,
col.low = col.low,
col.high = col.high,
indicate = FALSE,
n = n,
trim = trim,
verbose = verbose)
redraw.map(map)
}
#' @title Redraw Map
#' @author Thomas Bryce Kelly
#' @param map a map object as returned by make.map()
#' @export
redraw.map = function(map) {
for (i in 1:length(map$coastline$projected)) {
polygon(map$coastline$projected[[i]][,1], map$coastline$projected[[i]][,2], col = map$land.col)
}
if (map$grid$draw.grid) {
for (i in map$grid$lon.axis$label) {
add.map.line(map, lon = c(i,i), lat = c(-80, 80), col = 'dark grey')
}
for (i in map$grid$lat.axis$label) {
add.map.line(map, lon = c(-180,180), lat = c(i,i), col = 'dark grey')
}
}
box()
}
tbrycekelly/TheSource documentation built on Nov. 7, 2023, 12:48 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.
Defines functions redraw.map add.map.bathy add.map.eez add.map.layer add.map.contour add.map.quiver add.map.scale add.map.grid add.map.axis add.map.axis.label add.map.polygon add.map.text add.map.arrows add.map.line add.map.points get.map.line add.map.coastline
Documented in add.map.arrows add.map.axis add.map.axis.label add.map.bathy add.map.coastline add.map.contour add.map.grid add.map.layer add.map.line add.map.points add.map.polygon add.map.quiver add.map.scale add.map.text get.map.line redraw.map
#' @title Make Map
#' @description Default settings are for Arctic, lowest resolution coastline.
#' Exaples of projection arguments include: `+proj=stere +lat_0=90`, `+proj=merc`, `+proj=aea +lat_1=30 +lat_2=45 +lon_0=-120`.
#' @param coast Should be the name of a coastline data object. A value of NULL sets the default cosatline to 'coastline2'.
#' @param lon.min The minimum longitude displayed on the map.
#' @param lon.max The maximum longitude on the map
#' @param lat.min The minimum latitude shown on the map
#' @param lat.max The maximum latitude shown on the map
#' @param p a projection string, such as those generated by 'make.proj()'
#' @param land.col A color string or value for the land polygon
#' @param grid Boolean, draw lat/lon grid?
#' @param dlon The spacing for the longitude grid (in degrees)
#' @param dlat The spacing for the latitude grid (in degrees)
#' @import oce
#' @import ocedata
#' @author Laura Whitmore
#' @author Thomas Bryce Kelly
#' @export
make.map = function (coast = NULL,
lon.min = -180,
lon.max = 180,
lat.min = -60,
lat.max = 80,
p = NULL,
land.col = 'lightgray',
draw.grid = T,
draw.axis = T,
dlon = 60,
dlat = 20,
verbose = T) {
## Apply Defaults
if (all(is.null(c(lon.min, lon.max, lat.min, lat.max)))) {
lon.min = -180
lon.max = 180
lat.min = -60
lat.max = 80
} else if (any(is.null(c(lon.min, lon.max, lat.min, lat.max)))) {
stop('All lon/lat bounds must be provided!')
}
if (is.null(coast)) { coast = 'coastline2' }
if (is.null(p)) {
p = make.proj(projection = 1,
lat = mean(c(lat.max, lat.min)),
lon = mean(c(lon.min, lon.max)))
}
## Determine lat/lon axis details
lons = seq(lon.min - 5*dlon, lon.max + 5*dlon, by = dlon)
lats = seq(lat.min - 5*dlat, lat.max + 5*dlon, by = dlat)
if (lon.min < 0) {
lons = lons[lons >= -180 & lons <= 180]
} else {
lons = lons[lons >= 0 & lons <= 360]
}
lats = lats[lats >= -90 & lats <= 90]
lims = expand.grid(lon = lons, lats = lats)
lon.axis = rgdal::project(cbind(lons, rep(lat.min, length(lons))), proj = p)[,1]
lat.axis = rgdal::project(cbind(rep(lon.min, length(lats)), lats), proj = p)[,2]
lims = rgdal::project(cbind(lims$lon, lims$lat), proj = p)
## Start plotting
# Calculate projections of plotting area for plot limits
field = expand.grid(lon = seq(lon.min, lon.max, length.out = 10),
lat = seq(lat.min, lat.max, length.out = 10))
field = rgdal::project(cbind(field$lon, field$lat), proj = p)
## make plot
plot(NULL,
NULL,
xlim = range(field[,1], na.rm = T),
ylim = range(field[,2], na.rm = T),
xaxt = 'n', yaxt = 'n',
xlab = '', ylab = '',
xaxs = 'i', yaxs = 'i')
coast = add.map.coastline(coast, p = p, land.col = land.col)
## Setup map object
map = list(coastline = coast,
lon = (lon.min + lon.max)/2,
lat = (lat.min + lat.max)/2,
scale = NA,
p = p,
land.col = land.col,
grid = list(
draw.grid = draw.grid,
draw.axis = draw.axis,
dlon = dlon,
dlat = dlat
),
meta = list(
Source.version = packageVersion('TheSource'),
R.version = R.version.string,
Sys.time = Sys.time()
))
if (draw.grid) { add.map.grid(map) }
if (draw.axis) { add.map.axis(map, sides = c(1,2)) }
box()
map
}
#' @title Add Coastline to Map
#' @author Thomas Bryce Kelly
#' @description Adds a coastline list to a map given a projection and land color.
#' @export
add.map.coastline = function(coast, p, land.col, field.inv = NULL) {
## helper function
break.polygon = function(projected.coast) {
for (i in 1:length(projected.coast)) {
k = c(which(diff(projected.coast[[i]][,1])^2 + diff(projected.coast[[i]][,2])^2 > diff(par('usr')[1:2])^2 / 10), length(projected.coast[[i]][,1]))
if (length(k) > 1) {
for (j in 2:length(k)) {
## split off parts of polygons with more than 5 points
if (k[j] - k[j-1] > 5) {
projected.coast[[length(projected.coast) + 1]] = projected.coast[[i]][(k[j-1]+1):k[j],]
projected.coast[[i]][(k[j-1]+1):k[j],] = NA ## Remove newly split areas
}
}
## Trim original down to just first likely contiguous area
projected.coast[[i]] = array(projected.coast[[i]][1:k[1],], dim = c(k[1], 2))
}
}
## Return
projected.coast
}
# Load data if necessary
if (typeof(coast) != "character") {
coastline = coast
} else {
do.call('data', list(coast))
coastline = eval(parse(text = coast))
}
## First pass at trimming domain
if (!is.null(field.inv)) {
keep = rep(T, length(coastline$data))
dx = median(diff(field.inv[,1])) * 2
dy = median(diff(field.inv[,2])) * 2
for (i in 1:length(coastline$data)) {
if (all(coastline$data[[i]][,1] + dx < min(field.inv[,1], na.rm = T)) | all(coastline$data[[i]][,1] - dx > max(field.inv[,1], na.rm = T)) |
all(coastline$data[[i]][,2] + dy < min(field.inv[,2], na.rm = T)) | all(coastline$data[[i]][,2] - dy > max(field.inv[,2]), na.rm = T)) {
keep[i] = F
}
}
if (sum(keep) == 0) {
return(list(data = NA, meta = coastline$meta, projected = NA))
}
coastline$data = coastline$data[keep]
coastline$meta$trim1 = sum(keep)
}
## Project coastline (takes a while!)
projected.coast = lapply(coastline$data, function(x) {
rgdal::project(as.matrix(x)+1e-16, proj = p)
})
## Second round of trimming
keep = rep(T, length(projected.coast))
usr = par('usr')
for (i in 1:length(projected.coast)) {
if (all(projected.coast[[i]][,1] < usr[1]) | all(projected.coast[[i]][,1] > usr[2]) |
all(projected.coast[[i]][,2] < usr[3]) | all(projected.coast[[i]][,2] > usr[4])) {
keep[i] = F
}
}
if (sum(keep) == 0) {
return(list(data = NA, meta = coastline$meta, projected = NA))
}
coastline$data = coastline$data[keep]
projected.coast = projected.coast[keep]
projected.coast = break.polygon(projected.coast)
coastline$meta$trim2 = sum(keep)
## Add coastlines (takes a while!)
for (i in 1:length(projected.coast)) {
polygon(projected.coast[[i]][,1], projected.coast[[i]][,2], col = land.col)
}
list(data = coastline$data, meta = coastline$meta, projected = projected.coast)
}
#' @title Make Map
#' @description Default settings are for Arctic, lowest resolution coastline.
#' Exaples of projection arguments include: `+proj=stere +lat_0=90`, `+proj=merc`, `+proj=aea +lat_1=30 +lat_2=45 +lon_0=-120`.
#' @param coast Should be the name of a coastline data object. A value of NULL sets the default cosatline to 'coastline2'.
#' @param lon.min The minimum longitude displayed on the map.
#' @param lon.max The maximum longitude on the map
#' @param lat.min The minimum latitude shown on the map
#' @param lat.max The maximum latitude shown on the map
#' @param p a projection string, such as those generated by 'make.proj()'
#' @param land.col A color string or value for the land polygon
#' @param grid Boolean, draw lat/lon grid?
#' @param dlon The spacing for the longitude grid (in degrees)
#' @param dlat The spacing for the latitude grid (in degrees)
#' @import oce
#' @import ocedata
#' @author Laura Whitmore
#' @author Thomas Bryce Kelly
#' @export
make.map2 = function (coast = NULL,
lon = 0,
lat = 0,
scale = 1000,
p = NULL,
land.col = 'lightgray',
draw.grid = T,
draw.axis = T,
dlon = 60,
dlat = 20,
verbose = T) {
## Apply Defaults
if (is.null(coast)) { coast = 'coastline2' }
if (is.null(p)) {
p = make.proj(projection = 1,
lat = lat,
lon = lon)
}
center.point = rgdal::project(cbind(lon, lat), proj = p)
x.y.ratio = par()$pin[1] / par()$pin[2]
## Start plotting
# Calculate projections of plotting area for plot limits
field = expand.grid(lon = seq(center.point[,1] - scale * x.y.ratio * 1e3, center.point[,1] + scale * x.y.ratio * 1e3, length.out = 10),
lat = seq(center.point[,2] - scale * 1e3, center.point[,2] + scale * 1e3, length.out = 10))
field.inv = rgdal::project(cbind(field$lon, field$lat), proj = p, inv = T)
## make plot
plot(NULL,
NULL,
xlim = range(field[,1], na.rm = T),
ylim = range(field[,2], na.rm = T),
xaxt = 'n', yaxt = 'n',
xlab = '', ylab = '',
xaxs = 'i', yaxs = 'i')
coast = add.map.coastline(coast, p = p, land.col = land.col, field.inv = field.inv)
## Setup map object
map = list(coastline = coast,
lon = lon,
lat = lat,
scale = scale,
p = p,
land.col = land.col,
grid = list(
draw.grid = draw.grid,
draw.axis = draw.axis,
dlon = dlon,
dlat = dlat
),
meta = list(
Source.version = packageVersion('TheSource'),
R.version = R.version.string,
Sys.time = Sys.time()
))
if (draw.grid) { add.map.grid(map) }
if (draw.axis) { add.map.axis(map, sides = c(1,2)) }
box()
map
}
#' @title Get Map Line
#' @author Thomas Bryce Kelly
#' @import rgdal
#' @param lon a set of longitudes to draw a line between
#' @param lat a set of latitudes to draw a lien between
#' @param col the color of the line to be drawn
#' @param lty the type of line to be drawn
#' @param lty the type of line to be drawn
#' @param lwd the width of the line to be drawn
#' @export
get.map.line = function(map,
lon,
lat,
greatCircle = F,
N = 1e3) {
if (length(lon) != length(lat)) {stop('get.map.line: length of lon/lat are not the same.')}
if (greatCircle) {
lons = approx(1:length(lon), lon, xout = seq(1, length(lon), length.out = max(length(lon), 1e4)))$y
lats = approx(1:length(lat), lat, xout = seq(1, length(lat), length.out = max(length(lon), 1e4)))$y
} else {
lons = lon
lats = lat
}
## project
rgdal::project(cbind(lons, lats), proj = map$p)
}
#' @title Add Map Points
#' @author Laura Whitmore
#' @description Add station points to a map
#' @param lon the longitudes of the points to be drawn
#' @param lat the latitude of the points to be drawn
#' @param stn.lon (depreciated) the longitudes of the points to be drawn
#' @param stn.lat (depreciated) the latitude of the points to be drawn
#' @param col the color of the points
#' @param cex the size of point to be drawn
#' @param pch the point character to be used
#' @param ... optional arguments passed to points().
#' @import rgdal
#' @export
add.map.points = function(map,
lon,
lat,
col = 'black',
cex = 1,
pch = 16, ...){
if (length(lon) != length(lat)) { stop('add.map.points: lengths of lon/lat are not the same.') }
xy = rgdal::project(cbind(lon, lat), proj = map$p)
points(xy[,1], xy[,2], col = col, cex = cex, pch = pch, ...)
}
#' @title Add Map Line
#' @author Thomas Bryce Kelly
#' @param lon a set of longitudes to draw a line between
#' @param lat a set of latitudes to draw a lien between
#' @param col the color of the line to be drawn
#' @param lty the type of line to be drawn
#' @param lty the type of line to be drawn
#' @param lwd the width of the line to be drawn
#' @export
add.map.line = function(map,
lon,
lat,
col = 'black',
lty = 1,
lwd = 1,
greatCircle = T,
N = 1e3) {
line = get.map.line(map = map, lon = lon, lat = lat, greatCircle = greatCircle, N = N)
## Plot
lines(line[,1], line[,2], col = col, lty = lty, lwd = lwd)
}
#' @title Add Map Arrows
#' @author Thomas Bryce Kelly
#' @description Add station points to a map
#' @param ... optional arguments passed to shape::Arrows().
#' @import rgdal shape
#' @export
add.map.arrows = function(map,
lon,
lat,
lon2,
lat2,
col = 'black',
cex = 1,
lwd = 1,
...){
if (length(lon) != length(lat) | length(lon) != length(lat2) | length(lon2) != length(lat2)) { stop('add.map.arrows: lengths of lon/lat are not the same.') }
xy = rgdal::project(cbind(lon, lat), proj = map$p)
xy2 = rgdal::project(cbind(lon2, lat2), proj = map$p)
shape::Arrows(x0 = xy[,1],
y0 = xy[,2],
x1 = xy2[,1],
y1 = xy2[,2],
col = col,
lwd = lwd,
arr.adj = 1,
arr.length = 0.4 * cex,
arr.width = 0.4 * cex,
...)
}
#' @title Add Map Text
#' @author Thomas Bryce Kelly
#' @author Laura Whitmore
#' @import rgdal
#' @param lon longitude position of the text or a vector of positions
#' @param lat latitude position of the text or a vector of positions
#' @param text a string or vector of strings for the text to be written
#' @param col text color
#' @param cex size of the text to be written
#' @param adj one or two values in [0, 1] which specify the x (and optionally y) adjustment of the labels.
#' @param pos a position specifier for the text. If specified this overrides any adj value given. Values of 1, 2, 3 and 4, respectively indicate positions below, to the left of, above and to the right of the specified coordinates.
#' @export
add.map.text = function(map,
lon,
lat,
text,
col = 'black',
cex = 1,
adj = NULL,
pos = NULL,
...){
xy = rgdal::project(cbind(lon, lat), proj = map$p)
text(xy[,1], xy[,2], text, col = col, cex = cex, adj = adj, pos = pos, ...)
}
#' @title Add Map Polygon
#' @author Laura Whitmore
#' @import rgdal
#' @export
add.map.polygon = function(map,
lon,
lat,
col = "#00000020",
lty = 1,
lwd = 1,
border = NULL,
density = NULL,
angle = 45,
fillOddEven = FALSE) {
## Project and plot normally
xy = rgdal::project(cbind(lon, lat), proj = map$p)
polygon(xy[,1], xy[,2], col = col, lty = lty, lwd = lwd, border = border, density = density, angle = angle, fillOddEven = fillOddEven)
}
#' @title Add Axis Label to Map
#' @author Thomas Bryce Kelly
#' @export
add.map.axis.label = function(map, temp, label, sides) {
usr = par('usr')
dx = 2e-2 * (usr[2] - usr[1])
dy = 2e-2 * (usr[4] - usr[3])
## bottom
if (1 %in% sides & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] > usr[3]) & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] <= usr[3])) {
if (min(abs(temp[,2] - usr[3])) < dx) {
## Check angle
k = which.min(abs(temp[-nrow(temp),2] - usr[3]))
dd = abs(diff(temp)[k,])
if (dd[1]/dd[2] < dy / dx) {
axis(1, at = temp[which.min((temp[,2] - usr[3])^2),1], labels = label)
}
}
}
## top
if (3 %in% sides & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] < usr[4]) & any(temp[,1] > usr[1] & temp[,1] < usr[2] & temp[,2] >= usr[4])) {
if (min(abs(temp[,2] - usr[4])) < dx) {
## Check angle
k = which.min(abs(temp[-nrow(temp),2] - usr[4]))
dd = abs(diff(temp)[k,])
if (dd[1]/dd[2] < dy / dx) {
axis(3, at = temp[which.min((temp[,2] - usr[4])^2),1], labels = label)
}
}
}
## Left
if (2 %in% sides & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] > usr[1]) & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] <= usr[1])) {
if (min(abs(temp[,1] - usr[1])) < dy) {
## Check angle
k = which.min(abs(temp[-nrow(temp),1] - usr[1]))
dd = abs(diff(temp)[k,])
if (dd[2]/dd[1] < dx / dy) {
axis(2, at = temp[which.min((temp[,1] - usr[1])^2),2], labels = label, las = 2)
}
}
}
## Right
if (4 %in% sides & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] < usr[2]) & any(temp[,2] > usr[3] & temp[,2] < usr[4] & temp[,1] >= usr[2])) {
if (min(abs(temp[,1] - usr[2])) < dy) {
## Check angle
k = which.min((temp[-nrow(temp),1] - usr[2])^2)
dd = abs(diff(temp)[k,])
if (dd[2]/dd[1] < dx / dy) {
axis(4, at = temp[which.min((temp[,1] - usr[2])^2),2], labels = label, las = 2)
}
}
}
}
#' @title Add Axis to Map
#' @author Thomas Bryce Kelly
#' @export
add.map.axis = function(map, lons = NULL, lats = NULL, sides = c(1:4), N = 100) {
if (is.null(lons)) { lons = seq(-180, 180, by = map$grid$dlon) }
if (is.null(lats)) { lats = seq(-90, 90, by = map$grid$dlat) }
## Add lon axes
for (ll in lons) {
temp = get.map.line(map, lon = rep(ll, 2), lat = c(-90,90), N = N)
#temp = temp[is.finite(temp[,1]) & is.finite(temp[,2]),]
s = 'E'
if (ll < 0 | ll %% 360 >= 180) { s = 'W'; ll = -ll }
add.map.axis.label(map, temp, label = paste0(ll, s), sides = sides)
}
## Add lats
for (ll in lats) {
temp = get.map.line(map, lon = c(-180, 180), lat = rep(ll, 2), N = N)
#temp = temp[is.finite(temp[,1]) & is.finite(temp[,2]),]
s = 'N'
if (ll < 0) { s = 'S'; ll = -ll }
add.map.axis.label(map, temp, label = paste0(ll, s), sides = sides)
}
}
#' @title Add Grid to Map
#' @author Thomas Bryce Kelly
#' @export
add.map.grid = function(map, lons = NULL, lats = NULL, col = 'darkgrey', lty = 1, lwd = 1) {
if (is.null(lons)) { lons = seq(-180, 180, by = map$grid$dlon) }
if (is.null(lats)) { lats = seq(-90, 90, by = map$grid$dlat) }
## Add lon grid
for (ll in lons) {
add.map.line(map, lon = rep(ll, 2), lat = c(-90,90), lwd = lwd, lty = lty, col = col)
}
## Add lats
for (ll in lats) {
add.map.line(map, lon = c(-180, 180), lat = rep(ll, 2), lwd = lwd, lty = lty, col = col)
}
}
#########################################3
############ Compound Functions ##########
##########################################
#' @title Add Map Scalebar
#' @author Thomas Bryce Kelly
#' @export
add.map.scale = function(pos = 1,
scale = NULL,
verbose = T,
col = 'black',
cex = 1) {
usr = par('usr')
if (is.null(scale)) {
scale = (usr[2] - usr[1]) * 0.001 #km
scale = 10^(round(log10(scale)-1))
}
sign = 1
if (pos == 1) {
x.origin = 0.95 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.05 * (usr[4] - usr[3]) + usr[3]
} else if (pos == 2) {
x.origin = 0.05 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.05 * (usr[4] - usr[3]) + usr[3]
sign = -1
} else if (pos == 3) {
x.origin = 0.05 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.925 * (usr[4] - usr[3]) + usr[3]
sign = -1
} else if (pos == 4) {
x.origin = 0.95 * (usr[2] - usr[1]) + usr[1]
y.origin = 0.925 * (usr[4] - usr[3]) + usr[3]
} else {
stop ('pos should be 1, 2, 3, or 4!')
}
lines(x = c(x.origin, x.origin - sign*scale*1e3), y = rep(y.origin,2), lwd = 3*cex, col = col)
points(x = c(x.origin, x.origin - sign*scale*1e3), y = rep(y.origin,2), pch = '|', cex = cex, col = col)
text(x.origin - sign*scale*1e3*0.5, y = y.origin, pos = 3, paste0(scale, ' km'), col = col, cex = cex)
}
#' @title Add Quiver Lines
#' @export
#' @author Thomas Bryce Kelly
add.map.quiver = function(map,
lon,
lat,
u,
v,
zscale = 1,
cex = 1,
col = 'black',
lwd = 1,
show.scale = T,
verbose = T,
...) {
if (length(lon) != length(lat) | length(lon) != length(u) | length(lon) != length(v)) {stop('add.map.quiver: length of lat, lon, u and v must be the same.')}
lon = as.numeric(lon)
lat = as.numeric(lat)
u = as.numeric(u)
v = as.numeric(v)
## Filter our NAs
k = !is.na(u) & !is.na(v)
add.map.arrows(map, lon[k], lat[k], lon[k] + u[k] * zscale, lat[k] + v[k] * zscale, col = col, cex = cex / 2, lwd = lwd, ...)
if (show.scale) {
usr = par('usr')
dx = median(diff(unique(lon)))
dy = median(diff(unique(lon)))
x.origin = usr[2] - (usr[2] - usr[1])/10
y.origin = usr[3] + (usr[4] - usr[3])/10
lines(x = c(x.origin, x.origin - sign * scale * 1e3), y = rep(y.origin,2), lwd = 3 * cex, col = col)
points(x = c(x.origin, x.origin - sign * scale * 1e3), y = rep(y.origin,2), pch = '|', cex = cex, col = col)
text(x.origin - sign * scale * 1e3 * 0.5, y = y.origin, pos = 3, paste0(scale, ' km'), col = col, cex = cex)
}
}
#' @title Add Map Contours
#' @export
#' @author Thomas Bryce Kelly
add.map.contour = function(map,
lon,
lat,
z,
trim = T,
col = 'black',
levels = NULL,
zlim = NULL,
n = 3,
labels = TRUE,
lty = 1,
lwd = 1,
verbose = T) {
if (verbose) { message('Add Map Contours:')}
lon = lon %% 360
map$lon.min = map$lon.min %% 360
map$lon.max = map$lon.max %% 360
lon = as.array(lon)
lat = as.array(lat)
z = as.array(z)
if (is.null(zlim)) {
zlim = range(z, na.rm = T)
}
## Set default values
if (is.null(levels)) {
levels = pretty.default(zlim, n = n + 2)[-c(1, n+2)]
if (verbose) { message(' Levels: ', paste(levels, collapse = ', '))}
}
n = length(levels)
if (length(col) != n) { col = rep(col, n) }
if (length(lty) != n) { lty = rep(lty, n) }
if (length(lwd) != n) { lwd = rep(lwd, n) }
## make sure everything is a full grid
if(is.na(dim(lon)[2]) & is.na(dim(lat)[2])) {
if (length(z) == length(lon) * length(lat)) {
z = array(z, dim = c(length(lon), length(lat)))
lon = matrix(lon, nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(lat, nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
} else {
z = array(z, dim = c(length(unique(lon)), length(unique(lat))))
lon = matrix(unique(lon), nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(unique(lat), nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
}
}
## Trim
if (trim) {
nz = length(z)
if (verbose) { message(' Starting domain trimming... ', appendLF = F)}
corners = expand.grid(lon = c(par('usr')[1], par('usr')[2]),
lat = c(par('usr')[3], par('usr')[4]))
corners = rgdal::project(cbind(corners$lon, corners$lat), proj = map$p, inv = T)
field = expand.grid(lon = seq(par('usr')[1], par('usr')[2], length.out = 360),
lat = seq(par('usr')[3], par('usr')[4], length.out = 10))
field = rgdal::project(cbind(field$lon, field$lat), proj = map$p, inv = T)
field[,1] = field[,1] %% 360
field.lon = range(field[,1], na.rm = T)
field.lat = range(field[,2], na.rm = T)
## Trim longitude
if (corners[1,1] > corners[2,1]) { ## antimeridian
if (verbose) { message(' antimeridian... ', appendLF = F)}
k = apply(lon, 1, function(x) {any(x >= field.lon[1] & x <= field.lon[2])})
if (length(k) > 0) {
z = z[k,]
lon = lon[k,]
lat = lat[k,]
}
} else {
if (verbose) { message(' longitude... ', appendLF = F)}
k = apply(lon, 1, function(x) {any(x <= field.lon[2] & x >= field.lon[1])})
if (length(k) > 0) {
z = z[k,]
lon = lon[k,]
lat = lat[k,]
}
}
## Trim latitude
if (verbose) { message(' latitude... ', appendLF = F) }
k = apply(lat, 2, function(x) {any(x >= field.lat[1] & x <= field.lat[2])})
if (length(k) > 0) {
z = z[,k]
lon = lon[,k]
lat = lat[,k]
}
if (verbose) { message(' complete, n = ', length(z), ' (', 100*round(1 - n/nz, digits = 3), ' %)')}
}
## Calculate contors and plot
contour = contourLines(x = c(1:dim(z)[1]),
y = c(1:dim(z)[2]),
z = z,
levels = levels)
if (verbose) { message(' Starting plotting... ', appendLF = F) }
for (i in 1:length(contour)) {
n = which(contour[[i]]$level == levels)
if (length(contour[[i]]$x) > 5) {
xx = grid.interp(lon, contour[[i]]$x, contour[[i]]$y)
yy = grid.interp(lat, contour[[i]]$x, contour[[i]]$y)
add.map.line(map,
xx,
yy,
col = col[n],
lty = lty[n],
lwd = lwd[n],
greatCircle = F)
}
}
redraw.map(map)
if (verbose) { message(' done.') }
}
#' @title Add Map Layer
#' @description Add a image layer to the map!
#' @author Thomas Bryce Kelly
#' @param lon longitude of the data
#' @param lat latitude of the data
#' @param z the data
#' @param zlim the limits of the z-axis. A value of NULL indicates that zlim should equal the range of z values.
#' @param pal the name of a palette generating function. Should be a string.
#' @param col.na the color with which any NAs should be drawn with. A value of NA skips this step.
#' @param n the number of distinct colors to request from the palette function.
#' @param refinement the level of bilinear refinement to apply to the image grid
#' @export
add.map.layer = function(map,
lon,
lat,
z,
zlim = NULL,
pal = 'greyscale',
col.na = NA,
n = 255,
refine = 0,
trim = T,
col.low = '',
col.high = '',
rev = F,
indicate = T,
verbose = T) {
## Misc corrections
lon = lon %% 360
map$lon.min = map$lon.min %% 360
map$lon.max = map$lon.max %% 360
lon = as.array(lon)
lat = as.array(lat)
z = as.array(z)
if (length(dim(lon)) < 2) {
if (verbose) { message(' Establishing Grid...', appendLF = F) }
if (length(z) == 1.0 * length(lon) * length(lat)) {
z = array(z, dim = c(length(lon), length(lat)))
lon = matrix(lon, nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(lat, nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
} else { ## TODO add check here that the dimentions are actually correct!
z = array(z, dim = c(length(unique(lon)), length(unique(lat))))
lon = matrix(unique(lon), nrow = dim(z)[1], ncol = dim(z)[2])
lat = matrix(unique(lat), nrow = dim(z)[1], ncol = dim(z)[2], byrow = T)
}
if (verbose) { message(' Done. ') }
}
nz = length(z)
## Trim
if (trim & length(z) > 100) {
if (verbose) { message(' Starting domain trimming... ', appendLF = F)}
usr = par('usr')
corners = expand.grid(lon = c(usr[1], usr[2]),
lat = c(usr[3], usr[4]))
corners = rgdal::project(cbind(corners$lon, corners$lat), proj = map$p, inv = T)
if (corners[1,1] > corners[2,1]) { ## antimeridian
if (verbose) { message(' antimeridian... ', appendLF = F)}
antimeridian = T
} else {
antimeridian = F
lon = lon %% 360
lon[lon > 180] = lon[lon > 180] - 360
}
field = expand.grid(lon = seq(usr[1], usr[2], length.out = 100),
lat = seq(usr[3], usr[4], length.out = 100))
field = rgdal::project(cbind(field$lon, field$lat), proj = map$p, inv = T)
field[,1] = field[,1] %% 360
if (!antimeridian & any(field[,1] > 180)) {
l = which(!is.na(field[,1]) & field[,1] > 180)
field[l,1] = field[l,1] - 360
}
field.lon = range(field[,1], na.rm = T) %% 360
field.lat = range(field[,2], na.rm = T)
## Trim longitude
if (field.lat[1] > -80 & field.lat[2] < 80) { ## only if a pole isn't visible!
if (verbose) { message(' longitude... ', appendLF = F)}
if (corners[1,1] > corners[2,1]) { ## antimeridian
k = apply(lon, 1, function(x) {any(x >= field.lon[1] & x <= field.lon[2])})
} else {
k = apply(lon, 1, function(x) {any(x <= field.lon[2] & x >= field.lon[1])})
}
if (sum(k) > 2) {
z = z[k,]
lon = lon[k,]
lat = lat[k,]
}
}
## Trim latitude
if (verbose) { message(' latitude... ', appendLF = F) }
k = apply(lat, 2, function(x) {any(x >= field.lat[1] & x <= field.lat[2])})
if (sum(k) > 2) {
z = z[,k]
lon = lon[,k]
lat = lat[,k]
}
if (verbose) { message(' complete, n = ', length(z), ' (', 100 - round(100 * length(z) / nz), '% trimmed)')}
}
if (refine > 0) {
if (!antimeridan) {
lon[lon > 180] = lon[lon>180] - 360
}
for (i in 1:refine) {
temp = grid.refinement(lon, lat, z)
lon = temp$x
lat = temp$y
z = temp$z
}
if (verbose) {message(' Refined grid by ', 2^refine,'x.')}
} else if (refine < 0) {
for (i in 1:abs(refine)) {
temp = grid.subsample(lon, lat, z)
lon = temp$x
lat = temp$y
z = temp$z
}
if (verbose) {message(' Refined grid by 1:', 2^-refine,'x.')}
}
if (is.null(zlim)) { zlim = range(pretty(as.numeric(z), na.rm = TRUE)) }
## Order: Not needed?
#l = order(rgdal::project(cbind(lon[,round(dim(lon)[2]/2)], lat[,round(dim(lon)[2]/2)]), p = map$p)[,1])
#l = order(lon[,round(dim(lon)[2]/2)])
#lon = lon[l,]
#lat = lat[l,]
#z = z[l,]
## Color scale
col = make.pal(z, pal = pal, rev = rev, n = n, min = zlim[1], max = zlim[2])
col = array(col, dim = dim(lon))
## Apply out of range colors if provided:
if (is.na(col.low)) { col[z < zlim[1]] = NA }
if (is.na(col.high)) { col[z > zlim[2]] = NA }
## Project and Plot
if (verbose) {message(' Projecting grid...')}
#vertex = calc.vertex(lon, lat)
#xy = rgdal::project(cbind(as.numeric(vertex$x), as.numeric(vertex$y)), p = map$p)
#xy = list(x = array(xy[,1], dim = dim(vertex$x)),
# y = array(xy[,2], dim = dim(vertex$x)))
xy = rgdal::project(cbind(as.numeric(lon), as.numeric(lat)), p = map$p)
xy = list(x = array(xy[,1], dim = dim(lon)),
y = array(xy[,2], dim = dim(lat)))
xy = calc.vertex(xy$x, xy$y)
if (verbose) { message(' Starting plotting... ', appendLF = F) }
for (i in 1:dim(col)[1]) {
for (j in 1:dim(col)[2]) {
if (!is.na(col[i,j])) {
polygon(x = c(xy$x[i,j], xy$x[i,j+1], xy$x[i+1,j+1], xy$x[i+1,j]),
y = c(xy$y[i,j], xy$y[i,j+1], xy$y[i+1,j+1], xy$y[i+1,j]),
col = col[i,j], border = NA)
}
}
}
if (verbose) { message(' complete.') }
## Extras
if (verbose) { message(' Final stats: \tN.low: ', sum(z < zlim[1]), '\tN.high: ', sum(z > zlim[2]), '\tN: ', length(z)) }
if (indicate) {
st = paste0('Data range: (', round(min(z, na.rm = TRUE), 3), ', ', round(max(z, na.rm = TRUE), 3),
') Z range: (', round(zlim[1], 3), ', ', round(zlim[2], 3), ')'
)
mtext(st, line = 0.25, adj = 1, cex = 0.7)
}
}
#' @export
add.map.eez = function(map, eez, names = NULL, col = 'black', lwd = 1, lty = 1, verbose = T) {
if (is.null(names)) {
names = names(eez)
if (verbose) {
message('No names provided, using all available: ', paste0(names, collapse = ', '))
}
}
for (n in names) {
l = which(names(eez) == n)
if (length(l) > 0) {
for (i in 1:length(eez[[n]])) {
add.map.line(map, eez[[n]][[i]]$lon, eez[[n]][[i]]$lat, col = col, lwd = lwd, lty = lty)
}
}
}
}
#' @title Add Map Bathymetry (shaded)
#' @author Thomas Bryce Kelly
#' @param map a map object returned by the function make.map()
#' @param bathy a bathymetric dataframe returned by get.map.bathy
#' @param pal the color palette used for coloring the shading. Should consist of a function name as a string
#' @param n the number of distinct colors to use
#' @param zlim the z-axis range that the colormap is projected upon
#' @param rev a boolean flag to flip the color axis
#' @param filled a flag to turn on color smoothing and filling from the OCE package. Works for some datasets but may cause issues.
#' @param col.low the color used when painting data that is below the zlim range. A value of NA causes the data not to be drawn, and a value of '' uses the lowest color value for out of range data.
#' @param col.high same as \emph{col.low} but for data that is out of range above the zlim.
#' @param refinement the level of bilinear refinement to apply to the image grid
#' @export
add.map.bathy = function(map,
bathy,
refine = 0,
trim = T,
zlim = NULL,
pal = 'greyscale',
n = 255,
rev = FALSE,
col.low = '',
col.high = '',
verbose = T) {
add.map.layer(map,
lon = bathy$Lon,
lat = bathy$Lat,
z = bathy$Z,
refine = refine,
pal = pal,
rev = rev,
zlim = zlim,
col.low = col.low,
col.high = col.high,
indicate = FALSE,
n = n,
trim = trim,
verbose = verbose)
redraw.map(map)
}
#' @title Redraw Map
#' @author Thomas Bryce Kelly
#' @param map a map object as returned by make.map()
#' @export
redraw.map = function(map) {
for (i in 1:length(map$coastline$projected)) {
polygon(map$coastline$projected[[i]][,1], map$coastline$projected[[i]][,2], col = map$land.col)
}
if (map$grid$draw.grid) {
for (i in map$grid$lon.axis$label) {
add.map.line(map, lon = c(i,i), lat = c(-80, 80), col = 'dark grey')
}
for (i in map$grid$lat.axis$label) {
add.map.line(map, lon = c(-180,180), lat = c(i,i), col = 'dark grey')
}
}
box()
}
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.