Importing bathymetry and coastline data in R

After noticing some frustrating inaccuracies with the high-resolution world coastlines and national boundaries database found in worldHires from the package mapdata (based on CIA World Data Bank II data), I decided to look into other options. Although listed as "depreciated", the data found in NOAAs online "Coastline Extractor" is a big step forward. There seem to be more up-to-date products, but this served my needs for the moment, and I thought I'd pass along the address to other R users. I exported the data in ASCII "Matlab" format, which is basically just a 2-column text file (.dat) with NaN's in the rows that separate line segments.

I've also discovered the bathymetry / topography data from GEBCO. Again, very easy to import into R from the netCDF files.

The above map of the Galapagos Archipelago illustrates the quality of both datasets. It also shows the comparison of coastline accuracy between World Vector Shoreline (1:250,000), world (map package), and worldHires (mapdata package) datasets. Obviously, the low-resolution world data only makes sense for quick plotting at large scales, but the high-resolution data is as much as 1/10° off in some locations. I noticed these errors for the first time when trying to map some data for smaller coastal bays. It drove me crazy trying to figure out where the errors were - in my data locations or the map itself. Bathymetry used in the map was 30 arc-second resolution GEBCO data.

[EDIT: The comparison of coastline data now includes the high resolution data from the rworldmap package.]

A more detailed description export settings:

• Coastline data (from 'Coastline Extractor') :
• Coastline database: World Vector Shoreline (1:250,000)
• Compression method for extracted ASCII data: None
• Coast Format options: Matlab
• Coast Preview options: GMT Plot
• Bathymetry / topography data [link]:
• General Bathymetric Chart of the Oceans (GEBCO) : GEBCO_08 Grid (30 arc-second resolution)

Here's another example of bathymetry / topography data for the western Pacific (1 minute resolution GEBCO data):

For both maps, I took inspiration for the color palettes from GMT. The rgb color levels of these palettes have got to be documented somwhere, but I gave up looking after a while and managed to hack their levels from color scales contained in .png files [link].

Below is the R code to reproduce the figures.

To reproduce Galapagos coastline example:

###required packages
library(RNetCDF)
library(maps)
library(mapdata)
library(rworldmap)
library(rworldxtra)
 
###Data
#data locations
bathy_fname <- "galapagos_gebco_08_-92_-2_-88_2.nc" # from https://www.bodc.ac.uk/data/online_delivery/gebco/gebco_08_grid/
coast_fname <- "galapagos_18563.dat" # from http://www.ngdc.noaa.gov/mgg/coast/getcoast.html
 
#load bathy data
nc <- open.nc(bathy_fname)
print.nc(nc)
tmp <- read.nc(nc)
z <- array(tmp$z, dim=tmp$dim)
#z[which(z > 0)] <- NaN
z <- z[,rev(seq(ncol(z)))]
xran <- tmp$x_range
yran <- tmp$y_range
zran <- tmp$z_range
lon <- seq(tmp$x[1], tmp$x[2], tmp$spac[1])
lat <- seq(tmp$y[1], tmp$y[2], tmp$spac[1])
rm(tmp)
close.nc(nc)
 
#load coast data
coast <- read.table(coast_fname)
names(coast) <- c("lon", "lat")
 
###Plot
#make palette
ocean.pal <- colorRampPalette(
 c("#000000", "#000413", "#000728", "#002650", "#005E8C",
 "#0096C8", "#45BCBB", "#8AE2AE", "#BCF8B9", "#DBFBDC")
)
 
land.pal <- colorRampPalette(
 c("#467832", "#887438", "#B19D48", "#DBC758", "#FAE769",
 "#FAEB7E", "#FCED93", "#FCF1A7", "#FCF6C1", "#FDFAE0")
)
 
zbreaks <- seq(-8000, 8000, by=100)
cols <-c(ocean.pal(sum(zbreaks<=0)-1), land.pal(sum(zbreaks>0)))
 
