[R-sig-Geo] Indicator kriging map creation

Edzer Pebesma edzer.pebesma at uni-muenster.de
Sun Jul 24 21:30:45 CEST 2011 

You have to dive a little bit into what factor() does. One or more of
your classes never has the highest probability, it seems:

> factor(c(1,2,3,5))
[1] 1 2 3 5
Levels: 1 2 3 5
> factor(c(1,2,3,5), labels = c("1", "2", "3", "4", "5"))
Error in factor(c(1, 2, 3, 5), labels = c("1", "2", "3", "4", "5")) :
  invalid labels; length 5 should be 1 or 4
> levels(factor(c(1,2,3,5)))
[1] "1" "2" "3" "5"
> factor(c(1,2,3,5), levels=1:5, labels = c("1", "2", "3", "4", "5"))
[1] 1 2 3 5
Levels: 1 2 3 4 5
> factor(c(1,2,3,5), levels=1:5, labels = c("a", "b", "c", "d", "e"))
[1] a b c e
Levels: a b c d e


On 07/24/2011 08:52 PM, Matevž Pavlič wrote:
> 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

-- 
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

More information about the R-sig-Geo mailing list