Hi,

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 based on 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 example 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)?


Many thanks!

Trevor

Ph.D. Student
Scripps Institution of Oceanography
9500 Gilman Dr. MC 0208
La Jolla, CA 92093

	[[alternative HTML version deleted]]

