data(quakes)
set.seed(2)
quakes <- quakes[ ss <- sample(1:nrow(quakes), 50), ]
# Artificially split data to create prediction data set
quakes.pred <- quakes[ -ss, ]
quakes.geogam <- geoGAM(response = "mag",
covariates = c("depth", "stations"),
data = quakes,
max.stop = 5,
cores = 1)
predicted <- predict(quakes.geogam, newdata = quakes.pred, type = "response" )
# \donttest{
## Use soil data set of soil mapping study area near Berne
data(berne)
data(berne.grid)
# Split data sets and
# remove rows with missing values in response and covariates
d.cal <- berne[ berne$dataset == "calibration" & complete.cases(berne), ]
### Model selection for numeric response
ph10.geogam <- geoGAM(response = "ph.0.10",
covariates = names(d.cal)[14:ncol(d.cal)],
coords = c("x", "y"),
data = d.cal,
seed = 1,
cores = 1)
# Create GRID output with predictions
sp.grid <- berne.grid[, c("x", "y")]
sp.grid$pred.ph.0.10 <- predict(ph10.geogam, newdata = berne.grid)
if(requireNamespace("raster")){
require("sp")
# transform to sp object
coordinates(sp.grid) <- ~ x + y
# assign Swiss CH1903 / LV03 projection
proj4string(sp.grid) <- CRS("EPSG:21781")
# transform to grid
gridded(sp.grid) <- TRUE
plot(sp.grid)
# optionally save result to GeoTiff
# writeRaster(raster(sp.grid, layer = "pred.ph.0.10"),
# filename= "raspH10.tif", datatype = "FLT4S", format ="GTiff")
}
# }
Run the code above in your browser using DataLab