[R-sig-Geo] Simulating spatial data using gstat

[Surya Kr]

Hi All -

I'm experimenting with the gstat package trying to do something basic. I
want to simulate data with an isotropic covariance structure and turn
around and fit a variogram to the simulated data. I get several warnings
while estimating the variogram indicating failed convergence -

Warning messages:
1: In fit.variogram(o, m, fit.kappa = FALSE, fit.method = fit.method,  ... :
   No convergence after 200 iterations: try different initial values?

My trend surface has the following linear term:
longitude+latitude+covariate_x. Note there is no intercept. The error is
from the Matern model. I'm new to spatial data analysis in R. Very much
appreciate your patience and help if you could correct me if I missed out
on anything. Here's my code -

### simulate the data
# Matern with partial sill, range, and nugget.
v.sim <- vgm(200, "Mat", 30, 20, kappa = 0.1)
g.sim <- gstat(NULL, "z", z~longitude+latitude+covariate_x,
locations=~longitude_dgr+latitude_dgr, beta = c(0.10, 0.20, 0.40),
dummy=TRUE, model = v.sim, nmax=50)
yy <- predict(g.sim, newdata = mydata, nsim = 1) # mydata is a data.table
and has the longitude, latitude, and covariate_x
coordinates(yy) <- c("longitude", "latitude")

### Do an empirical variogram and fit a Matern model
## Add covariate_x to the simulated dataset
yy$covariate_x = mydata$covariate_x
variog_gstat <- variogram(sim1 ~ longitude+latitude+covariate_x, data = yy)
fit.g <- fit.variogram(variog_gstat, vgm(200, "Mat", 30, 20),
fit.kappa=TRUE)
> fit.g
   model    psill    range kappa
1   Nug  95.4979  0.00000   0.0
2   Mat 210.2117 46.92647   0.3
#plot(variog_gstat, fit.g)
## confirming that the above is equivalent to analyzing residuals from
ordinary least squares
lmObj <- lm(sim1 ~ longitude + latitude + covariate_x, data=yy)
yy$lmresids <- lmObj$residuals
variog_gstat1 <- variogram(lmresids ~ 1, data = yy)
#plot(variog_gstat1)
fit.g1 <- fit.variogram(variog_gstat1, vgm(200, "Mat", 30, 20),
fit.kappa=TRUE)
> fit.g1
   model    psill    range kappa
1   Nug  95.4979  0.00000   0.0
2   Mat 210.2117 46.92647   0.3
#plot(variog_gstat1, fit.g1)

My questions:
------------------
1. Did I simulate the dataset with the spatial trend and covariance right?
2. How to correctly estimate variogram from empirical variogram without
getting those warnings. Note that my initial values for the parameters are
the same as the true values. How to make the variogram fitting exercise
more reliable/robust?
2.1. How to choose initial values in an automated fashion? I need to fit
variograms for several spatial datasets (possibly in parallel on a cluster)
and I need to automate choosing the initial parameter values.
3. Is there an example showing similar code for variogram fitting and
simulation of spatial data using the "Fields" package?

Thanks much,
Surya

ps: can I attached plots when I message R-sig-geo and have them displayed
on the message?