Producing a global map of Sea Surface Oxygen in R

2 min readDec 13, 2021

Plotting spatial data is a common task in research, but this requires a set of technical / design tools, which sometimes are not available to all. Here I provide a step-by-step guide to help producing an equal area map of sea surface oxygen in R.

The code provided can be used to map additional spatial information in the same raster format. The adopted Robinson projection finds a good compromise to the problem of readily showing the whole globe as a flat image.

# 1. The first think is to get a shapefile defining global landmasses and a raster file with oxygen conditions. These are available at:

https://github.com/jorgeassis/dataVisualization/tree/master/Projects/seaSurfaceOxygen/Data

# 2. Then set the working directory in R and load the packages.

setwd("Projects/seaSurfaceOxygen")library(ggplot2)
library(raster)
library(rnaturalearth)
library(sf)
library(geosphere)
library(rgeos)

# 3. Choose the projection and load the spatial data.

The Robinson projection is a projection of the world that attempts to find a good compromise of showing the whole globe as a flat image.

projection <- CRS("+proj=robin +over")oxygenLayer <- raster("Data/DissolvedMolecularOxygen.tif")
worldMap <- shapefile("Data/LandMass Polygon.shp")

# 4. Generate a bounding box and clip spatial data.

bb <- sf::st_union(sf::st_make_grid(
  st_bbox(c(xmin = -180,
            xmax = 180,
            ymax = 90,
            ymin = -90), crs = st_crs(4326)), n = 100))
bb <- st_transform(bb, projection)oxygenLayer <- projectRaster(oxygenLayer, crs = projection)
oxygenLayer <- mask(oxygenLayer, as(bb, "Spatial"))
oxygenLayer <- as.data.frame(oxygenLayer,xy=TRUE,na.rm=T)
colnames(oxygenLayer) <- c("Lon","Lat","Val")worldMap <- spTransform(worldMap, CRSobj = projection)
worldMap <- gBuffer(worldMap, byid=TRUE, width=0.001)
worldMap <- crop(worldMap, as(bb, "Spatial"))

# 5. Define the theme for the final plot.

theme_map <- 
  theme(
    text = element_text(family = "Helvetica", color = "#22211d"),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.major = element_line(color = "black", size = 0.1),
    panel.grid.minor = element_blank(),
    plot.background = element_rect(fill = "#FFFFFF", color = NA),
    panel.background = element_rect(fill = "#FFFFFF", color = NA),
    panel.border = element_blank(),
    legend.background = element_rect(fill = "#FFFFFF", color = NA),
    legend.position="bottom", 
    legend.box = "horizontal",
    legend.margin=margin(0,0,0,0),
    legend.box.margin=margin(-10,-10,-10,-10),
    legend.key.height= unit(0.25, 'cm'),
    legend.key.width= unit(0.75, 'cm') )

# 6. Define a color ramp for the oxygen gradient and plot the map.

myColors <- c("#6FBBE8","#A1ECD8","#F6F9AB","#FCB46D","#B21414","#D278E4","#9914B3")plot <- ggplot() +
  geom_tile(data = oxygenLayer, aes(x=Lon,y=Lat,fill=Val)) +
  scale_fill_gradientn(guide = guide_legend(title="Sea surface dissolved oxygen [mmol.m3]", direction = "horizontal", title.position = "top", title.hjust = 0.5), colours=myColors, na.value='transparent') +
  geom_polygon(data = worldMap, aes(x = long, y = lat, group = group), fill="#A1A1A1", colour = "#A1A1A1" , size=0.25 ) +
  geom_sf(data = bb,fill=NA, colour = "white" , linetype='solid', size=2 ) +
  theme_map
plot

Producing a global map of Sea Surface Oxygen in R

Written by biodiversityDS.