Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 1 hour 36 min ago

Re: regression-kriging and co-kriging (Edzer Pebesma)

Thu, 08/15/2019 - 05:50


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 [hidden email] ha scritto:
>> Send R-sig-Geo mailing list submissions to
>>     [hidden email]
>>
>> 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
>>     [hidden email]
>>
>> You can reach the person managing the list at
>>     [hidden email]
>>
>> 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 <[hidden email]>
>> To: [hidden email]
>> Subject: [R-sig-Geo] Error running codes
>> Message-ID:
>>     <[hidden email]>
>> 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: [hidden email] <[hidden email]>*
>>
>>
>> *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 <[hidden email]>
>> To: [hidden email]
>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
>> Message-ID: <[hidden email]>
>> 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 <[hidden email]>
>>>> To: [hidden email]
>>>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
>>>> Message-ID: <[hidden email]>
>>>> 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
>>>>> [hidden email]
>>>>> 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
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>>
>>>> ------------------------------
>>>>
>>>> End of R-sig-Geo Digest, Vol 192, Issue 7
>>>> *****************************************
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> 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
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>> ------------------------------
>>
>> End of R-sig-Geo Digest, Vol 192, Issue 12
>> ******************************************
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (2K) Download Attachment

Re: regression-kriging and co-kriging (Edzer Pebesma)

Thu, 08/15/2019 - 05:31
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 [hidden email] ha scritto:
> Send R-sig-Geo mailing list submissions to
> [hidden email]
>
> 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
> [hidden email]
>
> You can reach the person managing the list at
> [hidden email]
>
> 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 <[hidden email]>
> To: [hidden email]
> Subject: [R-sig-Geo] Error running codes
> Message-ID:
> <[hidden email]>
> 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: [hidden email] <[hidden email]>*
>
>
> *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 <[hidden email]>
> To: [hidden email]
> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
> Message-ID: <[hidden email]>
> 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 <[hidden email]>
>>> To: [hidden email]
>>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
>>> Message-ID: <[hidden email]>
>>> 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
>>>> [hidden email]
>>>> 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
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>> ------------------------------
>>>
>>> End of R-sig-Geo Digest, Vol 192, Issue 7
>>> *****************************************
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> ------------------------------
>
> End of R-sig-Geo Digest, Vol 192, Issue 12
> ******************************************
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: regression-kriging and co-kriging

Thu, 08/15/2019 - 01:33


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 <[hidden email]>
>> To: [hidden email]
>> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
>> Message-ID: <[hidden email]>
>> 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
>>> [hidden email]
>>> 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
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>> ------------------------------
>>
>> End of R-sig-Geo Digest, Vol 192, Issue 7
>> *****************************************
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (2K) Download Attachment

Error running codes

Wed, 08/14/2019 - 23:51
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: [hidden email] <[hidden email]>*


*skype: enoch.ampadu*
*The highest evidence of nobility is self-control*.

*A simple act of kindness creates an endless ripple*.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Random forest for forest cover classification

Tue, 08/13/2019 - 19:56
Hi Enoch,
Here is a book chapter that provides R scripts for using random forest to
classify categorical data: 'Predicting seabed hardness using random forest
in R' 2013.

   - DOI:
   - 10.1016/B978-0-12-411511-8.00011-6.
   - In book: Data Mining Applications with R.
   -
   - Publisher: Elsevier.
   - Editors: Y. Zhao and Y. Cen

If you are looking for further materials, the R package, spm, provides more
advanced functions, including RFcv with an example for categorical data.

Hope this helps.
Kind regards,
Jin


On Tue, 13 Aug. 2019, 18:50 Enoch Gyamfi Ampadu, <[hidden email]> wrote:

