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)?