0

I try to do a kriging in the Jakarta Bay. I have a set of measurement points with appropriated coordinates and attributes (pH, salinity,...)

In order to do a kriging I first need to find a model for my variogram. When I use the "variogram" function the output is not perfect but is should be ok, but then when I try to fit the variogram I get a waring message: In fit.variogram(ph.vgm, model = vgm(0.12, "Sph", 0.1, 0.01)) : Warning: singular model in variogram fit and I have a singular model.

Here I read about singular models associated to variogram calculations. Can I do something to make it better ?

How can I obtain a better fit for my variogram ? Why do I obtain only little circles around the measurement points ? I would like to have my full map with prediction values.

I also tried the "automap" library which is even less flexible and I don't obtain good results.

 library(sp)
library(gstat)
library(automap)

x = c(11878417.51,11882987.17,11887690.42,11892582.91,11897119.18,11902527.08,11879348.14,11884237.29,11888933.86,11893819.67,11898835.73,11903940.84,11908386.94,11885529.71,11889836.66,11900118.13,11905765.37,11896037.16,11901234.67,11906244.04,11892136.86,11900822.56,11904493.1,11907692.42,11910346.05,11888709,11887268.41,11885237.28,11883450.38,11880668.5)
y = c(-668537.7429,-667290.838,-666043.9586,-663943.1247,-663992.3709,-662612.3726,-672878.6036,-672014.4364,-671960.7062,-669604.4601,-668541.1009,-667203.5333,-666181.6289,-676933.1896,-676566.0044,-673095.7667,-671736.8309,-679340.0992,-677788.4711,-676606.3051,-682542.446,-680607.5158,-680131.1539,-679733.0503,-662307.2774,-680754.1755,-681408.3272,-680494.7783,-680491.4197,-679426.19)
ph = c(7.1,7.76,7.14,7.19,7.56,7.56,7.11,8.14,7.22,7.17,7.33,7.37,7.36,7.23,7.12,7.54,7.96,7.98,7.96,7.2,7.44,7.36,7.71,7.71,8.01,7.73,8.11,7.03,7.26,7.77)
TSS = c(13.7,21,17.7,18.8,4.7,12.4,17.3,18.8,20.2,18.3,5.6,NA,NA,NA,21.9,11.1,NA,NA,21.2,29.1,31.3,29.3,21.3,25.4,31.8,14.5,2.9,11.7,8.4,NA)

df = data.frame(x,y,ph,TSS)
coordinates(df) = ~x+y
proj4string(df) <- CRS("+init=epsg:3857")
spplot(df)

grid <- data.frame(x=c(11877328.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11909828.43,11912328.43,11877328.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11877328.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11909828.43,11912328.43,11877328.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11909828.43,11912328.43,11877328.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11909828.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11909828.43,11912328.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11879828.43,11882328.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11902328.43,11904828.43,11907328.43,11909828.43,11884828.43,11887328.43,11889828.43,11892328.43,11894828.43,11897328.43,11899828.43,11904828.43,11894828.43),y=c(-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-659719.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-662219.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-664719.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-667219.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-669719.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-672219.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-674719.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-677219.4518,-679719.4518,-679719.4518,-679719.4518,-679719.4518,-679719.4518,-679719.4518,-679719.4518,-679719.4518,-682219.4518))

coordinates(grid)=~x+y
proj4string(grid) <- CRS("+init=epsg:3857")
gridded(grid)=T
spplot(grid)

ph.vgm <- variogram(ph~1, df[!is.na(df@data$ph),]); plot(ph.vgm)
ph.fit = fit.variogram(ph.vgm, model = vgm(0.12, "Sph", 4000, 0.01), warn.if.neg = T); ph.fit
plot(ph.vgm, ph.fit)
ph.kriged = krige(ph~1, df[!is.na(df@data$ph),], grid, model = ph.fit)
spplot(ph.kriged["var1.pred"])
ZKB
  • 101
  • 2

1 Answers1

0

And small sample.

If you want to see $0.5(Z(x)-Z(x+h))^2$ for every point pair, plotted against $h$, then use

plot(variogram(ph~1, df[!is.na(df@data$ph),], cloud=TRUE))

but this may also not make you happy. The bottomline is that you have very few observations (30), and with so few observations you essentially never get variograms that look nice.

Edzer Pebesma
  • 3,814
  • 16
  • 26