[R-sig-Geo] Three-dimensional interpolation of temperature profiles using gstat?

Wed Jun 4 09:20:38 CEST 2014

On 06/04/2014 03:00 AM, Trevor Joyce wrote:
> 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?

You seem to use an anisotropic variogram model, meaning equal decay of
autocorrelation in all directions. Even if your dimensions had
comparable units (e.g., meters), it is very unlikely that this would be
satisfactory, as nearly everything on earth changes faster in the
vertical as in horizontal. Also, most phenomena will show a trend in the
vertical (deeper = colder?) that kriging does not take into account when
it is not being told to do so.

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

that should not be a problem.

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

If you replace the word "accurate" with "meaningful": yes. After that,
you can start thinking about anisotropy in a meaningful way (i.e.
unit-less), I would say.

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

I would start thinking about a trend first, then about anisotropy of the
residual process.

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
83 33081 http://ifgi.uni-muenster.de GPG key ID 0xAC227795

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 490 bytes
Desc: OpenPGP digital signature
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140604/a25282fb/attachment.bin>