Plotting spatial data is a common task in research, but this requires a set of technical / design tools that are not always available to all. Here we provide step-by-step instructions for creating a map representing global sea surface oxygen.

The code provided can be used to map additional spatial information in the same raster format. The Robinson projection used is a good compromise for the problem of representing the entire globe as a flat image.

First, you need a shapefile defining global landmasses and a raster file containing the oxygen conditions. These are available at:

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

# 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)

# 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")

# 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') )

# 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