|
############################################### |
|
### Example Southern ocean map using ggplot ### |
|
############################################### |
|
|
|
# WJ Grecian 2022-06-23 |
|
|
|
# load libraries |
|
require(tidyverse) |
|
require(sf) |
|
require(raster) |
|
require(marmap) |
|
require(patchwork) |
|
|
|
# define projection |
|
prj <- "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=km +no_defs +ellps=WGS84 +towgs84=0,0,0" |
|
|
|
# download bathymetry data from marmap |
|
bathy <- marmap::getNOAA.bathy(lon1 = -180, |
|
lon2 = 180, |
|
lat1 = -90, |
|
lat2 = -30, |
|
resolution = 10) |
|
# convert bathmetry data to raster |
|
bat <- marmap::as.raster(bathy) |
|
bat[bat > 0] <- NA # not interested in values on land... |
|
|
|
# project bathymetry data |
|
bat_prj <- raster::projectRaster(bat, |
|
res = 50, |
|
method = "bilinear", |
|
crs = prj) |
|
|
|
# convert bathymetry data to xyz tibble for plotting in ggplot2 |
|
bat_df <- bat_prj %>% raster::rasterToPoints() %>% as_tibble() |
|
bat_df # check format looks right |
|
names(bat_df)[3] <- "depth" # change variable name |
|
|
|
# load shapefile from rnatural earth package |
|
world_shp <- rnaturalearth::ne_countries(scale = 110, returnclass = "sf") |
|
|
|
# define cropping polygon to clip shapefile to southern ocean |
|
CP <- sf::st_bbox(c(xmin = -180, xmax = 180, ymin = -90, ymax = 0), |
|
crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") %>% sf::st_as_sfc() |
|
|
|
# crop shapefile to southern hemisphere and project |
|
sf::sf_use_s2(FALSE) # switch off strange new sf warnings |
|
world_shp <- sf::st_crop(sf::st_buffer(world_shp, 0), CP) |
|
world_shp <- world_shp %>% sf::st_transform(prj) |
|
|
|
# quick example plot |
|
p1 <- ggplot() + |
|
theme_bw() + |
|
geom_raster(aes(x = x, y = y, fill = depth), data = bat_df) + |
|
geom_sf(aes(), data = world_shp, colour = "dark grey", fill = "dark grey") + |
|
coord_sf(xlim = c(-6500, 6500), ylim = c(-6500, 6500), crs = prj, expand = T) + |
|
xlab("") + ylab("") |
|
|
|
# make a 'pretty' circular plot |
|
# this requires making a masking shape that sits over the plot |
|
# the hole in the middle of the shape lets you see through |
|
# it's a bodge... |
|
|
|
# create a shape that is the same size as the visible hole you want |
|
# in this case I want to see up to 7000 km from the centre of the map |
|
# the projection has the south pole xy coordinate as 0,0 |
|
little <- st_point(c(0,0)) %>% st_sfc(crs = prj) %>% st_buffer(dist = 7000) |
|
|
|
# create enclosing rectangle to mask the map with |
|
xlim <- sf::st_bbox(little)[c("xmin", "xmax")]*1.5 |
|
ylim <- sf::st_bbox(little)[c("ymin", "ymax")]*1.5 |
|
# turn into enclosing rectangle |
|
encl_rect <- list(cbind(c(xlim[1], xlim[2], xlim[2], xlim[1], xlim[1]), |
|
c(ylim[1], ylim[1], ylim[2], ylim[2], ylim[1]))) %>% |
|
sf::st_polygon() %>% |
|
sf::st_sfc(crs = prj) |
|
|
|
# now cookie cut out the circle from the enclosing rectangle and the earth outline |
|
cookie <- sf::st_difference(encl_rect, little) |
|
|
|
p2 <- ggplot() + |
|
theme_void() + |
|
geom_raster(aes(x = x, y = y, fill = depth), data = bat_df) + |
|
geom_sf(aes(), data = world_shp, colour = "dark grey", fill = "dark grey") + |
|
geom_sf(aes(), data = cookie, fill = "white") + |
|
coord_sf(xlim = c(-6500, 6500), ylim = c(-6500, 6500), crs = prj, expand = T) + |
|
xlab("") + ylab("") |
|
|
|
p1 + p2 |
|
|
|
# ends |