> Dear List,
>
> Please I am now learning R to apply to my analysis. I am doing a forest
> cover classification using Random forest. So far I have stacked the image,
> done atmospheric correction, and clipped out my study area. I have created
> some training sites of the cover classes in ArcMap and imported into R as
> well as the clipped image of my study sites. I have been doing try and
> error with some online codes. I will be glad if you could provide me with
> some codes to apply.
>
> 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: [hidden email] <[hidden email]>*
>
>
> *skype: enoch.ampadu*
> *The highest evidence of nobility is self-control*.
>
> *A simple act of kindness creates an endless ripple*.
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

rgdal::writeOGR with driver='ESRI Shapefile' converts Polygon object into a hole

Tue, 08/13/2019 - 12:50
Dear list,

I have a single Polygons object containing multiple Polygon objects
that share a common border. When I output this using writeOGR() one of
the Polygon objects becomes a hole, as the following example shows.

Create a Polygons object containing two adjoining Polygon objects
> library(rgdal)
> r1 <- rbind(c(1,1),c(1,2),c(2,2),c(2,1),c(1,1))
> r2 <- r1; r2[,1] <- r2[,1]+1
> Ps1 = Polygons(list(Polygon(r1),Polygon(r2)),ID=1)
> SPDF = SpatialPolygonsDataFrame( SpatialPolygons(list(Ps1)), data.frame(Example=c("Minimal")))

Perform a write/readOGR() cycle
> fn <- tempfile()
> writeOGR(SPDF, fn, layer='test', driver='ESRI Shapefile')
> SPDF2 <- readOGR(dsn=fn,layer='test')

Second Polygon object is now a hole
> sapply(SPDF2@polygons[[1]]@Polygons,slot,"hole")
[1] FALSE  TRUE

I see from the sp documentation that "Polygon objects belonging to a
Polygons object should either not overlap one-other, or should be
fully included" but I am not sure how this relates to bordering
Polygon objects. I would welcome any advice as to whether what I am
asking of writeOGR is reasonable?

Thanks in advance for your time,
Phil

> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United
Kingdom.1252    LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C
                         LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgdal_1.4-4 sp_1.3-1

loaded via a namespace (and not attached):
[1] compiler_3.5.1  tools_3.5.1     yaml_2.2.0      grid_3.5.1
lattice_0.20-35

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Random forest for forest cover classification

Tue, 08/13/2019 - 10:30
Hello Enoch
I think segmentation is a good idea, it will group pixels by similar
characteristics, creating polygons (don't forget to validate the
generated polygons!).
then you can do classification on those polygons using polygons properties
(including shape indices, but also the values obtained from the colours of
the image) using random forest.

Carlos Alberto

On Tue, Aug 13, 2019 at 8:41 AM Juan Pablo Carranza <[hidden email]>
wrote:

> Dear Enoch, I strongly recomend to try an segmentation approach instead a
> classification. However, here's some code you can use:
>
> library(rgdal)
> library(raster)
> library(caret)
> library(kerndwd)
> library(sp)
> library(RColorBrewer)
> library(dplyr)
> library(leaflet)
> library(sf)
> img <- brick("******.tiff")
> names(img) <- paste0("B", c(1:*))
> plot(img)
> library(sf)
> sample <- st_read("******gpkg")
> sample  = as(sample, 'Spatial')
> sample$class = as.factor(sample$class) ; responseCol <- "class"
> pixels = data.frame(matrix(vector(), nrow = 0, ncol = length(names(img)) +
> 1))
>
> # A function to overlap de sample polygons on the image.
> 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)
>   }
> }
>
> # Define a number of pixels to use in the training of the model. More
> pixels, more accuracy but more computational time.
> n <- 1000
>
> # Take a subsample in order to reduce computational times.
> spixeles <- pixeles[sample(1:nrow(pixeles), n ), ]
>
> # Train model.
> modFit_rf <- train(as.factor(class) ~ ., method = "nnet", data = spixeles)
> modFit_rf
>
> # Predict.
> beginCluster()
> preds_rf <- clusterR(img, raster::predict, args = list(model = modFit_rf))
> endCluster()
>
> # Mapping results
> maxpixels =  1e+06
> m <- leaflet() %>%
>   addProviderTiles('Esri.WorldImagery') %>%
>   addRasterImage(preds_rf , opacity = 0.6) %>%
>   addScaleBar(position = "topright") %>%
>   addMiniMap( position = "bottomleft")
> m
>
> # Export raster
> writeRaster(preds_rf, "estimacion.tif")
>
>
>
> *Juan Pablo Carranza*
> Mgter. en Administración Pública
> Lic. en Economía
>
>
>
>
> El mar., 13 ago. 2019 a las 5:50, Enoch Gyamfi Ampadu (<[hidden email]
> >)
> escribió:
>
> > Dear List,
> >
> > Please I am now learning R to apply to my analysis. I am doing a forest
> > cover classification using Random forest. So far I have stacked the
> image,
> > done atmospheric correction, and clipped out my study area. I have
> created
> > some training sites of the cover classes in ArcMap and imported into R as
> > well as the clipped image of my study sites. I have been doing try and
> > error with some online codes. I will be glad if you could provide me with
> > some codes to apply.
> >
> > 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: [hidden email] <[hidden email]>*
> >
> >
> > *skype: enoch.ampadu*
> > *The highest evidence of nobility is self-control*.
> >
> > *A simple act of kindness creates an endless ripple*.
> >
> >         [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Random forest for forest cover classification

Tue, 08/13/2019 - 07:41
Dear Enoch, I strongly recomend to try an segmentation approach instead a
classification. However, here's some code you can use:

library(rgdal)
library(raster)
library(caret)
library(kerndwd)
library(sp)
library(RColorBrewer)
library(dplyr)
library(leaflet)
library(sf)
img <- brick("******.tiff")
names(img) <- paste0("B", c(1:*))
plot(img)
library(sf)
sample <- st_read("******gpkg")
sample  = as(sample, 'Spatial')
sample$class = as.factor(sample$class) ; responseCol <- "class"
pixels = data.frame(matrix(vector(), nrow = 0, ncol = length(names(img)) +
1))

# A function to overlap de sample polygons on the image.
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)
  }
}

# Define a number of pixels to use in the training of the model. More
pixels, more accuracy but more computational time.
n <- 1000

# Take a subsample in order to reduce computational times.
spixeles <- pixeles[sample(1:nrow(pixeles), n ), ]

# Train model.
modFit_rf <- train(as.factor(class) ~ ., method = "nnet", data = spixeles)
modFit_rf

# Predict.
beginCluster()
preds_rf <- clusterR(img, raster::predict, args = list(model = modFit_rf))
endCluster()

# Mapping results
maxpixels =  1e+06
m <- leaflet() %>%
  addProviderTiles('Esri.WorldImagery') %>%
  addRasterImage(preds_rf , opacity = 0.6) %>%
  addScaleBar(position = "topright") %>%
  addMiniMap( position = "bottomleft")
m

# Export raster
writeRaster(preds_rf, "estimacion.tif")



*Juan Pablo Carranza*
Mgter. en Administración Pública
Lic. en Economía




El mar., 13 ago. 2019 a las 5:50, Enoch Gyamfi Ampadu (<[hidden email]>)
escribió:

> Dear List,
>
> Please I am now learning R to apply to my analysis. I am doing a forest
> cover classification using Random forest. So far I have stacked the image,
> done atmospheric correction, and clipped out my study area. I have created
> some training sites of the cover classes in ArcMap and imported into R as
> well as the clipped image of my study sites. I have been doing try and
> error with some online codes. I will be glad if you could provide me with
> some codes to apply.
>
> 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: [hidden email] <[hidden email]>*
>
>
> *skype: enoch.ampadu*
> *The highest evidence of nobility is self-control*.
>
> *A simple act of kindness creates an endless ripple*.
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Random forest for forest cover classification

