[R-sig-Geo] regression-kriging and co-kriging (Edzer Pebesma)
Emanuele Barca
em@nue|e@b@rc@ @end|ng |rom b@@|r@@@cnr@|t
Thu Aug 15 12:31:04 CEST 2019
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?
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
> ******************************************
More information about the R-sig-Geo
mailing list