[R-sig-Geo] regression-kriging and co-kriging (Edzer Pebesma)
Edzer Pebesma
edzer@pebe@m@ @end|ng |rom un|-muen@ter@de
Thu Aug 15 12:50:12 CEST 2019
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-request using r-project.org ha scritto:
>> Send R-sig-Geo mailing list submissions to
>> r-sig-geo using 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-request using r-project.org
>>
>> You can reach the person managing the list at
>> r-sig-geo-owner using 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 <egampadu using gmail.com>
>> To: r-sig-geo using r-project.org
>> Subject: [R-sig-Geo] Error running codes
>> Message-ID:
>> <CAC4+n+yRavRHc_MCoAtV-ao_-SH7QyO8ZukZaoJqBNaBgJXmPQ using 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: egampadu using gmail.com <egampadu using 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.pebesma using uni-muenster.de>
>> To: r-sig-geo using r-project.org
>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
>> Message-ID: <bdae5226-7023-7e30-1927-5fb66ff413d1 using 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.pebesma using uni-muenster.de>
>>>> To: r-sig-geo using r-project.org
>>>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
>>>> Message-ID: <da76da51-e46b-8a8d-4952-7c3b85c19687 using 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 using 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 using data$Piezo.pred
>>>>> aaaa <- na.omit(g.cok.pred)
>>>>> g.cok.coords <- g.cok using 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 using 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 using 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 using 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 using 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 using 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/78eaf340/attachment.bin>
More information about the R-sig-Geo
mailing list