Tue, 08/13/2019 - 03:49
Dear List,

Please I am now learning R to apply to my analysis. I am doing a forest
cover classification using Random forest. So far I have stacked the image,
done atmospheric correction, and clipped out my study area. I have created
some training sites of the cover classes in ArcMap and imported into R as
well as the clipped image of my study sites. I have been doing try and
error with some online codes. I will be glad if you could provide me with
some codes to apply.

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: [hidden email] <[hidden email]>*


*skype: enoch.ampadu*
*The highest evidence of nobility is self-control*.

*A simple act of kindness creates an endless ripple*.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

regression-kriging and co-kriging

Mon, 08/12/2019 - 13:21
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?

regards

emanuele barca
------------------------------
>
> Message: 2
> Date: Sat, 10 Aug 2019 10:41:38 +0200
> From: Edzer Pebesma <[hidden email]>
> To: [hidden email]
> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging
> Message-ID: <[hidden email]>
> 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
>> [hidden email]
>> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> ------------------------------
>
> End of R-sig-Geo Digest, Vol 192, Issue 7
> *****************************************
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

geobr: easy access to official spatial data sets of Brazil

Mon, 08/12/2019 - 11:29
Dear all,

I'm glad to annouce geobr <https://github.com/ipeaGIT/geobr>, an R package
that provides easy and quick access to official spatia data sets of Brazil.
The package currently includes several data sets for various years such as
states, regions, municipalities, census tracts, statistical grid and
others.  geobr is available on CRAN, and the github repository is here:
https://github.com/ipeaGIT/geobr. I hope this can be useful to some of you
in the group.

The package was developed by my team and I at the Institute for Applied
Economic Research (Ipea). We are constantly working to expand and improve
the package, so if you have any suggestions/contributions, please feel free
to open an issue on Github.

best,

Rafael H. M. Pereira
Institute for Applied Economic Research (Ipea) - Brazil

Department of Urban, Regional and Environmental studies and policies (DIRUR)



Twitter - @UrbanDemog  <https://twitter.com/UrbanDemog>

Blog - urbandemographics.blogspot.com
<https://urbandemographics.blogspot.com.br/>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [FORGED] spatstat - user-supplied list of point patterns in envelope() error

Mon, 08/12/2019 - 07:01
Hi Rolf,

Thank you for the prompt reply.

Now that you mentioned it, it seems such a silly mistake for me to make!

I appreciate the help.

Kind regards,
Joe

On Mon, Aug 12, 2019 at 3:31 AM Rolf Turner <[hidden email]> wrote:
>
>
> Dear Joe,
>
> Further to my previous response to your post:  I have now been in touch
> with Adrian Baddeley and he has confirmed that the problem is as I
> conjectured:   The patterns in "pp_list" need to be compatible with
> "xx".  The error message was a bit misleading.  (Adrian has now adjusted
> the error message so as to make it more clear what the problem actually is.)
>
> However there was also a genuine bug (that was revealed by your
> question; thank you!!!) in the software, so even if pp_list had been a
> list of patterns of the appropriate type, you would still have got an
> error.  (A spurious error, caused by the bug; different message of course.)
>
> The bug has been been fixed in the development version of spatstat
> (1.60-1.022).  You will need to install that from github:
>
> library(remotes)
> install_github('spatstat/spatstat')
>
> or wait for the next "official release" of spatstat on CRAN.  Such
> releases happen about every three months; the last one was on
> 23/June/2019.
>
> cheers,
>
> Rolf
>
> > On 12/08/19 8:51 AM, Joe Lewis wrote:
> >
> >> Dear list,
> >>
> >> I'm trying to use a user-supplied list of point patterns in envelope()
> >> rather than test against the CSR. Page 400 of Spatial Point Patterns:
> >> Methodology and Applications in R states that "the argument simulate
> >> can be a list of point patterns". However, I get the following error
> >> when I try  to supply a list of ppp:
> >>
> >>   ekls <- envelope.lpp(Y =  xx, fun = linearK, nsim=5, simulate =
> >> pp_list)
> >> Error in envelopeEngine(X = X, fun = fun, simul = simrecipe, nsim =
> >> nsim,  :
> >>    ‘simulate’ should be an expression, or a list of point patterns
> >>
> >> for reference,
> >>
> >> xx
> >> Point pattern on linear network
> >> 680 points
> >> Linear network with 22 vertices and 21 lines
> >> Enclosing window: rectangle = [284086.69, 309740] x [709900, 726547.7]
> >> metres
> >>
> >> class(xx)
> >> [1] "lpp" "ppx"
> >>
> >> pp_list
> >> List of point patterns
> >>
> >> Component 1:
> >> Planar point pattern: 513 points
> >> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
> >>
> >> Component 2:
> >> Planar point pattern: 422 points
> >> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
> >>
> >> Component 3:
> >> Planar point pattern: 495 points
> >> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
> >>
> >> Component 4:
> >> Planar point pattern: 557 points
> >> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
> >>
> >> Component 5:
> >> Planar point pattern: 576 points
> >> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
> >>
> >> class(pp_list)
> >> [1] "ppplist" "solist"  "anylist" "listof"  "list"
> >>
> >> Any ideas why the error is occurring? Thanks.
> >
> > It's hard to say without having a *reproducible* example.  Since we
> > don't have access to "xx" or to "pp.list" we cannot experiment to see
> > what's going on.
> >
> > One problem that leaps out at me --- although it doesn't seem that this
> > should trigger the error message that you received --- is that the
> > entries of pp_list appear to be *planar point patterns* (of class "ppp")
> > whereas "xx" is a pattern on a linear network (of class "lpp").
> > Consequently there is a fundamental incompatibility here.
> >
> > However I don't see why you would get the error message that you did.  I
> > am CC-ing this email to Adrian Baddeley who has more insight than I, and
> > may be able to point you in the right direction.  Adrian is kind of
> > overwhelmed with work at the moment, so it may be a while till you hear
> > from him.
> >
> > If you provide a reproducible example I *may* be able to help.
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [FORGED] spatstat - user-supplied list of point patterns in envelope() error

Sun, 08/11/2019 - 21:31

Dear Joe,

Further to my previous response to your post:  I have now been in touch
with Adrian Baddeley and he has confirmed that the problem is as I
conjectured:   The patterns in "pp_list" need to be compatible with
"xx".  The error message was a bit misleading.  (Adrian has now adjusted
the error message so as to make it more clear what the problem actually is.)

However there was also a genuine bug (that was revealed by your
question; thank you!!!) in the software, so even if pp_list had been a
list of patterns of the appropriate type, you would still have got an
error.  (A spurious error, caused by the bug; different message of course.)

The bug has been been fixed in the development version of spatstat
(1.60-1.022).  You will need to install that from github:

library(remotes)
install_github('spatstat/spatstat')

or wait for the next "official release" of spatstat on CRAN.  Such
releases happen about every three months; the last one was on
23/June/2019.

cheers,

Rolf

> On 12/08/19 8:51 AM, Joe Lewis wrote:
>
>> Dear list,
>>
>> I'm trying to use a user-supplied list of point patterns in envelope()
>> rather than test against the CSR. Page 400 of Spatial Point Patterns:
>> Methodology and Applications in R states that "the argument simulate
>> can be a list of point patterns". However, I get the following error
>> when I try  to supply a list of ppp:
>>
>>   ekls <- envelope.lpp(Y =  xx, fun = linearK, nsim=5, simulate =
>> pp_list)
>> Error in envelopeEngine(X = X, fun = fun, simul = simrecipe, nsim =
>> nsim,  :
>>    ‘simulate’ should be an expression, or a list of point patterns
>>
>> for reference,
>>
>> xx
>> Point pattern on linear network
>> 680 points
>> Linear network with 22 vertices and 21 lines
>> Enclosing window: rectangle = [284086.69, 309740] x [709900, 726547.7]
>> metres
>>
>> class(xx)
>> [1] "lpp" "ppx"
>>
>> pp_list
>> List of point patterns
>>
>> Component 1:
>> Planar point pattern: 513 points
>> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>>
>> Component 2:
>> Planar point pattern: 422 points
>> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>>
>> Component 3:
>> Planar point pattern: 495 points
>> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>>
>> Component 4:
>> Planar point pattern: 557 points
>> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>>
>> Component 5:
>> Planar point pattern: 576 points
>> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>>
>> class(pp_list)
>> [1] "ppplist" "solist"  "anylist" "listof"  "list"
>>
>> Any ideas why the error is occurring? Thanks.
>
> It's hard to say without having a *reproducible* example.  Since we
> don't have access to "xx" or to "pp.list" we cannot experiment to see
> what's going on.
>
> One problem that leaps out at me --- although it doesn't seem that this
> should trigger the error message that you received --- is that the
> entries of pp_list appear to be *planar point patterns* (of class "ppp")
> whereas "xx" is a pattern on a linear network (of class "lpp").
> Consequently there is a fundamental incompatibility here.
>
> However I don't see why you would get the error message that you did.  I
> am CC-ing this email to Adrian Baddeley who has more insight than I, and
> may be able to point you in the right direction.  Adrian is kind of
> overwhelmed with work at the moment, so it may be a while till you hear
> from him.
>
> If you provide a reproducible example I *may* be able to help.
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [FORGED] spatstat - user-supplied list of point patterns in envelope() error

Sun, 08/11/2019 - 17:59

On 12/08/19 8:51 AM, Joe Lewis wrote:

> Dear list,
>
> I'm trying to use a user-supplied list of point patterns in envelope()
> rather than test against the CSR. Page 400 of Spatial Point Patterns:
> Methodology and Applications in R states that "the argument simulate
> can be a list of point patterns". However, I get the following error
> when I try  to supply a list of ppp:
>
>   ekls <- envelope.lpp(Y =  xx, fun = linearK, nsim=5, simulate = pp_list)
> Error in envelopeEngine(X = X, fun = fun, simul = simrecipe, nsim = nsim,  :
>    ‘simulate’ should be an expression, or a list of point patterns
>
> for reference,
>
> xx
> Point pattern on linear network
> 680 points
> Linear network with 22 vertices and 21 lines
> Enclosing window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>
> class(xx)
> [1] "lpp" "ppx"
>
> pp_list
> List of point patterns
>
> Component 1:
> Planar point pattern: 513 points
> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>
> Component 2:
> Planar point pattern: 422 points
> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>
> Component 3:
> Planar point pattern: 495 points
> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>
> Component 4:
> Planar point pattern: 557 points
> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>
> Component 5:
> Planar point pattern: 576 points
> window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres
>
> class(pp_list)
> [1] "ppplist" "solist"  "anylist" "listof"  "list"
>
> Any ideas why the error is occurring? Thanks.
It's hard to say without having a *reproducible* example.  Since we
don't have access to "xx" or to "pp.list" we cannot experiment to see
what's going on.

One problem that leaps out at me --- although it doesn't seem that this
should trigger the error message that you received --- is that the
entries of pp_list appear to be *planar point patterns* (of class "ppp")
whereas "xx" is a pattern on a linear network (of class "lpp").
Consequently there is a fundamental incompatibility here.

However I don't see why you would get the error message that you did.  I
am CC-ing this email to Adrian Baddeley who has more insight than I, and
may be able to point you in the right direction.  Adrian is kind of
overwhelmed with work at the moment, so it may be a while till you hear
from him.

If you provide a reproducible example I *may* be able to help.

cheers,

Rolf Turner

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

spatstat - user-supplied list of point patterns in envelope() error

Sun, 08/11/2019 - 15:51
Dear list,

I'm trying to use a user-supplied list of point patterns in envelope()
rather than test against the CSR. Page 400 of Spatial Point Patterns:
Methodology and Applications in R states that "the argument simulate
can be a list of point patterns". However, I get the following error
when I try  to supply a list of ppp:

 ekls <- envelope.lpp(Y =  xx, fun = linearK, nsim=5, simulate = pp_list)
Error in envelopeEngine(X = X, fun = fun, simul = simrecipe, nsim = nsim,  :
  ‘simulate’ should be an expression, or a list of point patterns

for reference,

xx
Point pattern on linear network
680 points
Linear network with 22 vertices and 21 lines
Enclosing window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres

class(xx)
[1] "lpp" "ppx"

pp_list
List of point patterns

Component 1:
Planar point pattern: 513 points
window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres

Component 2:
Planar point pattern: 422 points
window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres

Component 3:
Planar point pattern: 495 points
window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres

Component 4:
Planar point pattern: 557 points
window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres

Component 5:
Planar point pattern: 576 points
window: rectangle = [284086.69, 309740] x [709900, 726547.7] metres

class(pp_list)
[1] "ppplist" "solist"  "anylist" "listof"  "list"

Any ideas why the error is occurring? Thanks.

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Error when running a Spatial Durbin Error Model without intercept

Sat, 08/10/2019 - 09:34
And for completeness (https://github.com/r-spatial/spatialreg/issues/7):

> dem <- lmSLX(eq, data = COL.OLD, listw = W)
Error in lmSLX(eq, data = COL.OLD, listw = W) :
  object 'dirImps' not found
In addition: Warning message:
In lmSLX(eq, data = COL.OLD, listw = W) :
  model configuration issue: no total impacts
> traceback()
1: lmSLX(eq, data = COL.OLD, listw = W)

Please follow up on github if you are present there.

Roger



--
Roger Bivand
Norwegian School of Economics
Helleveien 30, 5045 Bergen, Norway
[hidden email]


________________________________________
Fra: R-sig-Geo <[hidden email]> på vegne av Roger Bivand <[hidden email]>
Sendt: lørdag 10. august 2019 14.05
Til: Rolando Valdez
Kopi: [hidden email]
Emne: Re: [R-sig-Geo]  Error when running a Spatial Durbin Error Model without intercept

On Fri, 9 Aug 2019, Rolando Valdez wrote:

> Dear list,
>
> I am trying to run a Spatial Durbin Error Model without intercept, however,
> I do get an error:
>
> In spatialreg::errorsarlm(formula = formula, data = data, listw = listw,  :
>  model configuration issue: no total impacts

I can reproduce the problem, with error messages:

> sdem <- errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed",
+                    method = "eigen")
Error in estimable.default(lm.target, cm) :
   `cm' argument must be of type vector, list, or matrix.
In addition: Warning messages:
1: Function errorsarlm moved to the spatialreg package
2: In spatialreg::errorsarlm(formula = formula, data = data, listw =
listw,  :
   model configuration issue: no total impacts
> traceback()
6: stop("`cm' argument must be of type vector, list, or matrix.")
5: estimable.default(lm.target, cm)
4: estimable(lm.target, cm)
3: as.matrix(estimable(lm.target, cm)[, 1:2, drop = FALSE])
2: spatialreg::errorsarlm(formula = formula, data = data, listw = listw,
        na.action = na.action, Durbin = Durbin, etype = etype, method =
method,
        quiet = quiet, zero.policy = zero.policy, interval = interval,
        tol.solve = tol.solve, trs = trs, control = control)
1: errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed", method =
"eigen")

I think that the construction of the specification matrix for
gmodels::estimable() is faulty; will explore whether no-intercept is
supportable.

Roger

>
> This is what I am doing:
>
> library(spatialreg)
> library(spdep)
> data(oldcol, package = "spdep")
> W <- spdep::nb2listw(COL.nb, style = "W")
> eq <- CRIME ~ 0 + INC + HOVAL
> sdem <- errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed",
>                   method = "eigen")
>
> Thanks for any advice.
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: Error when running a Spatial Durbin Error Model without intercept

Sat, 08/10/2019 - 07:05
On Fri, 9 Aug 2019, Rolando Valdez wrote:

> Dear list,
>
> I am trying to run a Spatial Durbin Error Model without intercept, however,
> I do get an error:
>
> In spatialreg::errorsarlm(formula = formula, data = data, listw = listw,  :
>  model configuration issue: no total impacts

I can reproduce the problem, with error messages:

> sdem <- errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed",
+                    method = "eigen")
Error in estimable.default(lm.target, cm) :
   `cm' argument must be of type vector, list, or matrix.
In addition: Warning messages:
1: Function errorsarlm moved to the spatialreg package
2: In spatialreg::errorsarlm(formula = formula, data = data, listw =
listw,  :
   model configuration issue: no total impacts
> traceback()
6: stop("`cm' argument must be of type vector, list, or matrix.")
5: estimable.default(lm.target, cm)
4: estimable(lm.target, cm)
3: as.matrix(estimable(lm.target, cm)[, 1:2, drop = FALSE])
2: spatialreg::errorsarlm(formula = formula, data = data, listw = listw,
        na.action = na.action, Durbin = Durbin, etype = etype, method =
method,
        quiet = quiet, zero.policy = zero.policy, interval = interval,
        tol.solve = tol.solve, trs = trs, control = control)
1: errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed", method =
"eigen")

I think that the construction of the specification matrix for
gmodels::estimable() is faulty; will explore whether no-intercept is
supportable.

Roger

>
> This is what I am doing:
>
> library(spatialreg)
> library(spdep)
> data(oldcol, package = "spdep")
> W <- spdep::nb2listw(COL.nb, style = "W")
> eq <- CRIME ~ 0 + INC + HOVAL
> sdem <- errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed",
>                   method = "eigen")
>
> Thanks for any advice.
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: regression-kriging and co-kriging

Sat, 08/10/2019 - 03:41
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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (2K) Download Attachment

Re: st_sample with raster

Sat, 08/10/2019 - 03:38
sf::st_sample has a size argument, which can be specified as a vector to
operate per feature. If you give that a value equal to

st_area(obj) * obj$required_density

and maybe set exact = FALSE, you should get the densities you want.

On 7/30/19 12:25 AM, Tyler Frazier via R-sig-Geo wrote:
> Hi All,
>
> Spatstat::rpoint has an argument f that accepts a pixel image object (class .im) to generate points based on the proportionate probability density of that raster.  Just wondering, does sf::st_sample or perhaps another function have a similar available method?
>
> Thanks so much!
> Ty
>
> Tyler Frazier, Ph.D.
> Lecturer of Interdisciplinary Studies
> Data Science Program
> College of William and Mary
> Webpage: http://tjfrazier.people.wm.edu
> Phone: +01 (757) 386-1269
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (2K) Download Attachment

Error when running a Spatial Durbin Error Model without intercept

Fri, 08/09/2019 - 16:16
Dear list,

I am trying to run a Spatial Durbin Error Model without intercept, however,
I do get an error:

 In spatialreg::errorsarlm(formula = formula, data = data, listw = listw,  :
  model configuration issue: no total impacts

This is what I am doing:

library(spatialreg)
library(spdep)
data(oldcol, package = "spdep")
W <- spdep::nb2listw(COL.nb, style = "W")
eq <- CRIME ~ 0 + INC + HOVAL
sdem <- errorsarlm(eq, data = COL.OLD, listw = W, etype = "emixed",
                   method = "eigen")

Thanks for any advice.
--
Rol~

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Pages