Hi list users,
i am trying to do indicator kriging on a set of soild samples. I manage to do the predction of several classes into a "SpatialPixelsDataFrame" using the gstat object and predict() function. First a created a gstat object of several soil classes which I later fitted with a linear model of coregonalization. After that i used predict() to generate a SpatialPixelsDataFrame which hold information about predictions, variance and covariances. I am mostly interested in predoctions of each class. I know how to create a spplot or levelplot of each of the predicted classes on its own, but what i would like to achieve is to create a map of this predictions, so that the class with highest value (from 0-1) would prevail on each of the grid cell. In this way i would get a map of soil classes. Attached is the reproducible sample and a picture of prediction that i would like covert to a map of Soil classes (SoilFactor). But i get stuck when i try to create a actual map. That means getting the max value of each pixels and assign a name to that pixel. I hope i explained ok, if not i'll try again. Thanks for any ideas or help, m
<<attachment: temp.png>>
library(lattice) library(rgdal) library(gstat) #####CREATE DATASET n<-100 DF <- data.frame(X = rnorm(n, mean=100, sd=200), Y = rnorm(n, mean=500, sd=500), SoilFactor = c("nic", "ena","dva","tri", "stiri")) print(xyplot(X~Y, data=DF, panel=function(x, y,...){ ltext(x, y, labels=DF[,"SoilFactor"], cex=0.7, col="red")})) ## create table Soil.table<-table(DF [, "SoilFactor"]) Soil.table tmp<-as.vector(DF[,"SoilFactor"]) DF[,"SoilFactor"] <-factor(tmp) lvls<-levels(DF[,"SoilFactor"]) print(lvls) ### CREATE GSTAT OBJECT -01 classes nmx<-10 dat.c <- gstat(id = lvls[1], formula = I(SoilFactor == lvls[1]) ~ 1, loc=~X+Y, data = DF, nmax = nmx, set = list(order = 2, zero = 0.01)) dat.c <- gstat(dat.c, lvls[2], I(SoilFactor == lvls[2]) ~ 1, ~X + Y, DF, nmax = nmx) dat.c <- gstat(dat.c, lvls[3], I(SoilFactor == lvls[3]) ~ 1, ~X + Y, DF, nmax = nmx) dat.c <- gstat(dat.c, lvls[4], I(SoilFactor == lvls[4]) ~ 1, ~X + Y, DF, nmax = nmx) dat.c <- gstat(dat.c, lvls[5], I(SoilFactor == lvls[5]) ~ 1, ~X + Y, DF, nmax = nmx) #fit variogram plot(variogram(dat.c)) dat.c <- gstat(dat.c, model = vgm(.2, "Cir", 400, 0), fill.all = T) x <- variogram(dat.c) # fit lmc dat.fit = fit.lmc(x, dat.c, fit.ranges=F, fit.sills=c(T,T,T)) plot(x, model = dat.fit) ####CREATE PREDICTION GRID GridRes<-10 coordinates(DF)<-~X+Y x.range <- as.integer(range(DF@coords[,1])) y.range <- as.integer(range(DF@coords[,2])) grd<-expand.grid(x=seq(from=x.range[1], to=x.range[2], by=GridRes), y=seq(from=y.range[1], to=y.range[2], by=GridRes)) coordinates(grd)<-~y+x gridded(grd)<-T plot(grd) #CHECK THAT GRID AND POINTS HAVE SAME COORDINATE NAMES dimnames(coordinates(DF)) dimnames(coordinates(grd)) dimnames(grd@coords) <-list(NULL,c("X", "Y")) dimnames(grd@bbox) <-list(c("X", "Y")) #PREDICT z<-predict(dat.fit, newdata=grd) str(z) print(spplot(z, z=seq(1,10, by=2), main="Prediction", col.regions=heat.colors(64)), split=c(1,1,1,2), more=T) print(spplot(z, z=seq(2,19, by=2), main="Variance",col.regions=heat.colors(64)), split=c(1,2,1,2))
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo