Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'm trying to interpolate the values of an attribute (called "disease_rate") that i've sampled at a set of spatial points and i'd like to add 'weights' to the resulting kriging estimates.

Basically, I have about 1500 points (a random sample of 50 is shown below). I have an attribute vector (the rate of a particular disease) and a vector of absolute population sizes ("population"). I'd like to interpolate about 10,000 new points for "disease_rate", and at the same time weight the estimates by absolute population size (as areas with higher absolute populations have absolutely more people with the disease, and should probably have greater influence on the kriging estimates).

I'm using the packages "automap" and "gstat", with the following code:

library(automap)
vgm <- autofitVariogram(disease_rate ~ 1, , input_data = dat)
plot(vgm)

kr.ok <- autoKrige(disease_rate ~ 1, input_data = dat, model = "Exp",
    verbose = TRUE, maxdist = Inf) 
plot(kr.ok)

I'm stuck on two points:

1) how to add the weights? Is there an argument for this within the autoKrige() or krige() functions? I don't see a way to do it.

2) when using autokrige() from the package "automap" a set of 5000 point locations for interpolated values is generated by default, by using I believe a convex hull. In "gstat", how do I get a similar sample of point locations? At the moment I'm doing this:

prediction_locations <- spsample(SpatialPolygonsDataFrame, n = 10000,
    type = "regular")
prediction_locations@coords <- jitter(prediction_locations@coords)

the object "SpatialPolygonsDataFrame" is too big to include here, but provides a reference area within which to sample the points. However, I always get the following error:

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

Any help/suggestions for the above two problems is much appreciated.

# sample of data (50 points)
dat <- new("SpatialPointsDataFrame"
, data = structure(list(population = structure(c(13L, 37L, 9L, 25L, 5L, 
2L, 40L, 31L, 10L, 43L, 22L, 7L, 32L, 19L, 12L, 28L, 23L, 30L, 
44L, 6L, 39L, 1L, 4L, 35L, 15L, 16L, 20L, 21L, 27L, 24L, 18L, 
41L, 42L, 26L, 14L, 45L, 47L, 33L, 36L, 3L, 46L, 4L, 4L, 11L, 
29L, 38L, 48L, 34L, 8L, 17L), .Label = c("10323", "1070", "1142", 
"1218", "1226", "1228", "12858", "1583", "1619", "1645", "1737", 
"1832", "1944", "2003", "2125", "2140", "2248", "273", "2816", 
"2868", "2986", "3093", "322", "3338", "3707", "371", "3841", 
"4151", "4184", "4189", "4198", "434", "4528", "491", "546", 
"560", "65", "654", "7138", "735", "776", "82", "845", "879", 
"91027", "943", "948", "989"), class = "factor"), disease_rate = c(0, 
0, 154.4163064, 140.2751551, 228.3849918, 0, 0, 71.46260124, 
0, 0, 226.3174911, 107.5091726, 0, 159.8011364, 0, 120.4529029, 
0, 45.35688709, 227.5312856, 0, 77.59122357, 33.5322916, 0, 0, 
0, 93.45794393, 239.0914525, 0, 0, 0, 0, 0, 0, 0, 0, 189.7113215, 
0, 0, 0, 175.1313485, 0, 369.4581281, 0, 0, 103.5691523, 0, 0, 
0, 0, 111.2099644)), .Names = c("population", "disease_rate"), class = "data.frame",     row.names = c(351L, 
251L, 477L, 349L, 56L, 434L, 299L, 248L, 301L, 262L, 115L, 325L, 
402L, 238L, 73L, 376L, 454L, 470L, 6L, 359L, 473L, 233L, 448L, 
 121L, 35L, 253L, 324L, 369L, 492L, 463L, 400L, 365L, 245L, 157L, 
264L, 101L, 367L, 142L, 109L, 389L, 158L, 254L, 457L, 68L, 219L, 
395L, 450L, 380L, 474L, 199L))
, coords.nrs = c(4L, 3L)
, coords = structure(c(634873.420766408, 541790.625187116, 649805.450219927, 
536333.978570239, 581957.093009378, 607202.361088042, 617240.615621174, 
505627.220834258, 620543.999367496, 592414.016643991, 579995.994396474, 
636851.153664358, 593891.033349101, 594650.473384273, 607103.681876253, 
659509.454474378, 604222.534339966, 592766.490473358, 625515.430661865, 
545043.782766906, 663659.275576225, 626566.467170621, 687406.871345609, 
579464.084578232, 594577.343169983, 650081.936891944, 661853.556699763, 
620028.294881262, 622881.614147092, 714608.702621971, 620991.249432193, 
531188.877201674, 654835.610827768, 527241.805893339, 604428.976031172, 
688457.08368477, 613888.674693582, 616783.067305461, 500082.617682936, 
609515.575391672, 611609.283708876, 689412.437797476, 675561.054104031, 
601548.525491549, 656484.375080421, 677341.365311871, 616862.413097344, 
644065.67844776, 475677.631939326, 617193.359923705, 5144526.77927643, 
5328308.33322789, 5216468.95632244, 5278378.21530125, 5242407.33311755, 
5107204.72128753, 5104424.13808425, 5346313.09064684, 5291125.71228228, 
5260872.67506553, 5296414.86191041, 5153372.04174738, 5107141.20462901, 
5126046.0142604, 5188835.84829657, 5241495.91305329, 5120799.1233737, 
5077612.68066607, 5173630.04242822, 5351027.28834978, 5231260.93265272, 
5321245.1139323, 5239845.72123907, 5257533.19157413, 5093198.30043019, 
5176787.1924563, 5208692.02318177, 5301924.71404294, 5307513.36239273, 
5170246.64046221, 5129845.77701215, 5284862.55527842, 5147875.41139468, 
5295710.09388113, 5293345.39302848, 5189299.58825685, 5165031.11113162, 
5136185.92439977, 5334176.39641511, 5166895.32812857, 5279273.17144778, 
5236974.38111571, 5222457.79639944, 5279521.25697787, 5171514.42644342, 
5244570.31251435, 5113000.12179841, 5170833.89650194, 5346485.72226854, 
5068270.63886515), .Dim = c(50L, 2L), .Dimnames = list(NULL, 
c("long.proj", "lat.proj")))
, bbox = structure(c(475677.631939326, 5068270.63886515, 714608.702621971, 
5351027.28834978), .Dim = c(2L, 2L), .Dimnames = list(c("long.proj", 
"lat.proj"), c("min", "max")))
, proj4string = new("CRS"
, projargs = NA_character_
)
)
share|improve this question
1  
There seem to be some misconceptions here. First, why should population have any correspondence at all with disease rate? Second, the BLUP criteria for kriging determine the weights. If you were to force different weights, you wouldn't be doing kriging. Thus it makes no sense for kriging to accept weights as an argument. It sounds like you might need a slightly more sophisticated approach, such as a kriging-based GLM like that offered in GeoRGLM –  whuber Mar 21 '12 at 16:44
add comment

1 Answer 1

Maybe it is too late but this error:

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

It is sometimes related with some points at nearly the same position.
Did you try zerodist(dat) to check that?

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.