# Simulating spatial data using gstat Classic List Threaded 1 message Reply | Threaded
Open this post in threaded view
|

## Simulating spatial data using gstat

 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? _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo