On 8/15/19 12:31 PM, Emanuele Barca wrote: > Dear Edzer, > > sorry for bothering you once more but I need to be sure about my R script. > > In summary, i'm comparing the performance of regression-kriging and > collocated co-kriging. > > Regression-kriging was based on an unique covariate, the elevation Z. > > I use Z as unique ancillary variable in the Co-kriging. > > As first attempt, the final raster maps were completely different. It > appeared that it was due > > to the fact that the dataset was non-stationary and only > regression-kriging overcomes this issue, while co-kriging not. > > But if I pass to universal co-kriging introducing Z as covariate, it > bacomes useless as ancillary variable! > > What is my mistake?
You assume that someone else can be sure about your analysis by looking at your R script without having access to your data. And you are starting a personal dialogue on a list, essentially uninviting everyone else to get involved. > > emanuele > > > > > Il 2019-08-15 12:00 r-sig-geo-requ...@r-project.org ha scritto: >> Send R-sig-Geo mailing list submissions to >> r-sig-geo@r-project.org >> >> To subscribe or unsubscribe via the World Wide Web, visit >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> or, via email, send a message with subject or body 'help' to >> r-sig-geo-requ...@r-project.org >> >> You can reach the person managing the list at >> r-sig-geo-ow...@r-project.org >> >> When replying, please edit your Subject line so it is more specific >> than "Re: Contents of R-sig-Geo digest..." >> >> >> Today's Topics: >> >> 1. Error running codes (Enoch Gyamfi Ampadu) >> 2. Re: regression-kriging and co-kriging (Edzer Pebesma) >> >> ---------------------------------------------------------------------- >> >> Message: 1 >> Date: Thu, 15 Aug 2019 06:51:49 +0200 >> From: Enoch Gyamfi Ampadu <egamp...@gmail.com> >> To: r-sig-geo@r-project.org >> Subject: [R-sig-Geo] Error running codes >> Message-ID: >> <cac4+n+yravrhc_mcoatv-ao_-sh7qyo8zukzaojqbnabgjx...@mail.gmail.com> >> Content-Type: text/plain; charset="utf-8" >> >> Dear List, >> >> Please I have been running a certain the codes below to read my image, >> shapeffile to enable me classify for cover with Random forest. I have >> gotten to the point of extracting the raster values and the raster >> information. However I keep getting errors. though I have tried different >> combination and other codes like readOGR for importing the training >> polygon. I will be glad if you could be of help as I am new to using R. >> >> #import clip image of study area >> >> TestImg3 <- brick("C:\\Users\\hp\\Desktop\\Data collection\\Nkandla >> Images\\Landsat\\2019\\08052019\\Corrected\\New >> Folder\\Image08052019_Clip.tif") >> >> #Assign name of bands >> names(TestImg3) <- c(paste0("B", 1:7, coll=""), "B9") >> plotRGB(TestImg3, r=4, g=3, b=2, stretch ="lin") >> >> #Import training shapefile >> sample <- shapefile("C:/Users/hp/Desktop/Data collection/Nkandla >> Images/Landsat/2019/08052019/Corrected/Train08052019_2.shp") >> >> responseCol <- "class" >> >> #(I have tried options of changing "class" to "classname" as reflect >> actaul >> name assigned in ArcMap) >> >> # Overlap the sample polygons on the image (combine the class information >> with extracted values). >> >> pixels = data.frame(matrix(vector(), nrow = 0, ncol = >> length(names(img)) + >> 1)) >> for (i in 1:length(unique(sample[[responseCol]]))){ >> category <- unique(sample[[responseCol]])[i] >> categorymap <- sample[sample[[responseCol]] == category,] >> dataSet <- extract(img, categorymap) >> dataSet <- dataSet[!unlist(lapply(dataSet, is.null))] >> if(is(sample, "SpatialPointsDataFrame")){ >> dataSet <- cbind(dataSet, class = as.numeric(category)) >> pixeles <- rbind(pixeles, dataSet) >> } >> if(is(sample, "SpatialPolygonsDataFrame")){ >> dataSet <- lapply(dataSet, function(x){cbind(x, class = >> as.numeric(rep(category, >> >> nrow(x))))}) >> df <- do.call("rbind", dataSet) >> pixels <- rbind(pixeles, df) >> } >> >> } >> >> THIS IS THE ERROR I GET FROM RUNNING THE ABOVE CODES >> >>> for (i in 1:length(unique(samples[[responseCol]]))){ >> + category <- unique(samples[[responseCol]])[i] >> + categorymap <- samples[samples[[responseCol]] == category,] >> + dataSet <- extract(img, categorymap) >> + >> + if(is(sample, "SpatialPointsDataFrame")){ >> + dataSet <- cbind(dataSet, class = as.numeric(category)) >> + pixeles <- rbind(pixeles, dataSet) >> + } >> + if(is(sample, "SpatialPolygonsDataFrame")){ >> + dataSet <- lapply(dataSet, function(x){cbind(x, class = >> as.numeric(rep(category, >> + >> nrow(x))))}) >> + df <- do.call("rbind", dataSet) >> + pixels <- rbind(pixeles, df) >> + } >> + >> + } >> Error in y[i, ] : >> cannot get a slot ("Polygons") from an object of type "NULL" >> >> >> Please help me out. >> Thank you. >> >> Best regards, >> >> Enoch >> >> -- >> *Enoch Gyamfi - Ampadu* >> >> *Geography & Environmental Sciences* >> >> *College of Agriculture, Engineering & Science* >> >> *University of KwaZulu-Natal, Westville Campus* >> >> *Private Bag X54001* >> *Durban, South Africa **– 4000**.* >> *Phone: +27 835 828255* >> >> *email: egamp...@gmail.com <egamp...@gmail.com>* >> >> >> *skype: enoch.ampadu* >> *The highest evidence of nobility is self-control*. >> >> *A simple act of kindness creates an endless ripple*. >> >> [[alternative HTML version deleted]] >> >> >> >> >> ------------------------------ >> >> Message: 2 >> Date: Thu, 15 Aug 2019 08:33:59 +0200 >> From: Edzer Pebesma <edzer.pebe...@uni-muenster.de> >> To: r-sig-geo@r-project.org >> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging >> Message-ID: <bdae5226-7023-7e30-1927-5fb66ff41...@uni-muenster.de> >> Content-Type: text/plain; charset="utf-8" >> >> >> >> On 8/12/19 8:21 PM, Emanuele Barca wrote: >>> Dear Edzer, >>> >>> maybe I found the solution. I found this in the predict function help: >>> "When a non-stationary (i.e., non-constant) mean is used, both for >>> simulation and prediction purposes the variogram model defined should be >>> that of the residual process, and not that of the raw observations" >>> Since my data were, actually, non-stationary, I applied the universal >>> co-kriging instead usual co-kriging. >>> now the maps of regression-kring and co-kriging are actually similar s >>> expected. >>> did I understand correctly the quoted sentence? >> >> I think so, but hard to be sure given the information you provide. >> >>> >>> regards >>> >>> emanuele barca >>> ------------------------------ >>>> >>>> Message: 2 >>>> Date: Sat, 10 Aug 2019 10:41:38 +0200 >>>> From: Edzer Pebesma <edzer.pebe...@uni-muenster.de> >>>> To: r-sig-geo@r-project.org >>>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging >>>> Message-ID: <da76da51-e46b-8a8d-4952-7c3b85c19...@uni-muenster.de> >>>> Content-Type: text/plain; charset="utf-8" >>>> >>>> Hard to tell from your script. Maybe give a reproducible example? >>>> >>>> On 8/6/19 1:07 PM, Emanuele Barca wrote: >>>>> Dear r-sig-geo friends, >>>>> >>>>> I produced two maps garnered in the following way: >>>>> >>>>> # for regression-kriging >>>>> Piezo.map <-autoKrige(LivStat ~ Z, input_data = mydata.sp, new_data >>>>> = covariates, model = "Ste") >>>>> >>>>> Piezork.pred <- Piezo.map$krige_output$var1.pred >>>>> Piezork.coords <- Piezo.map$krige_output@coords >>>>> Piezork.out <- as.data.frame(cbind(Piezork.coords, Piezork.pred)) >>>>> colnames(Piezork.out)[1:2] <- c("X", "Y") >>>>> coordinates(Piezork.out) = ~ X + Y >>>>> gridded(Piezork.out) <- TRUE >>>>> >>>>> spplot(Piezork.out, main = list(label = "R-k Hydraulic head", cex = >>>>> 1.5)) >>>>> >>>>> #for co-kriging >>>>> g <- gstat(id = "Piezo", formula = LivStat ~ 1, data = mydata.sp, set >>>>> = list(nocheck = 1)) >>>>> g <- gstat(g, id = "Z", formula = Z ~ 1, data = mydata.sp, set = >>>>> list(nocheck = 1)) >>>>> >>>>> v.g <- variogram(g) >>>>> >>>>> #g <- gstat(g, id = "Piezo", model = vgm(150, "Mat", 1350, 0.0, kappa >>>>> = 1.9), fill.all = T)# >>>>> g <- gstat(g, id = "Piezo", model = vgm(0.7, "Ste", 1300, 18, kappa = >>>>> 1.9), fill.all = T)# >>>>> g.fit <- fit.lmc(v.g, g, fit.lmc = TRUE, correct.diagonal = 1.01) # >>>>> fit multivariable variogram model , fit.lmc = TRUE, correct.diagonal >>>>> = 1.01 >>>>> g.fit >>>>> plot(v.g, model = g.fit, main = "Fitted Variogram Models - Raw Data")# >>>>> #gridded(covariates) <- TRUE >>>>> g.cok <- predict(g.fit, newdata = covariates)#grid >>>>> >>>>> g.cok.pred <- g.cok@data$Piezo.pred >>>>> aaaa <- na.omit(g.cok.pred) >>>>> g.cok.coords <- g.cok@coords >>>>> g.cok.out <- as.data.frame(cbind(g.cok.coords, g.cok.pred)) >>>>> colnames(g.cok.out)[1:2] <- c("X", "Y") >>>>> coordinates(g.cok.out) = ~ X + Y >>>>> gridded(g.cok.out) <- TRUE >>>>> spplot(g.cok.out, main = list(label = "Hydraulic head with >>>>> Co-kriging", cex = 1.5)) >>>>> >>>>> ########################################################################################################################### >>>>> >>>>> >>>>> >>>>> I am unable to understand why the first map appears as a raster and >>>>> the second not, notwithstanding the fact that they are both computed >>>>> on the same "covariates" DEM??? >>>>> >>>>> where is the mistake??? >>>>> >>>>> regards >>>>> >>>>> emanuele >>>>> >>>>> ________________________________________________________ >>>>> Emanuele Barca Researcher >>>>> Water Research Institute (IRSA-CNR) >>>>> Via De Blasio, 5 70123 Bari (Italy) >>>>> Phone +39 080 5820535 Fax +39 080 5313365 >>>>> Mobile +39 340 3420689 >>>>> _________________________________________________________ >>>>> >>>>> >>>>> >>>>> --- >>>>> Questa e-mail è stata controllata per individuare virus con Avast >>>>> antivirus. >>>>> https://www.avast.com/antivirus >>>>> >>>>> [[alternative HTML version deleted]] >>>>> >>>>> _______________________________________________ >>>>> R-sig-Geo mailing list >>>>> R-sig-Geo@r-project.org >>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>> >>>> >>>> -- >>>> Edzer Pebesma >>>> Institute for Geoinformatics >>>> Heisenbergstrasse 2, 48151 Muenster, Germany >>>> Phone: +49 251 8333081 >>>> >>>> -------------- next part -------------- >>>> A non-text attachment was scrubbed... >>>> Name: pEpkey.asc >>>> Type: application/pgp-keys >>>> Size: 2472 bytes >>>> Desc: not available >>>> URL: >>>> <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190810/b5cb9d77/attachment-0001.bin> >>>> >>>> >>>> >>>> >>>> ------------------------------ >>>> >>>> Subject: Digest Footer >>>> >>>> _______________________________________________ >>>> R-sig-Geo mailing list >>>> R-sig-Geo@r-project.org >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>> >>>> >>>> ------------------------------ >>>> >>>> End of R-sig-Geo Digest, Vol 192, Issue 7 >>>> ***************************************** >>> >>> _______________________________________________ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> -- >> Edzer Pebesma >> Institute for Geoinformatics >> Heisenbergstrasse 2, 48151 Muenster, Germany >> Phone: +49 251 8333081 >> >> -------------- next part -------------- >> A non-text attachment was scrubbed... >> Name: pEpkey.asc >> Type: application/pgp-keys >> Size: 2472 bytes >> Desc: not available >> URL: >> <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190815/7a4eef09/attachment-0001.bin> >> >> >> >> ------------------------------ >> >> Subject: Digest Footer >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> >> >> ------------------------------ >> >> End of R-sig-Geo Digest, Vol 192, Issue 12 >> ****************************************** > > _______________________________________________ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081
pEpkey.asc
Description: application/pgp-keys
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo