[R-sig-Geo] Indicator kriging map creation

Edzer Pebesma edzer.pebesma at uni-muenster.de
Sun Jul 24 23:28:38 CEST 2011 


On 07/24/2011 10:30 PM, Matevž Pavlič wrote:
> Hi, 
> 
> Thanks for the help! I realize I have a lot to learn , but sometimes you just get stuck. 
> 
> What i don't understad is this :
> when I use the first line of your code :
>> 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))
> 
> and check levels :
>> levels(z$class) 
> [1] "1" "2" "3" "4" "5" "6" "7" "8"
> 
> the levels are numerical.  I read in the help that " optional vector of the values that x might have taken. The default is the unique set of values taken by as.character(x), sorted into increasing order of x". 
> 
> I would like to create several IK maps with depth (every meter of depth). With doing that there will be areas that will not include all of the soil classes (or in this case factor levels). Will this be a problem with factor labels and map creating, if some of the soil classes will be missing? 

It might be annoying, and after all, it is your sole responsibility to
make sure the right labels end up in the right place. R can only prevent
you from making mistakes that are too obvious.

> 
> Also, you wrote " One or more of your classes never has the highest probability, it seems ". 
> As i understand the apply functions it checks the rows in grid and assignes the maximum values of pred to the cell. So  i do not understand how some predictions never have the highest probability since i have prediction in a SpatialGridDataFrame for each of the 10 classes? Could it be because the values in some cells are equal? 

Not so likely; I think it comes with simple kriging -- rare classes
often show weak autocorrelation, resulting in low estimates everywhere.
Without having seen your data, I suspect your classes 9 and 10 are like
that.

Please do read ?which.max

> 
> Thanks again for the answer and your time, 
> 
> matevz
>  
> 
> -----Original Message-----
> From: Edzer Pebesma [mailto:edzer.pebesma at uni-muenster.de] 
> Sent: Sunday, July 24, 2011 9:31 PM
> To: Matevž Pavlič
> Cc: r-sig-geo at r-project.org
> Subject: Re: [R-sig-Geo] Indicator kriging map creation
> 
> 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

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