[R-sig-Geo] GSTAT Co-Kriging: How can i fix a non positive definte correlation matrix

Wed Jan 7 14:37:26 CET 2015

Hi Blake

You should do the prediction using the 'g.fit' object instead of the 'g' object.

This works:

See below for a dput'ed version of the data 'dtrd'.

## R script

library(gstat)
library(sp)

## create a grid onto which i will interpolate # first get a range in data
bbox <- slot(dtrd, "bbox")
x.range <- bbox["x", ]
y.range <- bbox["y", ]

## now expand to a grid with desired spacing
grid <- expand.grid(x=seq(from=x.range[1], to=x.range[2], by=0.4), y=seq(from=y.range[1], to=y.range[2], by=1.6))
plot(grid)

# Convert to spatial data frame
coordinates(grid) = c("x", "y")

# Compute cross-variogram
g <- gstat(NULL, id = "n.ae", formula = n.ae ~ + 1, data = dtrd, nmax=10)
g <- gstat(g, "n.ts", n.ts ~ 1, dtrd, nmax=10)
g <- gstat(g, model = vgm(psill=-0.65, model = "Sph", range = 15, nugget = 0.05), fill.all = T)

## fit direct and cross variograms
g.fit <- fit.lmc(v, g, correct.diagonal = 1.01)

## plot direct and cross variograms
plot(variogram(g.fit), model=g.fit$model)

## Perform co-kriging predictions
ck <- predict.gstat(g.fit, grid) # Use g.fit and not g!!!!!!!!!

str(ck)

##
gridded(ck) <- TRUE
spplot(ck, c("n.ae.pred", "n.ts.pred"))

## data
dtrd <- new("SpatialPointsDataFrame"
    , data = structure(list(n.ts = c(-0.052514693, -0.375234492, 1.362810557, 
-0.333087778, 0.699916885, 0.605552409, 0.406868496, -0.529754386, 
-0.416509276, -0.450700276, -0.732937251, 0.377457633, -0.989029229, 
-1.259551904, -0.470220669, 1.046426438, -1.440639861, -0.899563555, 
-1.287567164, -0.816931437, -0.782655308, 0.048342822, -0.5581308, 
-0.217874377, -0.272702769, -0.996064, -1.97520507, -0.718163458, 
0.337266312, -0.565668608, 0.651736313, -0.715751462, 0.650820529, 
-1.184380204, 0.132388601, 1.08502868, -0.85196341, 0.099887288, 
-0.233091874, -0.989196908, -0.025054063, -0.394999955, 0.049839033, 
1.327518037, -0.421157849, -0.178717505, -0.065413063, 1.167410569, 
0.428968564, 1.550074254, 2.931334909, 4.126464345, 0.44561262, 
0.543591219, 0.716228164, 1.856183532, 0.774644882, 0.425336383, 
0.294908064, 2.38260471, -0.315391214, 1.136916243, -0.233091874, 
-0.352017426, -0.933099317, -0.865945243, -0.594339105, -0.653328511, 
0.076670222, -0.14022361, -0.545492977, -0.084386566, 2.107283841, 
-1.645006796, 0.512684145, 0.517794479, -0.359495901, 0.548288806, 
-0.597741695, -1.003803022, 0.060157729, 0.251793973, -1.148625922, 
-0.957779058, 0.267364885, 0.171656399, -0.912461925, -0.645228334, 
0.83148026, -1.187664129, -0.467591981, -0.121964677, -0.158281328, 
-0.235932095, -0.436274739, -1.136811015, 0.87349541, -0.661067533, 
0.301101862, 1.265698729, -0.806040053, 1.006815543, -0.0132082, 
-1.214201235, -0.323473333, -0.003250658, -0.168398809, -0.854793312, 
0.507014021, 3.104090518, -0.193785381, -0.71632415, 2.480632323, 
-0.534663506, -0.124544351), n.ae = c(0.055104755, -0.861381682, 
-0.23577909, 0.134741672, -2.348188189, -0.75492334, 0.329847607, 
0.251549504, 1.094807283, 0.988002955, 1.344413512, -0.542187169, 
2.002433474, -0.895859924, 0.370809323, 0.168121785, 1.423433672, 
0.83844681, -0.086478673, 1.156633451, 0.299641538, -0.763091614, 
-0.130509134, 0.815852429, 0.906621066, 0.274459785, 2.017626766, 
0.444173367, 0.831151021, 0.513114809, -0.055114304, 0.205488257, 
0.558032798, 0.612006591, 0.462059331, 0.399240335, 0.113079949, 
-1.922971577, -0.287556625, 0.10047403, 0.384754057, -1.272051846, 
-2.568551091, -0.974203156, -0.0101211, 0.408401438, -0.195765074, 
-1.626236059, -1.287922066, -1.011990828, -1.872352342, -2.254471152, 
-0.709554064, 1.03999109, -0.671315106, -1.403466298, -0.952601604, 
-0.776073605, -0.00019281, -0.778811406, -0.223383771, -0.238938091, 
-0.520269718, -0.556041646, -0.441354858, -1.011088256, -0.399550741, 
1.213119402, 0.819568017, 0.166527242, -0.613535469, -0.423754708, 
-0.891753223, 0.566968259, 0.515250895, -1.144578624, 0.191964723, 
-1.011825356, 0.046560408, 0.319904274, 0.824035747, 0.632946256, 
1.257360463, 0.30201831, 1.139304073, 0.623288738, 0.820049388, 
1.041375034, 2.914000894, 2.756035789, 0.467880919, 0.878460826, 
-0.266451488, 1.440612621, -0.298237058, -0.688794913, 0.449964869, 
-0.08882536, -0.255259598, 0.948786211, 0.352201302, 1.12999254, 
0.581935908, 1.431842632, -0.083379843, -0.248324838, -0.916859761, 
-0.231928117, -0.49107152, -2.654897126, 0.067244345, -1.46047875, 
-0.336852087, -1.093553232, -0.603005465)), .Names = c("n.ts", 
"n.ae"), class = "data.frame", row.names = c(NA, -115L))
    , coords.nrs = 1:2
    , coords = structure(c(3.2, 6.4, 9.6, 12.8, 16, 22.4, 25.6, 3.2, 6.4, 9.6, 
12.8, 16, 19.2, 22.4, 25.6, 3.2, 6.4, 9.6, 12.8, 16, 19.2, 22.4, 
3.2, 6.4, 9.6, 12.8, 16, 19.2, 22.4, 25.6, 3.2, 6.4, 9.6, 12.8, 
16, 19.2, 22.4, 25.6, 3.2, 6.4, 9.6, 12.8, 16, 19.2, 22.4, 25.6, 
3.2, 6.4, 9.6, 12.8, 16, 19.2, 22.4, 25.6, 3.2, 6.4, 9.6, 12.8, 
16, 19.2, 22.4, 25.6, 3.2, 6.4, 9.6, 16, 19.2, 22.4, 25.6, 3.2, 
6.4, 9.6, 12.8, 16, 22.4, 25.6, 3.2, 6.4, 9.6, 12.8, 16, 19.2, 
22.4, 25.6, 3.2, 6.4, 9.6, 12.8, 16, 19.2, 22.4, 25.6, 3.2, 6.4, 
12.8, 16, 19.2, 22.4, 25.6, 3.2, 6.4, 9.6, 12.8, 16, 19.2, 22.4, 
25.6, 3.2, 6.4, 9.6, 12.8, 16, 19.2, 22.4, 25.6, 6.4, 6.4, 6.4, 
6.4, 6.4, 6.4, 6.4, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 12.8, 
12.8, 12.8, 12.8, 12.8, 12.8, 12.8, 16, 16, 16, 16, 16, 16, 16, 
16, 19.2, 19.2, 19.2, 19.2, 19.2, 19.2, 19.2, 19.2, 22.4, 22.4, 
22.4, 22.4, 22.4, 22.4, 22.4, 22.4, 25.6, 25.6, 25.6, 25.6, 25.6, 
25.6, 25.6, 25.6, 28.8, 28.8, 28.8, 28.8, 28.8, 28.8, 28.8, 28.8, 
32, 32, 32, 32, 32, 32, 32, 35.2, 35.2, 35.2, 35.2, 35.2, 35.2, 
35.2, 38.4, 38.4, 38.4, 38.4, 38.4, 38.4, 38.4, 38.4, 41.6, 41.6, 
41.6, 41.6, 41.6, 41.6, 41.6, 41.6, 44.8, 44.8, 44.8, 44.8, 44.8, 
44.8, 44.8, 48, 48, 48, 48, 48, 48, 48, 48, 51.2, 51.2, 51.2, 
51.2, 51.2, 51.2, 51.2, 51.2), .Dim = c(115L, 2L), .Dimnames = list(
    NULL, c("x", "y")))
    , bbox = structure(c(3.2, 6.4, 25.6, 51.2), .Dim = c(2L, 2L), .Dimnames = list(
    c("x", "y"), c("min", "max")))
    , proj4string = new("CRS"
    , projargs = NA_character_
)
)

Yours sincerely / Med venlig hilsen

Frede Aakmann Tøgersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender. 

< snipped >