[R-sig-Geo] Indicator kriging map creation

Matevž Pavlič matevz.pavlic at gi-zrmk.si
Sun Jul 24 20:52:59 CEST 2011 

Hi Edzer, 

thanks for the help. It works ok on the sampleI provided, even if I enlarge the number of classes to 10. But when i try to use it on my dataset which has identical names of predictions and number of predisctions and covarinaces I get an error :

Error in factor(apply(cbind(z$nic.pred, z$ena.pred, z$dva.pred, z$tri.pred,  : 
  invalid labels; length 10 should be 1 or 8

This is the oputput of the prediction whit

> str(z):
Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
  ..@ data       :'data.frame':	8736 obs. of  65 variables:
  .. ..$ devet.pred     : num [1:8736] 0.001153 0.000876 0.100325 0.101093 0.102047 ...
  .. ..$ devet.var      : num [1:8736] 0.00548 0.00548 0.00549 0.0055 0.00551 ...
  .. ..$ dva.pred       : num [1:8736] 0.182 0.185 0.177 0.194 0.191 ...
  .. ..$ dva.var        : num [1:8736] 0.0744 0.0744 0.0745 0.0752 0.0753 ...
  .. ..$ ena.pred       : num [1:8736] 0.2009 0.1895 0.0476 0.0522 0.0482 ...
  .. ..$ ena.var        : num [1:8736] 0.149 0.149 0.15 0.155 0.155 ...
  .. ..$ nic.pred       : num [1:8736] 0.139 0.142 0.206 0.0252 0.0234 ...
  .. ..$ nic.var        : num [1:8736] 0.258 0.259 0.26 0.267 0.268 ...
  .. ..$ osem.pred      : num [1:8736] 9.15e-05 0.00 5.56e-05 6.82e-04 1.13e-03 ...
  .. ..$ osem.var       : num [1:8736] 0.00381 0.00382 0.00382 0.00386 0.00386 ...
  .. ..$ pet.pred       : num [1:8736] 0.000115 0.000303 0.000875 0.000423 0.000345 ...
  .. ..$ pet.var        : num [1:8736] 0.0126 0.0126 0.0126 0.0133 0.0133 ...
  .. ..$ sedem.pred     : num [1:8736] 0.0934 0.0924 0.0756 0.1084 0.1157 ...
  .. ..$ sedem.var      : num [1:8736] 0.0848 0.085 0.0852 0.0874 0.0877 ...
  .. ..$ set.pred       : num [1:8736] 0.271 0.278 0.27 0.395 0.389 ...
  .. ..$ set.var        : num [1:8736] 0.0881 0.0882 0.0885 0.0902 0.0904 ...
  .. ..$ stiri.pred     : num [1:8736] 0.001931 0.000777 0.0046 0 0 ...
  .. ..$ stiri.var      : num [1:8736] 0.101 0.101 0.101 0.106 0.106 ...
  .. ..$ tri.pred       : num [1:8736] 0.11 0.111 0.118 0.123 0.129 ...
  .. ..$ tri.var        : num [1:8736] 0.134 0.134 0.134 0.136 0.136 ...

I would like to create a map from predictions of course. 

What i found out is that your latest code works until labels line gets called.But funny thing is , when  I use this line :
z$class = factor(apply(cbind(z$nic.pred,z$ena.pred, z$dva.pred, z$tri.pred, z$stiri.pred,   z$pet.pred, z$set.pred, z$sedem.pred, z$osem.pred, z$devet.pred), 1, which.max))

i get :

> levels(z$class)
[1] "1" "2" "3" "4" "5" "6" "7" "8"

Any ideas? 

tnx, m
-----Original Message-----
From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-project.org] On Behalf Of Edzer Pebesma
Sent: Sunday, July 24, 2011 8:04 PM
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] Indicator kriging map creation

Matevz, try

z$class = factor(apply(cbind(z$dva.pred, z$ena.pred, z$nic.pred,
	z$stiri.pred, z$tri.pred), 1, which.max),
	labels=c("dva", "ena", "nic", "stiri", "tri"))
spplot(z["class"])


On 07/24/2011 06:22 PM, Matevž Pavlič wrote:
> 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
> 
> 
> 
> 
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

More information about the R-sig-Geo mailing list