#compare coastlines to package 'mapdata' and 'rworldxtra'
LWD <- 1
COL <- c(1:4)
png("coastline_compare2.png", width=7.5, height=6, units="in", res=600)
#quartz(width=7.5, height=6)
layout(matrix(1:2, 1,2), widths=c(6,1.5), heights=c(6))
 
par(mar=c(2,2,1,1), ps=10)
image(lon, lat, z=z, col=cols, breaks=zbreaks, useRaster=TRUE, ylim=c(-1.5,0.5), xlim=c(-92,-90))
lines(coast, col=COL[1], lwd=LWD)
map("world", col=COL[2], ylim=c(-2,2), xlim=c(-93,-88), add=TRUE, lwd=LWD)
map("worldHires", col=COL[3], ylim=c(-2,2), xlim=c(-93,-88), lwd=LWD, add=TRUE)
plot(getMap(resolution='high'), border=COL[4], lwd=LWD, add=TRUE)
 
legend("bottomleft", legend=c("World Vector Shoreline", "maps: world", "mapdata: worldHires", "rworldmap: high res."), lty=1, lwd=LWD, col=COL, bg="white") 
 
par(mar=c(2,0,1,5))
image(x=1, y=zbreaks, z=matrix(zbreaks, 1, length(zbreaks)), col=cols, breaks=zbreaks, useRaster=TRUE, xlab="", ylab="", axes=FALSE)
axis(4, at=seq(-8000, 8000, 1000), las=2)
mtext("[meters]", side=4, line=3)
box()
 
dev.off()

Created by Pretty R at inside-R.org





To reproduce western Pacific map:

###required packages
library(RNetCDF)
 
###Data
#data location
bathy_fname <- "west_pac_gebco_1min_100_-20_170_50.nc" # from https://www.bodc.ac.uk/data/online_delivery/gebco/gebco_08_grid/
 
#load bathy data
nc <- open.nc(bathy_fname)
print.nc(nc)
tmp <- read.nc(nc)
z <- array(tmp$z, dim=tmp$dim)
#z[which(z > 0)] <- NaN
z <- z[,rev(seq(ncol(z)))]
xran <- tmp$x_range
yran <- tmp$y_range
zran <- tmp$z_range
lon <- seq(tmp$x[1], tmp$x[2], tmp$spac[1])
lat <- seq(tmp$y[1], tmp$y[2], tmp$spac[1])
rm(tmp)
close.nc(nc)
 
###Plot
#make palette
ocean.pal <- colorRampPalette(
 c("#000000", "#000209", "#000413", "#00061E", "#000728", "#000932", "#002650", 
 "#00426E", "#005E8C", "#007AAA", "#0096C8", "#22A9C2", "#45BCBB", 
 "#67CFB5", "#8AE2AE", "#ACF6A8", "#BCF8B9", "#CBF9CA", "#DBFBDC", 
 "#EBFDED")
)
 
land.pal <- colorRampPalette(
 c("#336600", "#F3CA89", "#D9A627", 
 "#A49019", "#9F7B0D", "#996600", "#B27676", "#C2B0B0", "#E5E5E5", 
 "#FFFFFF")
)
 
zbreaks <- seq(-11000, 7000, by=10)
cols <-c(ocean.pal(sum(zbreaks<=0)-1), land.pal(sum(zbreaks>0)))
 
#compare coastlines to package 'mapdata'
png("west_pac.png", width=7.5, height=6, units="in", res=200)
#quartz(width=7.5, height=6)
layout(matrix(1:2, 1,2), widths=c(6,1.5), heights=c(6))
 
par(mar=c(2,2,1,1), ps=10)
image(lon, lat, z=z, col=cols, breaks=zbreaks, useRaster=TRUE, ylim=c(-20,50), xlim=c(100,170))
 
par(mar=c(2,0,1,5))
image(x=1, y=zbreaks, z=matrix(zbreaks, 1, length(zbreaks)), col=cols, breaks=zbreaks, useRaster=TRUE, xlab="", ylab="", axes=FALSE)
axis(4, at=seq(-11000, 7000, 1000), las=2)
mtext("[meters]", side=4, line=3)
box()
 
dev.off()

Created by Pretty R at inside-R.org