library(sp)
data(meuse, package= "sp")
coordinates(meuse)= ~x + y
proj4string(meuse) <- CRS("+init=epsg:28992")

library(rgdal)
writeOGR(meuse, ".", "meuse_points", driver = "ESRI Shapefile")

library(gstat)
vgm1 <- variogram(log(zinc)~1, meuse)
fit.variogram(vgm1, vgm("Sph"))       # Nugget= 0.214, Range= 514.4
plot(vgm1)                            # Sill ~ 0.85 by visual