Bathymetric maps in R

There are a number of different ways to make basic maps in R. In the last year or so I’ve become a big fan of leaflet and the R leaflet package that makes these maps a breeze to build in R. Leaflet makes very nice online interactive maps, but doesn’t provide a great option for a static map like you would put in a publication or presentation. I’ve jumped around between a few different methods for making static maps, but for the purposes of this post I’ll focus on the two that I like and use most.

The first method uses the oce package, developed primarily by Dan Kelley and Clark Richards at my native Dalhousie. The second uses the ggmap package, which is the brainchild of Hadley Wickham and other developers of the tidyverse. My plan is to present simple reproducible examples for each map making method to show how to use them to display basic types of spatial data: contours (bathymetry), lines, points, and polygons.

Get data

Step one is to get the data that we’ll be plotting. A very handy tool for this is the getNOAA.bathy() function from the marmap package, which, perhaps not surprisingly, gets user-selected bathymetric data from the NOAA ETOPO database. There are a number of other great functions within marmap, but, for reasons I won’t get into here, I’ve never quite been satisfied with using it to make maps. Anyway, here’s the data I’ll be using:


# get bathymetry data
b = getNOAA.bathy(lon1 = -68, lon2 = -62, lat1 = 40, lat2 = 45, resolution = 1)
## Querying NOAA database ...
## This may take seconds to minutes, depending on grid size
## Building bathy matrix ...

# make a simple track line
lon = c(-65.17536, -65.37423, -65.64541, -66.06122, -66.15161)
lat = c(43.30837, 42.94679, 42.87448, 42.92871, 42.72985)
lin =, lat)

# make a few points
lon = c(-65.3, -65.7, -64.1)
lat = c(43.4, 43, 42.9)
pts =, lat)

# build a polygon (in this case the 'Roseway Basin Area To Be Avoided')
ply =, 43.266667),
      c(-64.983333, 42.783333),
      c(-65.516667, 42.65),
      c(-66.083333, 42.866667))
colnames(ply) = c('lon', 'lat')


The following map with oce does not use a projection, which basically means that all plotting occurs as if the world were completely flat. Plotting in this way is quick and convenient, especially because it allows you to use all of R’s base graphic functions, but is not appropriate in many cases. I typically use this approach when plotting relatively small areas (maybe 10s of kilometers) and in circumstances where exact distances aren’t crticially important.


# convert bathymetry
bathyLon = as.numeric(rownames(b))
bathyLat = as.numeric(colnames(b))
bathyZ = as.numeric(b)
dim(bathyZ) = dim(b)

# define plotting region
mlon = mean(pts$lon)
mlat = mean(pts$lat)
span = 300

# --- unprojected --- #

# plot coastline
plot(coastlineWorldFine, clon = mlon, clat = mlat, span = span)

# plot bathymetry
        levels = c(-50, -100, -150, -200, -250),
        lwd = c(1, 1, 2, 2, 3),
        lty = c(3, 1, 3, 1, 3),
        drawlabels = F, add = TRUE, col = 'darkgray')

# add depth legend
legend("bottomright", seg.len = 3, cex = 0.8,
       lwd = c(1, 1, 2, 2, 3),
       lty = c(3, 1, 3, 1, 3),
       legend = c("50", "100", "150", "200", "250"),
       col = 'darkgray', title = "Depth [m]", bg= "white")

# add map data
points(pts, pch = 16, col = 'red')
lines(lin, col = 'blue')
polygon(ply, lty = 2)

When the need arises for a projection, as is often the case, the oce package has an impressive number available. Check out Dan’s blog post for a full list. Plotting projected data is slower and a little bit clunkier because you need to use oce custom functions (e.g., mapPoints(), mapLines(), etc), but in my experience this tends to work well. I can’t say the resulting maps are spectacularly pretty, but they’re very accurate and theoretically quite customizable if there’s something you’d like to change or add. Here’s what that looks like:

# --- projected ---#

# plot coastline
plot(coastlineWorldFine, clon = mlon, clat = mlat, span = 200, projection="+proj=merc", col = 'lightgrey')

# plot bathymetry
mapContour(bathyLon, bathyLat, bathyZ,
           levels = c(-50, -100, -150, -200, -250),
           lwd = c(1, 1, 2, 2, 3),
           lty = c(3, 1, 3, 1, 3),
           col = 'darkgray')

# add depth legend
legend("bottomright", seg.len = 3, cex = 0.7,
       lwd = c(1, 1, 2, 2, 3),
       lty = c(3, 1, 3, 1, 3),
       legend = c("50", "100", "150", "200", "250"),
       col = 'darkgray', title = "Depth [m]", bg = "white")

# add map data
mapPoints(longitude = pts$lon, latitude = pts$lat, pch = 16, col = 'red')
mapLines(longitude = lin$lon, latitude = lin$lat, col = 'blue')
mapPolygon(longitude = ply$lon, latitude = ply$lat, lty = 2)


ggmap is the tidyverse answer for plotting spatial data in R. As a result, it doesn’t use R’s basic graphical functions, but uses the grammer of graphics syntax of ggplot2. It offers a lot of powerful features, and the potential to make maps that look really, really nice. Here’s a quick example:


# define center of map
ll_mean = c(mean(pts$lon), mean(pts$lat))

# get raster background map data
ggm = get_map(location = ll_mean,  maptype = "satellite", source = "google", zoom = 8)

# convert bathymetry to data frame
bf = fortify.bathy(b)

# create map
gg = ggmap(ggm) + 
  geom_contour(data = bf, aes(x=x, y=y, z=z),
               colour="white", size=0.1)+
  geom_polygon(data = ply, color = "black", alpha = 0.3) +
  geom_line(data = lin, color = "blue", alpha = 0.8) +
  geom_point(data = pts, colour = "black", fill = "red", stroke = 1, size = 2, alpha = 0.7, shape = 21)+

# adjust theme for nicely plotting in markdown (white axis labs and transparent bg)
gg = gg + theme(axis.text=element_text(colour="gray"))+
  theme(plot.background = element_rect(fill = "transparent",colour = NA))

# show map

Pretty nice, eh?. I particularly like how easily it allows me to overlay data on Google satellite imagery. This is a relatively new discovery for me, and is something that I’ll likely start using more often.

My biggest issue with this map is that I haven’t been able to find an elegant way to add labels to the contour lines. There are ways to do this in ggplot2 (like with the direct.label package), but these fail with projected data used in ggmap. Perhaps I’ll file an issue with the developers to see if they can shed any light.


Both oce and ggmap offer simple, accurate, and elegent tools for plotting bathymetric data. I’m still new to ggmap, but I’m quite impressed with its capabilities. I suspect I’ll continue to use oce for quickly plotting results on the fly, but will definitely begin to use and explore ggmap for presentation and publication quality visuals.


comments powered by Disqus