Creating a landmask for the West Antarctic Peninsula in R

This is going to be a pretty niche topic, but probably useful for someone out there.  Lately I’ve been working with a lot of geospatial data for the West Antarctic Peninsula.  One of the things that I needed to do was krig the data (krigging is a form of 2D interpolation, I’m using the pracma library for this).  Krigging is a problem near coastlines because it assumes a contiguous space to work in.  If there happens to be an island or other landmass in the way there is no way to represent the resulting discontinuity in whatever parameter you’re looking at.  Because of this I needed to find a way to mask the output.  This doesn’t really solve the problem, but at least it allows me to identify areas of concern (for example interpolation that extends across an isthmus, if there are sample points only on one side).

I’m krigging and building maps entirely inside R, which has somewhat immature packages for dealing with geospatial data.  The easiest masking solution would be to use filled polygons from any polygon format shapefile that accurately represents the coastline.  Unfortunately I couldn’t find an R package that does this correctly with the shapefiles that I have access too.  In addition, because of my downstream analysis it was better to mask the data itself, and not just block out landmasses in the graphical output.

Sharon Stammerjohn at the Institute of Arctic and Alpine Research pointed me to the excellent Bathymetry and Global Relief dataset produced by NOAA.  This wasn’t a complete solution to the problem but it got me moving in the right direction.  From the custom grid extractor at http://maps.ngdc.noaa.gov/viewers/wcs-client/ I selected a ETOPO1 (bedrock) grid along the WAP, with xyz as the output format.  If you’re struggling with shapefiles the xyz format is like a cool drink of water, being a 3-column matrix of longitude, latitude, and height (or depth).  For the purpose of creating the mask I considered landmass as any lat-long combination with height > 0.

There is one more really, really big twist to what I was trying to do, however.  The Palmer LTER uses a custom 1 km pixel grid instead of latitude-longitude.  It’s a little easier to conceptualize than lat-long given the large longitude distortions at high latitude (and the inconvenient regional convergence of lat-long values on similar negative numbers).  It is also a little more ecologically relevant, being organized parallel to the coastline instead of north to south.  Unfortunately this makes the grid completely incompatible with other Euclidean reference systems such as UTM.  So before I could use my xyz file to construct a land mask I needed to convert it to the line-station grid system used by the Palmer LTER.  If you’re working in lat-long space you can skip over this part.

Many moons ago someone wrote a Matlab script to convert lat-long to line-station which you can find here.  Unfortunately I’m not a Matlab user, nor am I inclined to become one.  Fortunately it was pretty straightforward to copy-paste the code into R and fix the syntatic differences between the two languages.  Three functions in total are required:

## AUTHORS OF ORIGINAL MATLAB SCRIPT:
#   Richard A. Iannuzzi
#   Lamont-Doherty Earth Observatory
#   iannuzzi@ldeo.columbia.edu
#   based on: LTERGRID program written by Kirk Waters (NOAA Coastal Services Center), February 1997

## some functions that are used by the main function

SetStation <- function(e, n, CENTEREAST, CENTERNORTH, ANGLE){
  uu = e - CENTEREAST
  vv = n - CENTERNORTH
  z1 = cos(ANGLE)
  z2 = sin(ANGLE)
  NorthKm = (z1 * uu - z2 *vv) / 1000 + 600
  EastKm = (z2 * uu + z1 * vv) / 1000 + 40
  
  return(c(NorthKm, EastKm))
}

CentralMeridian <- function(iz){
  if(abs(iz) > 30){
    iutz = abs(iz) - 30
    cm = ((iutz * 6.0) -3.0) * -3600
  }
  else{
    iutz = 30 - abs(iz)
    cm = ((iutz * 6.0) +3.0) * +3600
  }
  return(cm)
}

GeoToUTM <- function(lat, lon, zone){

  axis = c(6378206.4,6378249.145,6377397.155,
          6378157.5,6378388.,6378135.,6377276.3452,
          6378145.,6378137.,6377563.396,6377304.063,
          6377341.89,6376896.0,6378155.0,6378160.,
          6378245.,6378270.,6378166.,6378150.)
  
  bxis = c(6356583.8,6356514.86955,6356078.96284,
          6356772.2,6356911.94613,6356750.519915,6356075.4133,
          6356759.769356,6356752.31414,6356256.91,6356103.039,
          6356036.143,6355834.8467,6356773.3205,6356774.719,
          6356863.0188,6356794.343479,6356784.283666,
          6356768.337303)
  
  ak0 = 0.9996
  radsec = 206264.8062470964  
  sphere = 9
  
  a = axis[sphere - 1]             # major axis size
  b = bxis[sphere - 1]             # minior axis size
  es = ((1-b^2/a^2)^(1/2))^2       # eccentricity squared
  slat = lat * 3600                # latitude in seconds
  slon = -lon * 3600               # longitude in seconds
  cm = 0                           # central meridian in sec
  iutz = 0
  
  cm = CentralMeridian(zone)       # call the function
  
  phi = slat/radsec
  dlam = -(slon - cm)/radsec
  epri = es/(1.0 - es)
  en = a/sqrt(1.0 - es * sin(phi)^2)
  t = tan(phi)^2
  c = epri * cos(phi)^2
  aa = dlam * cos(phi)
  s2 = sin(2.0 * phi)
  s4 = sin(4.0 * phi)
  s6 = sin(6.0 * phi)
  f1 = (1.0 - (es/4.)-(3.0*es*es/64.)-(5.0*es*es*es/256))
  f2 = ((3*es/8)+(3.0*es*es/32)+(45*es*es*es/1024))
  f3 = ((15*es*es/256)+(45*es*es*es/1024))
  f4 = (35*es*es*es/3072)
  em = a*(f1*phi - f2*s2 + f3*s4 - f4*s6)
  xx = ak0 * en * (aa + (1.-t+c) * aa^3/6 + (5 - 18*t + t*t + 72*c-58*epri)* aa^5/120) + 5e5
  yy = ak0 * (em + en * tan(phi) *((aa*aa/2) + (5-t+9*c+4*c*c)*aa^4/24 + (61-58*t +t*t +600*c - 330*epri)* aa^6/720))
  
  if(zone < 0 | slat < 0){
    yy = yy + 1e7
  }
  
  return(c(xx, yy))
}

## This function actually works with your data

ll2gridLter <- function(inlat, inlon){
  NorthKm = 0           # initialize
  EastKm = 0            # initialize
  zone = -20            # set zone (for LTER region, I think)
  ANGLE = -50 * pi / 180
  CENTEREAST = 433820.404        # eastings for station 600.040
  CENTERNORTH = 2798242.817      # northings for station 600.040
  
  # take latitude longitude and get station
  x.y = GeoToUTM(inlat, inlon, zone)
  NorthKm.EastKm = SetStation(x.y[1], x.y[2], CENTEREAST, CENTERNORTH, ANGLE)
  return(NorthKm.EastKm)
}

Once the functions are defined I used them to convert the lat/long coordinates in the xyz file to line-station.

## Read in xyz file.
lat.long.depth <- read.table('etopo1_bedrock.xyz', header = F, col.names = c('long', 'lat', 'depth'))

## Limit to points above sea level.
lat.long.land <- lat.long.depth[which(lat.long.depth$depth >= 0),]

## Create a matrix to hold the output.
line.station.land <- matrix(ncol = 3, nrow = length(lat.long.land$long))
colnames(line.station.depth) <- c('line', 'station', 'depth')

## Execute the ll2gridLter function on each point.  Yes, I'm using a loop to do this.
for(i in 1:length(lat.long.land$long)){
  line.station.land[i,] <- c(ll2gridLter(lat.long.land$lat[i], lat.long.land$long[i]), lat.long.land$depth[i])
  print(paste(c(i, line.station.land[i,])))
}

## Write out the matrix.
write.csv(line.station.land, 'palmer_grid_landmask.csv', row.names = F, quote = F)

At this point I had a nice csv file with line, station, and elevation.  I was able to read this into my existing krigging script and convert into a mask.

## Read in csv file.
landmask <- read.csv('palmer_grid_landmask.csv')

## Limit to the lines and stations that I'm interested in.
landmask <- landmask[which(landmask[,1] <= 825 & landmask[,1] >= -125),]
landmask <- landmask[which(landmask[,2] <= 285 & landmask[,2] >= -25),]

## Interpolated data is at 1 km resolution, need to round off
## to same resolution here.
landmask.expand <- cbind(ceiling(landmask[,1]), ceiling(landmask[,2]))

## Unfortunately this doesn't adequately mask the land.  Need to expand the size of each
## masked pixel 1 km in each direction.
landmask.expand <- rbind(landmask.expand, cbind(floor(landmask[,1]), floor(landmask[,2])))
landmask.expand <- rbind(landmask.expand, cbind(ceiling(landmask[,1]), floor(landmask[,2])))
landmask.expand <- rbind(landmask.expand, cbind(floor(landmask[,1]), ceiling(landmask[,2])))
landmask.expand <- unique(landmask.expand)

I’m not going to cover how I did the krigging in this post.  My krigged data is in matrix called temp.pred.matrix with colnames given by ‘x’ followed by ‘station’, as in x20 for station 20, and row names ‘y’ followed by ‘line’, as in y100 for line 100.  To convert interpolated points that are actually land to NA values I simply added this line to my code:

temp.pred.matrix[cbind(paste0('y', landmask.expand[,1]), paste0('x', landmask.expand[,2] * -1))]

Here’s what the krigged silicate data looks like after masking.

Excellent.  The masked area corresponds with known landmasses; that’s Adelaide Island (home of Rothera Station) coming in at the eastern portion of Line 300, and various islets and the western edge of the Antarctic Peninsula to the northeast.  At this point erroneous data has been eliminated from the matrix.  Annual inventories of silicate and other nutrients can be more accurately calculated form the data and our eyes are not drawn to interesting features in the interpolation that have no chance of reflecting reality because they are over land.  The white doesn’t look that appealing to me in this plot however, so I masked the land with black by adding points to the plot.  Again, I’m not going to show the whole plotting routine because some variables would require a lengthy explanation about how my larger dataset is structured.  The plot was created using imagep in the oce package.  One thing to be aware of is that this command automatically transposes the matrix, thus mask needs to be transposed as well.  I did this manually in the following line:

## Show masked points as black squares.
points(landmask.expand[,1] ~ {-1 * landmask.expand[,2]}, pch = 15, cex = 0.6)

And the final plot: