Three-dimensional interpolation of temperature profiles using gstat in R?

Question

I'm working on a problem in which I'd like to interpolate a set of oceanographic temperature profiles over a three-dimensional matrix of unsampled x, y, & z locations. The data to interpolate is housed in a SpatialPointsDataFrame (WOD) in which I've stored temperature data (WOD$TEMP) in association with coordinate information in the form of latitude (WOD$LAT,in degrees),longitude (WOD$LONG, in degrees), and depth (WOD$DEPTH, in meters). The data consists of multiple temperature measurements collected on different depth levels at each station (i.e. multiple TEMP and DEPTH values for each LAT and LONG pair). Here's an sample of the data which is drawn from the World Ocean Database (I can send along a more complete dataset as a .csv if you're interested)

   head(WOD[WOD$YEAR%in%2013 & WOD$TYPE%in%"PFL",c("LAT","LONG","DEPTH","TEMP")],100)[seq(1,100,2),]  
           LAT    LONG   DEPTH   TEMP  
4004511 27.644 -78.732    3.68 23.151  
4004711 27.644 -78.732   14.20 23.133  
4004911 27.644 -78.732   22.75 23.128  
4005111 27.644 -78.732   33.77 23.116  
4005311 27.644 -78.732   42.41 23.119  
4005511 27.644 -78.732   53.04 23.115  
4005711 27.644 -78.732   63.96 23.086  
4005911 27.644 -78.732   72.50 23.045  
4006111 27.644 -78.732   83.82 23.043  
4006311 27.644 -78.732   92.36 23.022  
4006511 27.644 -78.732  103.29 22.977  
4006711 27.644 -78.732  116.59 22.967  
4006911 27.644 -78.732  137.54 22.618  
4007111 27.644 -78.732  176.15 20.855  
4007311 27.644 -78.732  222.69 19.743  
4007511 27.644 -78.732  272.30 18.979  
4007711 27.644 -78.732  320.31 18.334  
4007911 27.644 -78.732  371.48 17.855  
4008111 27.644 -78.732  445.34 16.469  
4008311 27.644 -78.732  543.75 13.591  
4008511 27.644 -78.732  663.89  9.110  
4011811 26.952 -77.340   10.33 23.660  
4012011 26.952 -77.340   20.07 23.699  
4012211 26.952 -77.340   30.00 23.693  
4012411 26.952 -77.340   39.83 23.683  
4012611 26.952 -77.340   50.16 23.664  
4012811 26.952 -77.340   59.50 23.641  
4013011 26.952 -77.340   69.53 23.601  
4013211 26.952 -77.340   80.05 23.580  
4013411 26.952 -77.340   89.79 23.552  
401369  26.952 -77.340   99.52 23.625  
401382  26.952 -77.340  109.45 23.349  
4014011 26.952 -77.340  129.60 22.105  
4014211 26.952 -77.340  158.99 20.615  
4014411 26.952 -77.340  198.49 19.640  
4014611 26.952 -77.340  249.10 18.740  
4014811 26.952 -77.340  298.21 18.209  
4015011 26.952 -77.340  347.80 17.910  
4015211 26.952 -77.340  397.58 17.292  
4015411 26.952 -77.340  496.51 15.465  
4015611 26.952 -77.340  594.80 13.275  
4015811 26.952 -77.340  744.02  9.963  
4016011 26.952 -77.340  892.25  6.338  
4016211 26.952 -77.340 1090.28  5.297  
...  

   WOD=WOD[WOD$YEAR%in%2013 & WOD$TYPE%in%"PFL" & WOD$REGION%in%c("TOTO","SEPC","NWPN","NWPS"),]

   coordinates(WOD) = ~LONG+LAT+DEPTH


   #Based on the demo(gstat3D) example I have fitted the following variogram models:

   library(sp);library(gstat);library(raster);library(rgdal)

   WOD_vg=variogram(TEMP~1,WOD)
   WOD_vgm=fit.variogram(WOD_vg,mod=vgm(psill=2,"Exp",range=1))
   WOD_vgm1=fit.variogram(WOD_vg,mod=vgm(psill=1.8,"Sph",range=1.0,nugget=0.6))
   WOD_vgm2=fit.variogram(WOD_vg,mod=vgm(psill=2,"Sph",range=1.5,nugget=0.6))
   WOD_vgm3=fit.variogram(WOD_vg,mod=vgm(psill=2,"Mat",range=1.5,nugget=0.6))

   plot(WOD_vg$dist,WOD_vg$gamma)
   lines(variogramLine(WOD_vgm,1.5)$dist,variogramLine(WOD_vgm,1.5)$gamma)
   lines(variogramLine(WOD_vgm1,1.5)$dist,variogramLine(WOD_vgm1,1.5)$gamma)
   lines(variogramLine(WOD_vgm2,1.5)$dist,variogramLine(WOD_vgm2,1.5)$gamma)
   lines(variogramLine(WOD_vgm3,1.5)$dist,variogramLine(WOD_vgm3,1.5)$gamma)

   # and attempted to krige the data and variogram model over three dimensional SpatialPixels defined in terms of ~LONG+LAT+DEPTH

   grid3D <- expand.grid(LONG=seq(from = -79.00, to = -76.50, by=0.05),
                         LAT=seq(from = 23.25, to = 26.75, by=0.05),
                         DEPTH=seq(from = 0, to = 1200, by=50))
   gridded(grid3D) = ~LONG+LAT+DEPTH

   res3D <- krige(formula = TEMP ~ 1, WOD, grid3D, model = WOD_vgm2)

#This produces the following error:

> res3D <- krige(formula = TEMP ~ 1, WOD, grid3D, model = WOD_vgm2)
[using ordinary kriging]

"chfactor.c", line 131: singular matrix in function LDLfactor()
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim,  : 
  LDLfactor

This error I gather from reading lots of earlier posts on various forums probably relates to having multiple observations (at different depth levels) with same LAT and LONG coordinates (i.e. horizontal distance between these observations = 0).

I'm relatively new to the world of geostatistics so I have a number of questions:

1) Is 3D kriging using gstat an appropriate method of interpolation in this context? If not is there an alternative approach that I should pursue?

2) How should I deal with having observations on multiple depth levels at the same LAT and LONG coordinates?

3) Do latitude, longitude and depth need to be in the same units (e.g. meters) for accurate estimation of distances in variogram calculation and kriging interpolation?

4)How should I set initial parameters for the anisotropy term given that temperature varies by almost an order of magnitude (4-30 C) over 1000m in the vertical direction but <=1 degree C over >100 km in the horizontal dimensions (more with latitude than with longitude)?

Answer 1

Answer to question 1: it is possibly appropriate, and there are alternatives, and some of them have been implemented rather recently in gstat, see this vignette. Q 2: there is a function zerodist in sp; if the collocated data have different time stemps, it should not be a problem if you do proper ST kriging. Q 3: if gstat is informed that coordinates are long lat (i.e., is.projected(obj) returns FALSE) it will take care of this (and compute great circle distances, with units km). It is up to you to make sure the variogram models chosen are valid on the sphere, though!

Answer 2

gstat can do spatial-temporal kriging, as explained here: https://www.r-bloggers.com/spatio-temporal-kriging-in-r/ . Time is just a third dimension, and you could use "elevation" data instead of time, in the same algorithms, reaching the desired output.

As an alternative, although not completely native in R, Saga Gis 7.4 just implemented full 3d kriging. Saga can be accessed from R using the RSAGA package.