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 -------------- next part -------------- An HTML attachment was scrubbed... URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20110724/d65f6e00/attachment.html> -------------- next part -------------- A non-text attachment was scrubbed... Name: temp.png Type: image/png Size: 16254 bytes Desc: temp.png URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20110724/d65f6e00/attachment.png> -------------- next part -------------- An embedded and charset-unspecified text was scrubbed... Name: temp.txt URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20110724/d65f6e00/attachment.txt>
Indicator kriging map creation
7 messages · Matevz Pavlic, Edzer Pebesma
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
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
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
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? 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? 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
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
Hi Edzer, sorry for the long time to reply.
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.
What do you think would be the best thing to do about that? Predict the classes on 3D grid? I mean, i would like to create a code that would generate maps of several classes (variing numbers of them) each depth interval. I think that prediction on 3D grid would be the right thing?
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.
Attached are the co-variograms of ten classes. As i can see there are few classes (devet=level 10 , osem=level 9...) that have nugget (pure nugget effect). This is probably because there are only a few point with that class in the area :
devet dva ena nic osem pet sedem set stiri tri
5 50 85 209 3 7 83 71 49 89
Is there a way to fix this without discarding this classes?
Regards,
m
-------------- next part --------------
A non-text attachment was scrubbed...
Name: variogram.png
Type: image/png
Size: 26318 bytes
Desc: variogram.png
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20110725/487b5c82/attachment.png>