I am trying to interpolate a raster layer of chlorophyll-a data using kriging (gstat::krige) in R, however the resulting interpolation is showing NA values in a peculiar straight line through the raster (see last plot).

From the interpolated raster I have then extracted values corresponding to spatial points.

The chlorophyll-a raster (chl), colour palette (pal) and spatial points (subSP) relevant to this question can be found here.


#Load subsetted lat and lon coordinates
coordinates(subSP) <- c("lon", "lat")
projection(subSP) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

#set grid xmin, xmax, ymin, ymax of total dataset extent of ALL track points
parGrid <- raster(extent(-72, 52, -68, -21), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
res(parGrid) <-  1.851111 ## ~about 100km?? Not sure if this is correct.
#Crop the parent grid to the extent of the SUBSET of track points (plus a little buffer)
ex <- extent(plyr::round_any(extent(subSP)[1]-2, 1, floor), plyr::round_any(extent(subSP)[2]+2, 1, ceiling), plyr::round_any(extent(subSP)[3]-2, 1, floor), plyr::round_any(extent(subSP)[4]+2, 1, ceiling))
subGrid <- crop(parGrid, ex)

#Read in chlorophyll-a raster (same extent as cropped grid).

#Plot raster to identify missing values
plot(chl, col=pal$cols, breaks = pal$breaks, legend = FALSE, main="Chl-a", xlab="lat", ylab="lon")
mp <- crop(countriesLow, extent(chl))
plot(mp,col="black",border=FALSE, add=T)
plot(subSP,add=T, pch=16, cex=0.3, col="grey50")

enter image description here

#Krige raster layer to interpolate NA values
#Convert chl raster to dataframe
P <- data.frame(X=coordinates(chl)[,1], Y=coordinates(chl)[,2], chl=values(chl))
#Remove NAs from chl raster
P <- P[!is.na(P$chl),]
#Convert to SpatialPoints
coordinates(P) <- c("X", "Y")
#Set projection
projection(P) <- projection(chl)

#Convert empty raster to SpatialGrid
grd <- as.data.frame(coordinates(subGrid))
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE  # Create SpatialPixel object
fullgrid(grd) <- TRUE  # Create SpatialGrid object
proj4string(grd) <- proj4string(P)

# Replace chl point boundary extent with that of empty grid - perhaps this is where the problem orginates??
P@bbox <- grd@bbox

# Define the 1st order polynomial equation
f.1 <- as.formula(chl ~ X + Y) 
# Compute the sample variogram
var.smpl <- variogram(f.1, P)
# Compute the variogram model by passing the nugget, sill and range values.
dat.fit  <- fit.variogram(var.smpl, fit.ranges = FALSE, fit.sills = FALSE, vgm(psill=0.65, model="Exp", range=800, nugget=0.15))
# Plot variogram 
plot(var.smpl, dat.fit, xlim=range(var.smpl$dist), ylim=range(var.smpl$gamma))

enter image description here

# Perform the krige interpolation
dat.krg <- krige(f.1, P, grd, dat.fit)
# Convert kriged surface to a raster object
r <- raster(dat.krg)

#Plot interpolated data
plot(r, col=pal$cols, breaks = pal$breaks, legend = FALSE, main="Chl-a interpolated", xlab="lat", ylab="lon")

enter image description here

What is causing the straight line of NA values on the left side of the raster plot and how do I fix this?

#Extract chl from interpolated raster values
chl_ex <- raster::extract(r, subSP)
  • First thing I'd check is plot(r) and plot(is.na(r)) just to make sure the white pixels really are NA, and not just outside of the palette breaks. (Can't check your files rn, but can have a closer look later). – mdsumner Jul 3 at 2:53
  • Thanks @mdsumner - double checked and yep they are NAs. – JaimieC Jul 3 at 2:57
  • The returned values from krige are actually NaN rather than NA. So numerical problem rather than missing values... – Spacedman Jul 3 at 10:23

If you reduce the cell size it seems to fill the missing values, and a way to guesstimate roughly 100km * 100km is to check the sqrt of the area of pixels, here it varies from 80km to 128km. (They are always taller than they are wide so might depend if balancing area or dimension is more important)

I don't know why it gets missing values with lower res, however.

res(parGrid) <-  1.2 ## ~about 100km?? Not sure if this is correct.
class      : RasterLayer 
dimensions : 39, 103, 4017  (nrow, ncol, ncell)
resolution : 1.2, 1.2  (x, y)
extent     : -72, 51.6, -67.8, -21  (xmin, xmax, ymin, ymax)
source     : memory
names      : layer 
values     : 82.99557, 128.4049  (min, max)

Then I get what seems reasonable:

enter image description here

  • Thanks @mdsumner! It would be good to work out why this is happening because I am aggregating data over many seasons and it may arise again. Very bizarre that its occurring in a straight line. I thought maybe it had to do with the extent/bbox because its close to the edge?! – JaimieC Jul 3 at 4:05
  • 1
    Note that with this setup, gstat uses great circle distances to create the covariance matrix. Are you sure this variogram model is valid on spherical / ellipsoidal distance measures? – Edzer Pebesma Jul 3 at 12:12
  • @mdsumner - I'm still seeing the same problem with higher resolution. I don't understand why I see a straight line of NaNs if the variogram model isn't valid. This is my first time kriging in R. – JaimieC Jul 4 at 2:16
  • If anyone has any thoughts on what I should try next, that would be appreciated. Just a newbie here :) – JaimieC Jul 10 at 4:05
  • Ultimately I think kriging is unsuitable, the missing data extends in time as well as space so a good compromise requires going back to the daily data, and depends on the extraction required - getting nearest values etc. I don't think any interpolation can work here, and is why this isn't solved yet. (Requires conversation) – mdsumner Jul 10 at 8:23

Your Answer

