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: 55 min 15 sec ago

Re: Getting spgm Models to work at Substate Levels - SAMHSA NSDUH Data - with Added nb Links

Thu, 08/08/2019 - 14:41
On Thu, 8 Aug 2019, Stuart Reece wrote:

> Thankyou so much again Roger.
>
> Both the connectivity maps and regressions work fine with the state data
> so I hope you will not mind if I use the substate data from the link
> mentioned earlier at :
>
> https://www.samhsa.gov/data/sites/default/files/cbhsq-reports/NSDUHsubstateShapeFile2016/ShapeFile2016.zip
>
> and then just subset for one state.  New York would be good to choose as
> it has Richmond Island off the southern tip of Long Island which is not
> connected to any other area.  So it naturally forms a major point of
> interest.
>
> I have done my level best to follow your instructions very carefully.
Thanks for your patience - now there is something to work on. It is easier
to extract the code if you post plain text only, so my abbreviation is:

library(sf)
USC <- st_read(dsn="SubstateRegionData141516.shp")
USCReduced <- USC[,c(1:6,10,14,18,22,26,30,34,38,42,46,50,66)]
NYsf <- subset(USCReduced, ST_NAME=="New York")
NYsf$Ix <- 0:14
NYsf$YearDt <- as.POSIXct(as.Date("01/01/2015", format="%d/%m/%Y"),
tz="GMT")
NYsf <- NYsf[,c(19,20,1:5,9,6,7,10,14,15,12,8,13,16,17,18)]
library(spdep)
NYsfnb <- poly2nb(NYsf)
NYsfnb
NYsfnb[[8]] <- as.integer(5)
NYsfnb[[5]] <- as.integer(unique(sort(c(NYsfnb[[5]], 8))))
NYsflw  <- nb2listw(NYsfnb)
library(splm)
ml_obj <- spml(asinh(smiyr) ~ cigmon * asinh(mrjmon), data=NYsf,
  listw=NYsflw, lag=TRUE, effect="individual", model="random",
  spatial.error="b"))
gm_obj <- spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data=NYsf,
  listw=NYsflw, lag=TRUE, moments="fullweights", method="g2sls",
  model="random", spatial.error=TRUE)
# Error in x[, ii] : subscript out of bounds
gm_obj <- spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data=NYsf,
  listw=NYsflw, lag=TRUE, moments="weights", method="g2sls",
  model="random", spatial.error=TRUE)
# Error in x[, ii] : subscript out of bounds
gm_obj <- spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data=NYsf,
  listw=NYsflw, lag=TRUE, moments="initial", method="g2sls",
  model="random", spatial.error=TRUE)
# Error in x[, ii] : subscript out of bounds
gm_obj <- spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data=NYsf,
  listw=NYsflw, moments="fullweights", spatial.error=TRUE)
# Error in x[, ii] : subscript out of bounds

Could you please try to run the model through plm - this subset seems only
to have one time period? Perhaps plm is more forgiving than splm? Maybe
that is confusing the internals?

summary(lm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data=NYsf))
library(spatialreg)
summary(lagsarlm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data=NYsf,
  listw=NYsflw))

both complete without trouble?

Roger

>
> My code follows - obviously including the path to the file on my machine.
>
> The connectivity map - closely following your own model script of course -  is perfect.  The inserted link works just exactly correctly as you will see.
>
> The spml model also works OK albeit with some mysterious error which I do not well understand, but I do not think it matters much either.
>
>
>
> However the spgm model does not run at all.
>
> Of course this is also the most important model.
>
> This throws the usual error "Error in x[, ii] : subscript out of bounds".
>
> I have included the traceback notes, but I do not understand them at all.
>
>
>
> Thankyou again for all of your kind and gracious advice.
>
>
>
> Yours sincerely,
>
> Stuart.
>
>
>
>
>
> ################################################################################################################################################################################################
>
> ################## Using NYC as a Reproducible Example
>
> # Substate Shapefile - Single Year
>
>
>
> USC <- st_read(dsn=" [path to file]  /NSDUH/SubstateRegionData141516.shp")
>
> dim(USC)
>
> names(USC)
>
>
>
> USCReduced <- USC[,c(1:6,10,14,18,22,26,30,34,38,42,46,50,66)]
>
> dim(USCReduced)
>
> names(USCReduced)
>
>
>
> NYsf <- subset(USCReduced, ST_NAME=="New York")
>
> NYsf$Ix <- 0:14
>
> NYsf$YearDt <- as.POSIXct(as.Date("01/01/2015", format="%d/%m/%Y"), tz="GMT")
>
> class(NYsf$YearDt)
>
> dim(NYsf)
>
> names(NYsf)
>
> NYsf <- NYsf[,c(19,20,1:5,9,6,7,10,14,15,12,8,13,16,17,18)]
>
> dim(NYsf)
>
> names(NYsf)
>
> glimpse(NYsf)
>
> print(NYsf, n=15)
>
>
>
> NYsfnb <- poly2nb(NYsf)
>
> str(NYsfnb, max.level=0)
>
> str(NYsfnb, max.level=1)
>
> summary(NYsfnb)
>
>
>
> attr(NYsfnb, "region.id")[1:8]
>
>
>
> NYsfnbb <- NYsfnb
>
>
>
> attr(NYsfnbb, "region.id")
>
> NYsfnbb[[8]]                           # == Region 322 on the Original map
>
> NYsfnbb[[8]] <- as.integer(5)
>
> NYsfnbb[[8]]                           # == Region 322 on the Original map
>
> NYsfnbb[[5]]                           # == Region 319 on the Original map
>
> NYsfnbb[[5]] <- as.integer(c(6,7,8))
>
> NYsfnbb[[5]]
>
> card(NYsfnbb)
>
> summary(NYsfnbb)
>
>
>
> NYsflw  <- nb2listw(NYsfnb, zero.policy=TRUE)
>
> NYsflww <- nb2listw(NYsfnbb)
>
> str(NYsfnbb, max.level = 0)
>
> str(NYsflww, max.level = 0)
>
>
>
> can.be.simmed(NYsflww)
>
>
>
> ############### US Substate Test Mapping::
>
> coordsSubstateNY <- st_coordinates(st_centroid(st_geometry(NYsf), of_largest_polygon = TRUE))
>
>
>
> ## SubState_q1
>
> dxxx <- diffnb(NYsfnb, NYsfnbb)
>
> plot(st_geometry(NYsf), border = "grey50")
>
> plot(NYsfnb, coordsSubstateNY, add = TRUE, lwd=0.5)
>
> plot(dxxx, coordsSubstateNY, add = TRUE, col = "red", lwd=2)
>
> title(main=paste("Differences (red) in New York Sub-State GAL Queen Weights (Black) First Order - USA", cex.main=1))
>
>
>
> ### Checking w Regressions
>
> summary(spml(asinh(smiyr) ~ cigmon * asinh(mrjmon),
>
>             data=NYsf, listw=NYsflww,
>
>             lag=TRUE, effect="individual",
>
>             model="random", spatial.error="b"))
>
>
>
> summary(spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon),
>
>             data=NYsf, listw=NYsflww,
>
>             lag=TRUE, moments="fullweights", method="g2sls",
>
>             model="random", spatial.error=TRUE))
>
>
>
> #    Error in x[, ii] : subscript out of bounds
>
> #
>
> # > traceback()
>
> # 13: na.omit.data.frame(list(`yend[, i]` = c(0, 0, 0, 0, 0, 0, 0,
>
> #                                             0, 0, 0, 0, 0, 0, 0, 0), H = numeric(0)))
>
> # 12: na.omit(list(`yend[, i]` = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
>
> #                                  0, 0, 0, 0), H = numeric(0)))
>
> # 11: model.frame.default(formula = yend[, i] ~ H - 1, drop.unused.levels = TRUE)
>
> # 10: stats::model.frame(formula = yend[, i] ~ H - 1, drop.unused.levels = TRUE)
>
> # 9: eval(mf, parent.frame())
>
> # 8: eval(mf, parent.frame())
>
> # 7: lm(yend[, i] ~ H - 1)
>
> # 6: fitted.values(lm(yend[, i] ~ H - 1))
>
> # 5: spgm.tsls(ywithin, wywithin, Xwithin, Hwithin)
>
> # 4: ivplm.w2sls(Y = y, X = x, lag = TRUE, listw = Ws, listw2 = Ws2,
>
> #                twow = twow, lag.instruments = lag.instruments, T = T, N = N,
>
> #                NT = NT)
>
> # 3: spsarargm(formula = formula, data = data, index = index, listw = listw,
>
> #              listw2 = listw2, moments = moments, lag = lag, endog = endog,
>
> #              instruments = instruments, verbose = verbose, effects = effects,
>
> #              control = control, lag.instruments = lag.instruments, optim.method = optim.method,
>
> #              pars = pars, twow = twow)
>
> # 2: spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data = NYsf, listw = NYsflww,
>
> #         lag = TRUE, moments = "fullweights", model = "random", spatial.error = TRUE)
>
> # 1: summary(spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data = NYsf,
>
> #                 listw = NYsflww, lag = TRUE, moments = "fullweights", model = "random",
>
> #                 spatial.error = TRUE))
>
>
>
>
>
> -----Original Message-----
> From: Roger Bivand <[hidden email]>
> Sent: Thursday, 8 August, 2019 9:35 PM
> To: Stuart Reece <[hidden email]>
> Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
> Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
>
>
>
> On Thu, 8 Aug 2019, Stuart Reece wrote:
>
>
>
>> Thankyou for this advice Roger.
>
>> Happy to provide my work for your review but I am not sure how.
>
>> This list server has a limit of 50KB and these files are about 100MB....
>
>> Love to assist but I would need advice on how best to proceed....???
>
>
>
> Choose a built-in data set, such as Produc, change two states to ALASKA and HAWAII, drop the ones you chose from the map and any other states not on your map, and use that data set, putting the script needed at the head of your example. Then create the weights, initially with no links to ALASKA and HAWAII, and run the model. Next, edit the neighbour object in script form, and try again. Then anyone with your map object (you gave a link), your script, and the plm and splm packages can reproduce your problem.
>
>
>
> Do also run traceback() after errors, and post what you see.
>
>
>
> Roger
>
>
>
>> Many thanks again,
>
>> Stuart.
>
>>
>
>
--
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

Trying to reproduce rLandsat8 package codes

Thu, 08/08/2019 - 14:18
Dear folks,

Could anyone give a support in downloading LandSat images ?

I am trying to reproduce codes from rLandSat8 package.
However I could not move to the next step due to a problem with the file
downloaded.

*#First step: Download a file for a reproducible example like suggested*
#####################
DownloadLandsat <- function(url, output.name) {
  # todo: use RCurl instead of the system call

command.args <- paste0("-c cookies.txt -d 'username=",
usgs.username="myusername", "&password=", usgs.password="mypassword","'
https://earthexplorer.usgs.gov/login")

  # invoke the system call to curl.
  # I'd rather have done this with RCurl but couldn't get it working ;-(
  ret <- system2("curl", command.args, stdout=TRUE, stderr=TRUE)

  command.args <- paste0("-b cookies.txt -L ", url," -o ", output.name)
  ret <- system2("curl", command.args, stdout=TRUE, stderr=TRUE)
  #file.remove(cookies.txt)
  return(ret)
}
setwd("my folder")
curl_download(DownloadLandsat("
http://earthexplorer.usgs.gov/download/4923/LC82040322013219LGN00/STANDARD/INVSVC",
"LC82040322013219LGN00"))

*#Second Step: Use the file downloaded to reproduce the code.  *
####################
setwd("my folder")
product  <- "LC82040322013219LGN00"
l <- ReadLandsat8(product)
Error in file(con, "r") : cannot open the connection
In addition: Warning message:
In file(con, "r") :
  cannot open file 'LC82040322013219LGN00/LC82040322013219LGN00_MTL.txt':
No such file or directory

The file is in my folder and it seems to be downloaded properly.
However, this file has no extension in other words, it is
LC82040322013219LGN00 and not  LC82040322013219LGN00_MTL.txt.

What is wrong?
Is there any other easier way to download LandSat images using R? Long
series for example?

Best regards,

Jackson

        [[alternative HTML version deleted]]

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

Re: Getting spgm Models to work at Substate Levels - SAMHSA NSDUH Data - with Added nb Links

Thu, 08/08/2019 - 08:41
Thankyou so much again Roger.

 

Both the connectivity maps and regressions work fine with the state data so I hope you will not mind if I use the substate data from the link mentioned earlier at :

https://www.samhsa.gov/data/sites/default/files/cbhsq-reports/NSDUHsubstateShapeFile2016/ShapeFile2016.zip 

 

and then just subset for one state.  New York would be good to choose as it has Richmond Island off the southern tip of Long Island which is not connected to any other area.  So it naturally forms a major point of interest.

 

I have done my level best to follow your instructions very carefully.

 

My code follows - obviously including the path to the file on my machine.

The connectivity map - closely following your own model script of course -  is perfect.  The inserted link works just exactly correctly as you will see.

The spml model also works OK albeit with some mysterious error which I do not well understand, but I do not think it matters much either.

 

However the spgm model does not run at all.

Of course this is also the most important model.

This throws the usual error "Error in x[, ii] : subscript out of bounds".

I have included the traceback notes, but I do not understand them at all.

 

Thankyou again for all of your kind and gracious advice.

 

Yours sincerely,

Stuart.

 

 

################################################################################################################################################################################################

################## Using NYC as a Reproducible Example

# Substate Shapefile - Single Year



USC <- st_read(dsn=" [path to file]  /NSDUH/SubstateRegionData141516.shp")

dim(USC)

names(USC)



USCReduced <- USC[,c(1:6,10,14,18,22,26,30,34,38,42,46,50,66)]

dim(USCReduced)

names(USCReduced)



NYsf <- subset(USCReduced, ST_NAME=="New York")

NYsf$Ix <- 0:14

NYsf$YearDt <- as.POSIXct(as.Date("01/01/2015", format="%d/%m/%Y"), tz="GMT")

class(NYsf$YearDt)

dim(NYsf)

names(NYsf)

NYsf <- NYsf[,c(19,20,1:5,9,6,7,10,14,15,12,8,13,16,17,18)]

dim(NYsf)

names(NYsf)

glimpse(NYsf)

print(NYsf, n=15)



NYsfnb <- poly2nb(NYsf)

str(NYsfnb, max.level=0)

str(NYsfnb, max.level=1)

summary(NYsfnb)



attr(NYsfnb, "region.id")[1:8]



NYsfnbb <- NYsfnb



attr(NYsfnbb, "region.id")

NYsfnbb[[8]]                           # == Region 322 on the Original map

NYsfnbb[[8]] <- as.integer(5)

NYsfnbb[[8]]                           # == Region 322 on the Original map

NYsfnbb[[5]]                           # == Region 319 on the Original map

NYsfnbb[[5]] <- as.integer(c(6,7,8))

NYsfnbb[[5]]

card(NYsfnbb)

summary(NYsfnbb)



NYsflw  <- nb2listw(NYsfnb, zero.policy=TRUE)

NYsflww <- nb2listw(NYsfnbb)

str(NYsfnbb, max.level = 0)

str(NYsflww, max.level = 0)



can.be.simmed(NYsflww)



############### US Substate Test Mapping::

coordsSubstateNY <- st_coordinates(st_centroid(st_geometry(NYsf), of_largest_polygon = TRUE))



## SubState_q1

dxxx <- diffnb(NYsfnb, NYsfnbb)

plot(st_geometry(NYsf), border = "grey50")

plot(NYsfnb, coordsSubstateNY, add = TRUE, lwd=0.5)

plot(dxxx, coordsSubstateNY, add = TRUE, col = "red", lwd=2)

title(main=paste("Differences (red) in New York Sub-State GAL Queen Weights (Black) First Order - USA", cex.main=1))



### Checking w Regressions

summary(spml(asinh(smiyr) ~ cigmon * asinh(mrjmon),

             data=NYsf, listw=NYsflww,

             lag=TRUE, effect="individual",

             model="random", spatial.error="b"))



summary(spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon),

             data=NYsf, listw=NYsflww,

             lag=TRUE, moments="fullweights", method="g2sls",

             model="random", spatial.error=TRUE))



#    Error in x[, ii] : subscript out of bounds

#

# > traceback()

# 13: na.omit.data.frame(list(`yend[, i]` = c(0, 0, 0, 0, 0, 0, 0,

#                                             0, 0, 0, 0, 0, 0, 0, 0), H = numeric(0)))

# 12: na.omit(list(`yend[, i]` = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

#                                  0, 0, 0, 0), H = numeric(0)))

# 11: model.frame.default(formula = yend[, i] ~ H - 1, drop.unused.levels = TRUE)

# 10: stats::model.frame(formula = yend[, i] ~ H - 1, drop.unused.levels = TRUE)

# 9: eval(mf, parent.frame())

# 8: eval(mf, parent.frame())

# 7: lm(yend[, i] ~ H - 1)

# 6: fitted.values(lm(yend[, i] ~ H - 1))

# 5: spgm.tsls(ywithin, wywithin, Xwithin, Hwithin)

# 4: ivplm.w2sls(Y = y, X = x, lag = TRUE, listw = Ws, listw2 = Ws2,

#                twow = twow, lag.instruments = lag.instruments, T = T, N = N,

#                NT = NT)

# 3: spsarargm(formula = formula, data = data, index = index, listw = listw,

#              listw2 = listw2, moments = moments, lag = lag, endog = endog,

#              instruments = instruments, verbose = verbose, effects = effects,

#              control = control, lag.instruments = lag.instruments, optim.method = optim.method,

#              pars = pars, twow = twow)

# 2: spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data = NYsf, listw = NYsflww,

#         lag = TRUE, moments = "fullweights", model = "random", spatial.error = TRUE)

# 1: summary(spgm(asinh(smiyr) ~ cigmon * asinh(mrjmon), data = NYsf,

#                 listw = NYsflww, lag = TRUE, moments = "fullweights", model = "random",

#                 spatial.error = TRUE))





-----Original Message-----
From: Roger Bivand <[hidden email]>
Sent: Thursday, 8 August, 2019 9:35 PM
To: Stuart Reece <[hidden email]>
Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List



On Thu, 8 Aug 2019, Stuart Reece wrote:



> Thankyou for this advice Roger.

> Happy to provide my work for your review but I am not sure how.

> This list server has a limit of 50KB and these files are about 100MB....

> Love to assist but I would need advice on how best to proceed....???



Choose a built-in data set, such as Produc, change two states to ALASKA and HAWAII, drop the ones you chose from the map and any other states not on your map, and use that data set, putting the script needed at the head of your example. Then create the weights, initially with no links to ALASKA and HAWAII, and run the model. Next, edit the neighbour object in script form, and try again. Then anyone with your map object (you gave a link), your script, and the plm and splm packages can reproduce your problem.



Do also run traceback() after errors, and post what you see.



Roger



> Many thanks again,

> Stuart.

>


        [[alternative HTML version deleted]]

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

Re: Adding a Few Neighbour Relationships to a nb List

Thu, 08/08/2019 - 06:52
On Thu, 8 Aug 2019, Stuart Reece wrote:

> Thanks Roger.
> Happy to do as you have requested.
> Which server do you suggest??

You gave the link to the map as:

https://www.samhsa.gov/data/report/2014-2016-nsduh-substate-region-shapefile 
leading to
https://www.samhsa.gov/data/sites/default/files/cbhsq-reports/NSDUHsubstateShapeFile2016/ShapeFile2016.zip

I guess that is the input map. Produc is built into the plm package. So
all we need is a script. Your data are no help; creating a reproducible
example very often guides the person with the problem to a resolution even
without external help. I see that your map has 300+ entities - in that
case find another map that matches the setting (US with Alaska and Hawaii
states, drop to 48 to match Produc).

Roger

> Thanks again,
> Stuart.
>
> -----Original Message-----
> From: Roger Bivand <[hidden email]>
> Sent: Thursday, 8 August, 2019 9:35 PM
> To: Stuart Reece <[hidden email]>
> Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
> Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
>
> On Thu, 8 Aug 2019, Stuart Reece wrote:
>
>> Thankyou for this advice Roger.
>> Happy to provide my work for your review but I am not sure how.
>> This list server has a limit of 50KB and these files are about 100MB....
>> Love to assist but I would need advice on how best to proceed....???
>
> Choose a built-in data set, such as Produc, change two states to ALASKA and HAWAII, drop the ones you chose from the map and any other states not on your map, and use that data set, putting the script needed at the head of your example. Then create the weights, initially with no links to ALASKA and HAWAII, and run the model. Next, edit the neighbour object in script form, and try again. Then anyone with your map object (you gave a link), your script, and the plm and splm packages can reproduce your problem.
>
> Do also run traceback() after errors, and post what you see.
>
> Roger
>
>> Many thanks again,
>> Stuart.
>>
>> -----Original Message-----
>> From: Roger Bivand <[hidden email]>
>> Sent: Thursday, 8 August, 2019 9:10 PM
>> To: Stuart Reece <[hidden email]>
>> Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson'
>> <[hidden email]>; 'R-sig-geo Mailing List'
>> <[hidden email]>
>> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
>> List
>>
>> On Thu, 8 Aug 2019, Stuart Reece wrote:
>>
>>> Thankyou so much Vijay and Barry.
>>
>> Yes, thanks!
>>
>>>
>>> Yes I have found that using these techniques I can add nb
>>> relationships nicely.
>>>
>>> And they card() nicely too.
>>>
>>> And when I nb2listw they transform well and when I test them
>>> can.be.simmed(x) they pass this test also.
>>>
>>> When I draw maps with them Like Roger Bivand showed in his examples
>>> for poly2nb – the maps draw just perfectly.
>>
>> This is flying on hope not knowledge. If the assigned values are not
>> within the 1:n set, there will be trouble - this is an S3 class with
>> no validity check. Further, there may be trouble if the length of
>> attr(.,
>> "region.id") is incorrect, but this cannot be checked unless you run
>> traceback() after the error in spgm(). We still need a fully reproducible example ...
>>
>> Roger
>>
>>> However when I use them in the spml or spgm regressions in package
>>> splm they fail with an error message usually about indexes being out
>>> of bounds.
>>>
>>> Running this code ….
>>>
>>> summary(spgm(PC1MI ~ PC1Rx * log(mrjmon),
>>> +              data=MIdf9dfF2, listw=MIdf9sflww,
>>> +              lag=TRUE, moments="fullweights", method="g2sls",
>>> +              model="random", spatial.error=TRUE))
>>>
>>> Generates this error:
>>>
>>> Error in x[, ii] : subscript out of bounds
>>>
>>> I feel that I must be missing something here but am not able to put
>>> my finger on what it is???
>>>
>>> Thanks so much again,
>>>
>>> Stuart.
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> From: Dr Stuart Reece <[hidden email]>
>>> Sent: Thursday, 8 August, 2019 2:01 PM
>>> To: 'Barry Rowlingson' <[hidden email]>; 'Stuart Reece'
>>> <[hidden email]>
>>> Cc: 'Vijay Lulla' <[hidden email]>; 'R-sig-geo Mailing List'
>>> <[hidden email]>
>>> Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
>>> List
>>> Importance: High
>>>
>>>
>>>
>>> Thanks Barry.
>>>
>>> That is so perfect and so super helpful!!!
>>>
>>> And if I want to add three areas to an area – say I want to add 6,7 and 8 to the area 10??
>>>
>>> Please may I have the syntax for that to avoid the integer error??
>>>
>>> Is this also the root of the error about not being the correct index??
>>>
>>> Many thanks again,
>>>
>>> Stuart.
>>>
>>>
>>>
>>> From: Barry Rowlingson [mailto:[hidden email]]
>>> Sent: Thursday, 8 August 2019 8:51 AM
>>> To: Stuart Reece
>>> Cc: Vijay Lulla; R-sig-geo Mailing List; Stuart Reece
>>> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
>>> List
>>>
>>>
>>>
>>> I recently answered a similar question on Stack Overflow where someone needed to add detached polygons to their connected network by connecting them to their nearest neighbour:
>>>
>>>
>>>
>>> https://stackoverflow.com/questions/57269254/how-to-impute-missing-ne
>>> i
>>> ghbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredire
>>> c
>>> t=1#comment101246065_57378930
>>>
>>>
>>>
>>> in short, you can treat a `nb` object like a list of vectors: nb[[i]]
>>> is a vector of indexes of objects connected to object `i`
>>>
>>>
>>>
>>> BUT you have to make sure you store integers:
>>>
>>>
>>>
>>> Here's a `nb` object from that question which in summary has this manyneighbours for each region:
>>>
>>>
>>>
>>>> card(nb)
>>> [1] 2 3 4 3 2 0 0
>>>
>>>
>>>
>>> lets set the 6th feature to be a neighbour of the first:
>>>
>>>
>>>
>>>> nb[[6]] = 1
>>>
>>>
>>>
>>> then uh-oh...
>>>
>>>
>>>
>>>> card(nb)
>>> Error in card(nb) :
>>>  INTEGER() can only be applied to a 'integer', not a 'double'
>>>
>>>
>>>
>>> same again only `as.integer`:
>>>
>>>
>>>
>>>> nb[[6]] = as.integer(1)
>>>
>>>
>>>
>>> and its happy:
>>>
>>>
>>>
>>>> card(nb)
>>> [1] 2 3 4 3 2 1 0
>>>
>>>
>>>
>>> if you want to set the nighbours of 6 to several features:
>>>
>>>
>>>
>>>> nb[[6]] = as.integer(c(1,2,3))
>>>> card(nb)
>>> [1] 2 3 4 3 2 3 0
>>>
>>>
>>>
>>> Barry
>>>
>>>
>>>
>>>
>>>
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> 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.
>> 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
>>
>>
>
> --
> 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
>
> --
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: Adding a Few Neighbour Relationships to a nb List

Thu, 08/08/2019 - 06:42
Thanks Roger.
Happy to do as you have requested.
Which server do you suggest??
Thanks again,
Stuart.

-----Original Message-----
From: Roger Bivand <[hidden email]>
Sent: Thursday, 8 August, 2019 9:35 PM
To: Stuart Reece <[hidden email]>
Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List

On Thu, 8 Aug 2019, Stuart Reece wrote:

> Thankyou for this advice Roger.
> Happy to provide my work for your review but I am not sure how.
> This list server has a limit of 50KB and these files are about 100MB....
> Love to assist but I would need advice on how best to proceed....???

Choose a built-in data set, such as Produc, change two states to ALASKA and HAWAII, drop the ones you chose from the map and any other states not on your map, and use that data set, putting the script needed at the head of your example. Then create the weights, initially with no links to ALASKA and HAWAII, and run the model. Next, edit the neighbour object in script form, and try again. Then anyone with your map object (you gave a link), your script, and the plm and splm packages can reproduce your problem.

Do also run traceback() after errors, and post what you see.

Roger

> Many thanks again,
> Stuart.
>
> -----Original Message-----
> From: Roger Bivand <[hidden email]>
> Sent: Thursday, 8 August, 2019 9:10 PM
> To: Stuart Reece <[hidden email]>
> Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson'
> <[hidden email]>; 'R-sig-geo Mailing List'
> <[hidden email]>
> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
> List
>
> On Thu, 8 Aug 2019, Stuart Reece wrote:
>
>> Thankyou so much Vijay and Barry.
>
> Yes, thanks!
>
>>
>> Yes I have found that using these techniques I can add nb
>> relationships nicely.
>>
>> And they card() nicely too.
>>
>> And when I nb2listw they transform well and when I test them
>> can.be.simmed(x) they pass this test also.
>>
>> When I draw maps with them Like Roger Bivand showed in his examples
>> for poly2nb – the maps draw just perfectly.
>
> This is flying on hope not knowledge. If the assigned values are not
> within the 1:n set, there will be trouble - this is an S3 class with
> no validity check. Further, there may be trouble if the length of
> attr(.,
> "region.id") is incorrect, but this cannot be checked unless you run
> traceback() after the error in spgm(). We still need a fully reproducible example ...
>
> Roger
>
>> However when I use them in the spml or spgm regressions in package
>> splm they fail with an error message usually about indexes being out
>> of bounds.
>>
>> Running this code ….
>>
>> summary(spgm(PC1MI ~ PC1Rx * log(mrjmon),
>> +              data=MIdf9dfF2, listw=MIdf9sflww,
>> +              lag=TRUE, moments="fullweights", method="g2sls",
>> +              model="random", spatial.error=TRUE))
>>
>> Generates this error:
>>
>> Error in x[, ii] : subscript out of bounds
>>
>> I feel that I must be missing something here but am not able to put
>> my finger on what it is???
>>
>> Thanks so much again,
>>
>> Stuart.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> From: Dr Stuart Reece <[hidden email]>
>> Sent: Thursday, 8 August, 2019 2:01 PM
>> To: 'Barry Rowlingson' <[hidden email]>; 'Stuart Reece'
>> <[hidden email]>
>> Cc: 'Vijay Lulla' <[hidden email]>; 'R-sig-geo Mailing List'
>> <[hidden email]>
>> Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
>> List
>> Importance: High
>>
>>
>>
>> Thanks Barry.
>>
>> That is so perfect and so super helpful!!!
>>
>> And if I want to add three areas to an area – say I want to add 6,7 and 8 to the area 10??
>>
>> Please may I have the syntax for that to avoid the integer error??
>>
>> Is this also the root of the error about not being the correct index??
>>
>> Many thanks again,
>>
>> Stuart.
>>
>>
>>
>> From: Barry Rowlingson [mailto:[hidden email]]
>> Sent: Thursday, 8 August 2019 8:51 AM
>> To: Stuart Reece
>> Cc: Vijay Lulla; R-sig-geo Mailing List; Stuart Reece
>> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
>> List
>>
>>
>>
>> I recently answered a similar question on Stack Overflow where someone needed to add detached polygons to their connected network by connecting them to their nearest neighbour:
>>
>>
>>
>> https://stackoverflow.com/questions/57269254/how-to-impute-missing-ne
>> i
>> ghbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredire
>> c
>> t=1#comment101246065_57378930
>>
>>
>>
>> in short, you can treat a `nb` object like a list of vectors: nb[[i]]
>> is a vector of indexes of objects connected to object `i`
>>
>>
>>
>> BUT you have to make sure you store integers:
>>
>>
>>
>> Here's a `nb` object from that question which in summary has this manyneighbours for each region:
>>
>>
>>
>>> card(nb)
>> [1] 2 3 4 3 2 0 0
>>
>>
>>
>> lets set the 6th feature to be a neighbour of the first:
>>
>>
>>
>>> nb[[6]] = 1
>>
>>
>>
>> then uh-oh...
>>
>>
>>
>>> card(nb)
>> Error in card(nb) :
>>  INTEGER() can only be applied to a 'integer', not a 'double'
>>
>>
>>
>> same again only `as.integer`:
>>
>>
>>
>>> nb[[6]] = as.integer(1)
>>
>>
>>
>> and its happy:
>>
>>
>>
>>> card(nb)
>> [1] 2 3 4 3 2 1 0
>>
>>
>>
>> if you want to set the nighbours of 6 to several features:
>>
>>
>>
>>> nb[[6]] = as.integer(c(1,2,3))
>>> card(nb)
>> [1] 2 3 4 3 2 3 0
>>
>>
>>
>> Barry
>>
>>
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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.
> 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
>
>
--
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

Re: Adding a Few Neighbour Relationships to a nb List

Thu, 08/08/2019 - 06:34
On Thu, 8 Aug 2019, Stuart Reece wrote:

> Thankyou for this advice Roger.
> Happy to provide my work for your review but I am not sure how.
> This list server has a limit of 50KB and these files are about 100MB....
> Love to assist but I would need advice on how best to proceed....???

Choose a built-in data set, such as Produc, change two states to ALASKA
and HAWAII, drop the ones you chose from the map and any other states not
on your map, and use that data set, putting the script needed at the head
of your example. Then create the weights, initially with no links to
ALASKA and HAWAII, and run the model. Next, edit the neighbour object in
script form, and try again. Then anyone with your map object (you gave a
link), your script, and the plm and splm packages can reproduce your
problem.

Do also run traceback() after errors, and post what you see.

Roger

> Many thanks again,
> Stuart.
>
> -----Original Message-----
> From: Roger Bivand <[hidden email]>
> Sent: Thursday, 8 August, 2019 9:10 PM
> To: Stuart Reece <[hidden email]>
> Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
>
> On Thu, 8 Aug 2019, Stuart Reece wrote:
>
>> Thankyou so much Vijay and Barry.
>
> Yes, thanks!
>
>>
>> Yes I have found that using these techniques I can add nb
>> relationships nicely.
>>
>> And they card() nicely too.
>>
>> And when I nb2listw they transform well and when I test them
>> can.be.simmed(x) they pass this test also.
>>
>> When I draw maps with them Like Roger Bivand showed in his examples
>> for poly2nb – the maps draw just perfectly.
>
> This is flying on hope not knowledge. If the assigned values are not within the 1:n set, there will be trouble - this is an S3 class with no validity check. Further, there may be trouble if the length of attr(.,
> "region.id") is incorrect, but this cannot be checked unless you run
> traceback() after the error in spgm(). We still need a fully reproducible example ...
>
> Roger
>
>> However when I use them in the spml or spgm regressions in package
>> splm they fail with an error message usually about indexes being out
>> of bounds.
>>
>> Running this code ….
>>
>> summary(spgm(PC1MI ~ PC1Rx * log(mrjmon),
>> +              data=MIdf9dfF2, listw=MIdf9sflww,
>> +              lag=TRUE, moments="fullweights", method="g2sls",
>> +              model="random", spatial.error=TRUE))
>>
>> Generates this error:
>>
>> Error in x[, ii] : subscript out of bounds
>>
>> I feel that I must be missing something here but am not able to put my
>> finger on what it is???
>>
>> Thanks so much again,
>>
>> Stuart.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> From: Dr Stuart Reece <[hidden email]>
>> Sent: Thursday, 8 August, 2019 2:01 PM
>> To: 'Barry Rowlingson' <[hidden email]>; 'Stuart Reece'
>> <[hidden email]>
>> Cc: 'Vijay Lulla' <[hidden email]>; 'R-sig-geo Mailing List'
>> <[hidden email]>
>> Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
>> List
>> Importance: High
>>
>>
>>
>> Thanks Barry.
>>
>> That is so perfect and so super helpful!!!
>>
>> And if I want to add three areas to an area – say I want to add 6,7 and 8 to the area 10??
>>
>> Please may I have the syntax for that to avoid the integer error??
>>
>> Is this also the root of the error about not being the correct index??
>>
>> Many thanks again,
>>
>> Stuart.
>>
>>
>>
>> From: Barry Rowlingson [mailto:[hidden email]]
>> Sent: Thursday, 8 August 2019 8:51 AM
>> To: Stuart Reece
>> Cc: Vijay Lulla; R-sig-geo Mailing List; Stuart Reece
>> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
>> List
>>
>>
>>
>> I recently answered a similar question on Stack Overflow where someone needed to add detached polygons to their connected network by connecting them to their nearest neighbour:
>>
>>
>>
>> https://stackoverflow.com/questions/57269254/how-to-impute-missing-nei
>> ghbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredirec
>> t=1#comment101246065_57378930
>>
>>
>>
>> in short, you can treat a `nb` object like a list of vectors: nb[[i]]
>> is a vector of indexes of objects connected to object `i`
>>
>>
>>
>> BUT you have to make sure you store integers:
>>
>>
>>
>> Here's a `nb` object from that question which in summary has this manyneighbours for each region:
>>
>>
>>
>>> card(nb)
>> [1] 2 3 4 3 2 0 0
>>
>>
>>
>> lets set the 6th feature to be a neighbour of the first:
>>
>>
>>
>>> nb[[6]] = 1
>>
>>
>>
>> then uh-oh...
>>
>>
>>
>>> card(nb)
>> Error in card(nb) :
>>  INTEGER() can only be applied to a 'integer', not a 'double'
>>
>>
>>
>> same again only `as.integer`:
>>
>>
>>
>>> nb[[6]] = as.integer(1)
>>
>>
>>
>> and its happy:
>>
>>
>>
>>> card(nb)
>> [1] 2 3 4 3 2 1 0
>>
>>
>>
>> if you want to set the nighbours of 6 to several features:
>>
>>
>>
>>> nb[[6]] = as.integer(c(1,2,3))
>>> card(nb)
>> [1] 2 3 4 3 2 3 0
>>
>>
>>
>> Barry
>>
>>
>>
>>
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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.
> 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
>
> --
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: Adding a Few Neighbour Relationships to a nb List

Thu, 08/08/2019 - 06:23
Thankyou for this advice Roger.
Happy to provide my work for your review but I am not sure how.
This list server has a limit of 50KB and these files are about 100MB....
Love to assist but I would need advice on how best to proceed....???
Many thanks again,
Stuart.

-----Original Message-----
From: Roger Bivand <[hidden email]>
Sent: Thursday, 8 August, 2019 9:10 PM
To: Stuart Reece <[hidden email]>
Cc: 'Dr Stuart Reece' <[hidden email]>; 'Barry Rowlingson' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List

On Thu, 8 Aug 2019, Stuart Reece wrote:

> Thankyou so much Vijay and Barry.

Yes, thanks!

>
> Yes I have found that using these techniques I can add nb
> relationships nicely.
>
> And they card() nicely too.
>
> And when I nb2listw they transform well and when I test them
> can.be.simmed(x) they pass this test also.
>
> When I draw maps with them Like Roger Bivand showed in his examples
> for poly2nb – the maps draw just perfectly.
This is flying on hope not knowledge. If the assigned values are not within the 1:n set, there will be trouble - this is an S3 class with no validity check. Further, there may be trouble if the length of attr(.,
"region.id") is incorrect, but this cannot be checked unless you run
traceback() after the error in spgm(). We still need a fully reproducible example ...

Roger

> However when I use them in the spml or spgm regressions in package
> splm they fail with an error message usually about indexes being out
> of bounds.
>
> Running this code ….
>
> summary(spgm(PC1MI ~ PC1Rx * log(mrjmon),
> +              data=MIdf9dfF2, listw=MIdf9sflww,
> +              lag=TRUE, moments="fullweights", method="g2sls",
> +              model="random", spatial.error=TRUE))
>
> Generates this error:
>
> Error in x[, ii] : subscript out of bounds
>
> I feel that I must be missing something here but am not able to put my
> finger on what it is???
>
> Thanks so much again,
>
> Stuart.
>
>
>
>
>
>
>
>
>
> From: Dr Stuart Reece <[hidden email]>
> Sent: Thursday, 8 August, 2019 2:01 PM
> To: 'Barry Rowlingson' <[hidden email]>; 'Stuart Reece'
> <[hidden email]>
> Cc: 'Vijay Lulla' <[hidden email]>; 'R-sig-geo Mailing List'
> <[hidden email]>
> Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
> List
> Importance: High
>
>
>
> Thanks Barry.
>
> That is so perfect and so super helpful!!!
>
> And if I want to add three areas to an area – say I want to add 6,7 and 8 to the area 10??
>
> Please may I have the syntax for that to avoid the integer error??
>
> Is this also the root of the error about not being the correct index??
>
> Many thanks again,
>
> Stuart.
>
>
>
> From: Barry Rowlingson [mailto:[hidden email]]
> Sent: Thursday, 8 August 2019 8:51 AM
> To: Stuart Reece
> Cc: Vijay Lulla; R-sig-geo Mailing List; Stuart Reece
> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
> List
>
>
>
> I recently answered a similar question on Stack Overflow where someone needed to add detached polygons to their connected network by connecting them to their nearest neighbour:
>
>
>
> https://stackoverflow.com/questions/57269254/how-to-impute-missing-nei
> ghbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredirec
> t=1#comment101246065_57378930
>
>
>
> in short, you can treat a `nb` object like a list of vectors: nb[[i]]
> is a vector of indexes of objects connected to object `i`
>
>
>
> BUT you have to make sure you store integers:
>
>
>
> Here's a `nb` object from that question which in summary has this manyneighbours for each region:
>
>
>
>> card(nb)
> [1] 2 3 4 3 2 0 0
>
>
>
> lets set the 6th feature to be a neighbour of the first:
>
>
>
>> nb[[6]] = 1
>
>
>
> then uh-oh...
>
>
>
>> card(nb)
> Error in card(nb) :
>  INTEGER() can only be applied to a 'integer', not a 'double'
>
>
>
> same again only `as.integer`:
>
>
>
>> nb[[6]] = as.integer(1)
>
>
>
> and its happy:
>
>
>
>> card(nb)
> [1] 2 3 4 3 2 1 0
>
>
>
> if you want to set the nighbours of 6 to several features:
>
>
>
>> nb[[6]] = as.integer(c(1,2,3))
>> card(nb)
> [1] 2 3 4 3 2 3 0
>
>
>
> Barry
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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.
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

Re: Adding a Few Neighbour Relationships to a nb List

Thu, 08/08/2019 - 06:10
On Thu, 8 Aug 2019, Stuart Reece wrote:

> Thankyou so much Vijay and Barry.

Yes, thanks!

>
> Yes I have found that using these techniques I can add nb relationships
> nicely.
>
> And they card() nicely too.
>
> And when I nb2listw they transform well and when I test them
> can.be.simmed(x) they pass this test also.
>
> When I draw maps with them Like Roger Bivand showed in his examples for
> poly2nb – the maps draw just perfectly. This is flying on hope not knowledge. If the assigned values are not
within the 1:n set, there will be trouble - this is an S3 class with no
validity check. Further, there may be trouble if the length of attr(.,
"region.id") is incorrect, but this cannot be checked unless you run
traceback() after the error in spgm(). We still need a fully reproducible
example ...

Roger

> However when I use them in the spml or spgm regressions in package splm
> they fail with an error message usually about indexes being out of
> bounds.
>
> Running this code ….
>
> summary(spgm(PC1MI ~ PC1Rx * log(mrjmon),
> +              data=MIdf9dfF2, listw=MIdf9sflww,
> +              lag=TRUE, moments="fullweights", method="g2sls",
> +              model="random", spatial.error=TRUE))
>
> Generates this error:
>
> Error in x[, ii] : subscript out of bounds
>
> I feel that I must be missing something here but am not able to put my
> finger on what it is???
>
> Thanks so much again,
>
> Stuart.
>
>
>
>
>
>
>
>
>
> From: Dr Stuart Reece <[hidden email]>
> Sent: Thursday, 8 August, 2019 2:01 PM
> To: 'Barry Rowlingson' <[hidden email]>; 'Stuart Reece' <[hidden email]>
> Cc: 'Vijay Lulla' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
> Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
> Importance: High
>
>
>
> Thanks Barry.
>
> That is so perfect and so super helpful!!!
>
> And if I want to add three areas to an area – say I want to add 6,7 and 8 to the area 10??
>
> Please may I have the syntax for that to avoid the integer error??
>
> Is this also the root of the error about not being the correct index??
>
> Many thanks again,
>
> Stuart.
>
>
>
> From: Barry Rowlingson [mailto:[hidden email]]
> Sent: Thursday, 8 August 2019 8:51 AM
> To: Stuart Reece
> Cc: Vijay Lulla; R-sig-geo Mailing List; Stuart Reece
> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
>
>
>
> I recently answered a similar question on Stack Overflow where someone needed to add detached polygons to their connected network by connecting them to their nearest neighbour:
>
>
>
> https://stackoverflow.com/questions/57269254/how-to-impute-missing-neighbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredirect=1#comment101246065_57378930
>
>
>
> in short, you can treat a `nb` object like a list of vectors: nb[[i]] is a vector of indexes of objects connected to object `i`
>
>
>
> BUT you have to make sure you store integers:
>
>
>
> Here's a `nb` object from that question which in summary has this manyneighbours for each region:
>
>
>
>> card(nb)
> [1] 2 3 4 3 2 0 0
>
>
>
> lets set the 6th feature to be a neighbour of the first:
>
>
>
>> nb[[6]] = 1
>
>
>
> then uh-oh...
>
>
>
>> card(nb)
> Error in card(nb) :
>  INTEGER() can only be applied to a 'integer', not a 'double'
>
>
>
> same again only `as.integer`:
>
>
>
>> nb[[6]] = as.integer(1)
>
>
>
> and its happy:
>
>
>
>> card(nb)
> [1] 2 3 4 3 2 1 0
>
>
>
> if you want to set the nighbours of 6 to several features:
>
>
>
>> nb[[6]] = as.integer(c(1,2,3))
>> card(nb)
> [1] 2 3 4 3 2 3 0
>
>
>
> Barry
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> 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.
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: Adding a Few Neighbour Relationships to a nb List

Thu, 08/08/2019 - 05:48
Thankyou so much Vijay and Barry.

Yes I have found that using these techniques I can add nb relationships nicely.

And they card() nicely too.

And when I nb2listw they transform well and when I test them can.be.simmed(x) they pass this test also.

 

When I draw maps with them Like Roger Bivand showed in his examples for poly2nb – the maps draw just perfectly.

 

However when I use them in the spml or spgm regressions in package splm they fail with an error message usually about indexes being out of bounds.

 

Running this code ….

 

summary(spgm(PC1MI ~ PC1Rx * log(mrjmon),

+              data=MIdf9dfF2, listw=MIdf9sflww,

+              lag=TRUE, moments="fullweights", method="g2sls",

+              model="random", spatial.error=TRUE))

 

Generates this error:

Error in x[, ii] : subscript out of bounds

 

I feel that I must be missing something here but am not able to put my finger on what it is???

Thanks so much again,

Stuart.

 

 

 

 

From: Dr Stuart Reece <[hidden email]>
Sent: Thursday, 8 August, 2019 2:01 PM
To: 'Barry Rowlingson' <[hidden email]>; 'Stuart Reece' <[hidden email]>
Cc: 'Vijay Lulla' <[hidden email]>; 'R-sig-geo Mailing List' <[hidden email]>
Subject: RE: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
Importance: High

 

Thanks Barry.

That is so perfect and so super helpful!!!

And if I want to add three areas to an area – say I want to add 6,7 and 8 to the area 10??

Please may I have the syntax for that to avoid the integer error??

Is this also the root of the error about not being the correct index??

Many thanks again,

Stuart.

 

From: Barry Rowlingson [mailto:[hidden email]]
Sent: Thursday, 8 August 2019 8:51 AM
To: Stuart Reece
Cc: Vijay Lulla; R-sig-geo Mailing List; Stuart Reece
Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List

 

I recently answered a similar question on Stack Overflow where someone needed to add detached polygons to their connected network by connecting them to their nearest neighbour:

 

https://stackoverflow.com/questions/57269254/how-to-impute-missing-neighbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredirect=1#comment101246065_57378930

 

in short, you can treat a `nb` object like a list of vectors: nb[[i]] is a vector of indexes of objects connected to object `i`

 

BUT you have to make sure you store integers:

 

Here's a `nb` object from that question which in summary has this manyneighbours for each region:

 

> card(nb)
[1] 2 3 4 3 2 0 0

 

lets set the 6th feature to be a neighbour of the first:

 

> nb[[6]] = 1

 

then uh-oh...

 

> card(nb)
Error in card(nb) :
  INTEGER() can only be applied to a 'integer', not a 'double'

 

same again only `as.integer`:

 

> nb[[6]] = as.integer(1)

 

and its happy:

 

> card(nb)
[1] 2 3 4 3 2 1 0

 

if you want to set the nighbours of 6 to several features:

 

> nb[[6]] = as.integer(c(1,2,3))
> card(nb)
[1] 2 3 4 3 2 3 0

 

Barry

 

 


        [[alternative HTML version deleted]]

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

Re: Adding a Few Neighbour Relationships to a nb List

Thu, 08/08/2019 - 05:02
Wrap the assigned vector in `as.integer`:

Works:
> col.gal.nb[[1]] = as.integer(c(4,5))
> card(col.gal.nb)
 [1]  2  3  4  4  7  2  4  6  8  4  5  6  4  6  6  7  3  4  3 10  3  6  3
 7  7
[26]  6  4  9  7  4  2  4  4  4  7  5  6  6  2  5  3  2  6  5  4  2  2  4  3


Fails:

> col.gal.nb[[1]] = c(4,5)
> card(col.gal.nb)
Error in card(col.gal.nb) :
  INTEGER() can only be applied to a 'integer', not a 'double'

Barry


On Thu, Aug 8, 2019 at 5:01 AM Dr Stuart Reece <[hidden email]>
wrote:

> Thanks Barry.
>
> That is so perfect and so super helpful!!!
>
> And if I want to add three areas to an area – say I want to add 6,7 and 8
> to the area 10??
>
> Please may I have the syntax for that to avoid the integer error??
>
> Is this also the root of the error about not being the correct index??
>
> Many thanks again,
>
> Stuart.
>
>
>
> *From:* Barry Rowlingson [mailto:[hidden email]]
> *Sent:* Thursday, 8 August 2019 8:51 AM
> *To:* Stuart Reece
> *Cc:* Vijay Lulla; R-sig-geo Mailing List; Stuart Reece
> *Subject:* Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb
> List
>
>
>
> I recently answered a similar question on Stack Overflow where someone
> needed to add detached polygons to their connected network by connecting
> them to their nearest neighbour:
>
>
>
>
> https://stackoverflow.com/questions/57269254/how-to-impute-missing-neighbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredirect=1#comment101246065_57378930
>
>
>
> in short, you can treat a `nb` object like a list of vectors: nb[[i]] is a
> vector of indexes of objects connected to object `i`
>
>
>
> BUT you have to make sure you store integers:
>
>
>
> Here's a `nb` object from that question which in summary has this
> manyneighbours for each region:
>
>
>
> > card(nb)
> [1] 2 3 4 3 2 0 0
>
>
>
> lets set the 6th feature to be a neighbour of the first:
>
>
>
> > nb[[6]] = 1
>
>
>
> then uh-oh...
>
>
>
> > card(nb)
> Error in card(nb) :
>   INTEGER() can only be applied to a 'integer', not a 'double'
>
>
>
> same again only `as.integer`:
>
>
>
> > nb[[6]] = as.integer(1)
>
>
>
> and its happy:
>
>
>
> > card(nb)
> [1] 2 3 4 3 2 1 0
>
>
>
> if you want to set the nighbours of 6 to several features:
>
>
>
> > nb[[6]] = as.integer(c(1,2,3))
> > card(nb)
> [1] 2 3 4 3 2 3 0
>
>
>
> Barry
>
>
>
>
>
>
>
> On Wed, Aug 7, 2019 at 10:25 PM Stuart Reece <[hidden email]>
> wrote:
>
> Thanks Vijay.
>
> Big help.
>
> I will go through the recommended chapter in detail.
>
> Normally for a list I can just use single square brackets like this [] to
> access the elements and change them by assignment.
>
> But I cannot quite figure out how to do this with the nb lists.
>
>
>
> Doing this with [[]] works really well to create nicely corrected graphs.
>
>
>
> But fails due it “out of bounds index errors” with regression equations????
>
>
>
> I find this ever so confusing….???
>
>
>
> Thankyou so much again,
>
>
>
> Stuart Reece.
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> From: Vijay Lulla <[hidden email]>
> Sent: Thursday, 8 August, 2019 1:53 AM
> To: Stuart Reece <[hidden email]>
> Cc: Roger Bivand <[hidden email]>; R-sig-geo Mailing List <
> [hidden email]>
> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
>
>
>
> Maybe
> https://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html
> can help with all your questions.
> https://cran.r-project.org/web/packages/spdep/vignettes/nb_sf.html
> contains a little more detail about nb structure.  Finally, I encourage you
> to use `str` to study the structure of R objects.
>
> HTH,
>
> Vijay.
>
>
>
> On Wed, Aug 7, 2019 at 10:44 AM Stuart Reece <[hidden email]
> <mailto:[hidden email]> > wrote:
>
> Dear R-Sig-Geo list,
>
>
>
> I was wondering if it might be possible please to request assistance with
> adding some nb relationships to a .nb.gal list composed either by GeoDa or
> poly2nb in R????
>
>
>
> The shapefile at this URL
> <
> https://www.samhsa.gov/data/report/2014-2016-nsduh-substate-region-shapefil
> e>  divides USA into 395 substate regions.  For health and demographical
> reasons it is important to include both Hawaii and Alaska in the
> spatiotemporal analysis so I want to introduce these into the Southeast
> coast of California and the Pacific northwest respectively.
>
>
>
> This is just as Giovanni Millo added in spatial relationships for Sicily
> across the Strait of Messina for splm on page 7 of the splm pdf.
>
>
>
> I found edit.nb in spdep and operated it just as described in the
> instructions and here
> <https://github.com/r-spatial/spdep/blob/master/man/edit.nb.Rd> .  It
> crashed RStudio many times but ran well in R3.6.1.  However even though I
> assigned it to a new object it did not save well.  Although when I plotted
> the dxxx file as the difference between the old and modified files it
> plotted the changes beautifully in red and black respectively when plotted
> by themselves it introduced many long distance extraneous relationships.
> To
> get the edited nb list file out of R 3.6.1 and into RStudio I saved it as
> an
> RDS file.  However when opened in RStudio it was grossly erroneous and
> included extraneous links from Hawaii to Boston and New York.  When I
> opened
> the file in RStudio it again introduced these extraneous links.
>
>
>
> Saving it as a further new object in R 3.6.1 did not remedy these
> difficulties.
>
>
>
> The other problem I have is that the spdep poly2nb function excludes
> Richmond, an island off the southern tip of Long Island near New York as it
> is an island.  Also one of the areas - Region 10 in Washington DC - is also
> excluded for reasons of which I am unsure.
>
>
>
> I found some code here to just patch single areas
> <https://stat.ethz.ch/pipermail/r-sig-geo/2006-June/001073.html>  like
> this
> but when I run it, it throws an integer error
>
> "  INTEGER() can only be applied to a 'integer', not a 'double'
>
>
>
> No combination of bracketing around subscripts helps or works at all.  The
> link mentioned has these statements in it
>
>
>
> nb[[ij[1]]] <- sort(unique(c(nb[[ij[1]]], ij[2])))
> nb[[ij[2]]] <- sort(unique(c(nb[[ij[2]]], ij[1])))
>
>
>
> which makes me think that I should insert a vector "  c(i,j) "  where
> indicated.  Even using "c(as.integer(i),as.integer(j)) "
>
> or " as.integer(c(I,j)) "   doesn't work and still gives rise to the same
> error.
>
>
>
> I am sure I am not the only one to have encountered such difficulties but I
> have really tried everything I can think of.
>
>
>
> The other thing I would really like is some clear instructions as to the
> true underlying structure of the nb list.  If I could clearly understand
> this then I could just go into the affected lines of the list of lists and
> edit them directly.
>
>
>
> However I am quite unable to find any clear description of its structure on
> line.
>
>
>
> Similarly I cannot find the source code for drop.links online to try to
> translate this code into add.links directly, as was also suggested.
>
>
>
> But such a function would I think be enormously helpful and of invaluable
> assistance for final editing.
>
>
>
> Thankyou ever so much in advance for your kind and gracious assistance.
>
>
>
> Yours sincerely,
>
>
>
> Stuart Reece.
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email] <mailto:[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: Adding a Few Neighbour Relationships to a nb List

Wed, 08/07/2019 - 18:11
In addition to Barry's recommendation/suggestion I will just add that I
always use the suffix L to numbers so that I have integers whenever I need
them.  Here's an example:
> str(1)
 num 1
> str(1L)
 int 1
> str(c(1,2,3))
 num [1:3] 1 2 3
> str(c(1L,2L,3L))
 int [1:3] 1 2 3
>

This isn't possible when you're indexing based on values extracted from
other lists/vectors in which case you have to use `as.integer` to force
integers.  Again, I find `str` to be invaluable in figuring all this out.
HTH,
Vijay.

On Wed, Aug 7, 2019 at 6:50 PM Barry Rowlingson <[hidden email]>
wrote:

> I recently answered a similar question on Stack Overflow where someone
> needed to add detached polygons to their connected network by connecting
> them to their nearest neighbour:
>
>
> https://stackoverflow.com/questions/57269254/how-to-impute-missing-neighbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredirect=1#comment101246065_57378930
>
> in short, you can treat a `nb` object like a list of vectors: nb[[i]] is a
> vector of indexes of objects connected to object `i`
>
> BUT you have to make sure you store integers:
>
> Here's a `nb` object from that question which in summary has this
> manyneighbours for each region:
>
> > card(nb)
> [1] 2 3 4 3 2 0 0
>
> lets set the 6th feature to be a neighbour of the first:
>
> > nb[[6]] = 1
>
> then uh-oh...
>
> > card(nb)
> Error in card(nb) :
>   INTEGER() can only be applied to a 'integer', not a 'double'
>
> same again only `as.integer`:
>
> > nb[[6]] = as.integer(1)
>
> and its happy:
>
> > card(nb)
> [1] 2 3 4 3 2 1 0
>
> if you want to set the nighbours of 6 to several features:
>
> > nb[[6]] = as.integer(c(1,2,3))
> > card(nb)
> [1] 2 3 4 3 2 3 0
>
> Barry
>
>
>
> On Wed, Aug 7, 2019 at 10:25 PM Stuart Reece <[hidden email]>
> wrote:
>
>> Thanks Vijay.
>>
>> Big help.
>>
>> I will go through the recommended chapter in detail.
>>
>> Normally for a list I can just use single square brackets like this [] to
>> access the elements and change them by assignment.
>>
>> But I cannot quite figure out how to do this with the nb lists.
>>
>>
>>
>> Doing this with [[]] works really well to create nicely corrected graphs.
>>
>>
>>
>> But fails due it “out of bounds index errors” with regression
>> equations????
>>
>>
>>
>> I find this ever so confusing….???
>>
>>
>>
>> Thankyou so much again,
>>
>>
>>
>> Stuart Reece.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> From: Vijay Lulla <[hidden email]>
>> Sent: Thursday, 8 August, 2019 1:53 AM
>> To: Stuart Reece <[hidden email]>
>> Cc: Roger Bivand <[hidden email]>; R-sig-geo Mailing List <
>> [hidden email]>
>> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
>>
>>
>>
>> Maybe
>> https://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html
>> can help with all your questions.
>> https://cran.r-project.org/web/packages/spdep/vignettes/nb_sf.html
>> contains a little more detail about nb structure.  Finally, I encourage you
>> to use `str` to study the structure of R objects.
>>
>> HTH,
>>
>> Vijay.
>>
>>
>>
>> On Wed, Aug 7, 2019 at 10:44 AM Stuart Reece <[hidden email]
>> <mailto:[hidden email]> > wrote:
>>
>> Dear R-Sig-Geo list,
>>
>>
>>
>> I was wondering if it might be possible please to request assistance with
>> adding some nb relationships to a .nb.gal list composed either by GeoDa or
>> poly2nb in R????
>>
>>
>>
>> The shapefile at this URL
>> <
>> https://www.samhsa.gov/data/report/2014-2016-nsduh-substate-region-shapefil
>> e>  divides USA into 395 substate regions.  For health and demographical
>> reasons it is important to include both Hawaii and Alaska in the
>> spatiotemporal analysis so I want to introduce these into the Southeast
>> coast of California and the Pacific northwest respectively.
>>
>>
>>
>> This is just as Giovanni Millo added in spatial relationships for Sicily
>> across the Strait of Messina for splm on page 7 of the splm pdf.
>>
>>
>>
>> I found edit.nb in spdep and operated it just as described in the
>> instructions and here
>> <https://github.com/r-spatial/spdep/blob/master/man/edit.nb.Rd> .  It
>> crashed RStudio many times but ran well in R3.6.1.  However even though I
>> assigned it to a new object it did not save well.  Although when I plotted
>> the dxxx file as the difference between the old and modified files it
>> plotted the changes beautifully in red and black respectively when plotted
>> by themselves it introduced many long distance extraneous relationships.
>> To
>> get the edited nb list file out of R 3.6.1 and into RStudio I saved it as
>> an
>> RDS file.  However when opened in RStudio it was grossly erroneous and
>> included extraneous links from Hawaii to Boston and New York.  When I
>> opened
>> the file in RStudio it again introduced these extraneous links.
>>
>>
>>
>> Saving it as a further new object in R 3.6.1 did not remedy these
>> difficulties.
>>
>>
>>
>> The other problem I have is that the spdep poly2nb function excludes
>> Richmond, an island off the southern tip of Long Island near New York as
>> it
>> is an island.  Also one of the areas - Region 10 in Washington DC - is
>> also
>> excluded for reasons of which I am unsure.
>>
>>
>>
>> I found some code here to just patch single areas
>> <https://stat.ethz.ch/pipermail/r-sig-geo/2006-June/001073.html>  like
>> this
>> but when I run it, it throws an integer error
>>
>> "  INTEGER() can only be applied to a 'integer', not a 'double'
>>
>>
>>
>> No combination of bracketing around subscripts helps or works at all.  The
>> link mentioned has these statements in it
>>
>>
>>
>> nb[[ij[1]]] <- sort(unique(c(nb[[ij[1]]], ij[2])))
>> nb[[ij[2]]] <- sort(unique(c(nb[[ij[2]]], ij[1])))
>>
>>
>>
>> which makes me think that I should insert a vector "  c(i,j) "  where
>> indicated.  Even using "c(as.integer(i),as.integer(j)) "
>>
>> or " as.integer(c(I,j)) "   doesn't work and still gives rise to the same
>> error.
>>
>>
>>
>> I am sure I am not the only one to have encountered such difficulties but
>> I
>> have really tried everything I can think of.
>>
>>
>>
>> The other thing I would really like is some clear instructions as to the
>> true underlying structure of the nb list.  If I could clearly understand
>> this then I could just go into the affected lines of the list of lists and
>> edit them directly.
>>
>>
>>
>> However I am quite unable to find any clear description of its structure
>> on
>> line.
>>
>>
>>
>> Similarly I cannot find the source code for drop.links online to try to
>> translate this code into add.links directly, as was also suggested.
>>
>>
>>
>> But such a function would I think be enormously helpful and of invaluable
>> assistance for final editing.
>>
>>
>>
>> Thankyou ever so much in advance for your kind and gracious assistance.
>>
>>
>>
>> Yours sincerely,
>>
>>
>>
>> Stuart Reece.
>>
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email] <mailto:[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: Adding a Few Neighbour Relationships to a nb List

Wed, 08/07/2019 - 17:50
I recently answered a similar question on Stack Overflow where someone
needed to add detached polygons to their connected network by connecting
them to their nearest neighbour:

https://stackoverflow.com/questions/57269254/how-to-impute-missing-neighbours-of-a-spatial-weight-matrix-queen-contiguity/57378930?noredirect=1#comment101246065_57378930

in short, you can treat a `nb` object like a list of vectors: nb[[i]] is a
vector of indexes of objects connected to object `i`

BUT you have to make sure you store integers:

Here's a `nb` object from that question which in summary has this
manyneighbours for each region:

> card(nb)
[1] 2 3 4 3 2 0 0

lets set the 6th feature to be a neighbour of the first:

> nb[[6]] = 1

then uh-oh...

> card(nb)
Error in card(nb) :
  INTEGER() can only be applied to a 'integer', not a 'double'

same again only `as.integer`:

> nb[[6]] = as.integer(1)

and its happy:

> card(nb)
[1] 2 3 4 3 2 1 0

if you want to set the nighbours of 6 to several features:

> nb[[6]] = as.integer(c(1,2,3))
> card(nb)
[1] 2 3 4 3 2 3 0

Barry



On Wed, Aug 7, 2019 at 10:25 PM Stuart Reece <[hidden email]>
wrote:

> Thanks Vijay.
>
> Big help.
>
> I will go through the recommended chapter in detail.
>
> Normally for a list I can just use single square brackets like this [] to
> access the elements and change them by assignment.
>
> But I cannot quite figure out how to do this with the nb lists.
>
>
>
> Doing this with [[]] works really well to create nicely corrected graphs.
>
>
>
> But fails due it “out of bounds index errors” with regression equations????
>
>
>
> I find this ever so confusing….???
>
>
>
> Thankyou so much again,
>
>
>
> Stuart Reece.
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> From: Vijay Lulla <[hidden email]>
> Sent: Thursday, 8 August, 2019 1:53 AM
> To: Stuart Reece <[hidden email]>
> Cc: Roger Bivand <[hidden email]>; R-sig-geo Mailing List <
> [hidden email]>
> Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List
>
>
>
> Maybe
> https://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html
> can help with all your questions.
> https://cran.r-project.org/web/packages/spdep/vignettes/nb_sf.html
> contains a little more detail about nb structure.  Finally, I encourage you
> to use `str` to study the structure of R objects.
>
> HTH,
>
> Vijay.
>
>
>
> On Wed, Aug 7, 2019 at 10:44 AM Stuart Reece <[hidden email]
> <mailto:[hidden email]> > wrote:
>
> Dear R-Sig-Geo list,
>
>
>
> I was wondering if it might be possible please to request assistance with
> adding some nb relationships to a .nb.gal list composed either by GeoDa or
> poly2nb in R????
>
>
>
> The shapefile at this URL
> <
> https://www.samhsa.gov/data/report/2014-2016-nsduh-substate-region-shapefil
> e>  divides USA into 395 substate regions.  For health and demographical
> reasons it is important to include both Hawaii and Alaska in the
> spatiotemporal analysis so I want to introduce these into the Southeast
> coast of California and the Pacific northwest respectively.
>
>
>
> This is just as Giovanni Millo added in spatial relationships for Sicily
> across the Strait of Messina for splm on page 7 of the splm pdf.
>
>
>
> I found edit.nb in spdep and operated it just as described in the
> instructions and here
> <https://github.com/r-spatial/spdep/blob/master/man/edit.nb.Rd> .  It
> crashed RStudio many times but ran well in R3.6.1.  However even though I
> assigned it to a new object it did not save well.  Although when I plotted
> the dxxx file as the difference between the old and modified files it
> plotted the changes beautifully in red and black respectively when plotted
> by themselves it introduced many long distance extraneous relationships.
> To
> get the edited nb list file out of R 3.6.1 and into RStudio I saved it as
> an
> RDS file.  However when opened in RStudio it was grossly erroneous and
> included extraneous links from Hawaii to Boston and New York.  When I
> opened
> the file in RStudio it again introduced these extraneous links.
>
>
>
> Saving it as a further new object in R 3.6.1 did not remedy these
> difficulties.
>
>
>
> The other problem I have is that the spdep poly2nb function excludes
> Richmond, an island off the southern tip of Long Island near New York as it
> is an island.  Also one of the areas - Region 10 in Washington DC - is also
> excluded for reasons of which I am unsure.
>
>
>
> I found some code here to just patch single areas
> <https://stat.ethz.ch/pipermail/r-sig-geo/2006-June/001073.html>  like
> this
> but when I run it, it throws an integer error
>
> "  INTEGER() can only be applied to a 'integer', not a 'double'
>
>
>
> No combination of bracketing around subscripts helps or works at all.  The
> link mentioned has these statements in it
>
>
>
> nb[[ij[1]]] <- sort(unique(c(nb[[ij[1]]], ij[2])))
> nb[[ij[2]]] <- sort(unique(c(nb[[ij[2]]], ij[1])))
>
>
>
> which makes me think that I should insert a vector "  c(i,j) "  where
> indicated.  Even using "c(as.integer(i),as.integer(j)) "
>
> or " as.integer(c(I,j)) "   doesn't work and still gives rise to the same
> error.
>
>
>
> I am sure I am not the only one to have encountered such difficulties but I
> have really tried everything I can think of.
>
>
>
> The other thing I would really like is some clear instructions as to the
> true underlying structure of the nb list.  If I could clearly understand
> this then I could just go into the affected lines of the list of lists and
> edit them directly.
>
>
>
> However I am quite unable to find any clear description of its structure on
> line.
>
>
>
> Similarly I cannot find the source code for drop.links online to try to
> translate this code into add.links directly, as was also suggested.
>
>
>
> But such a function would I think be enormously helpful and of invaluable
> assistance for final editing.
>
>
>
> Thankyou ever so much in advance for your kind and gracious assistance.
>
>
>
> Yours sincerely,
>
>
>
> Stuart Reece.
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email] <mailto:[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: Adding a Few Neighbour Relationships to a nb List

Wed, 08/07/2019 - 16:24
Thanks Vijay.

Big help.

I will go through the recommended chapter in detail.

Normally for a list I can just use single square brackets like this [] to access the elements and change them by assignment.

But I cannot quite figure out how to do this with the nb lists.

 

Doing this with [[]] works really well to create nicely corrected graphs.

 

But fails due it “out of bounds index errors” with regression equations????

 

I find this ever so confusing….???

 

Thankyou so much again,

 

Stuart Reece.

 

 

 

 

 

 

 

 

 

 

From: Vijay Lulla <[hidden email]>
Sent: Thursday, 8 August, 2019 1:53 AM
To: Stuart Reece <[hidden email]>
Cc: Roger Bivand <[hidden email]>; R-sig-geo Mailing List <[hidden email]>
Subject: Re: [R-sig-Geo] Adding a Few Neighbour Relationships to a nb List

 

Maybe https://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html can help with all your questions. https://cran.r-project.org/web/packages/spdep/vignettes/nb_sf.html contains a little more detail about nb structure.  Finally, I encourage you to use `str` to study the structure of R objects.

HTH,

Vijay.

 

On Wed, Aug 7, 2019 at 10:44 AM Stuart Reece <[hidden email] <mailto:[hidden email]> > wrote:

Dear R-Sig-Geo list,



I was wondering if it might be possible please to request assistance with
adding some nb relationships to a .nb.gal list composed either by GeoDa or
poly2nb in R????



The shapefile at this URL
<https://www.samhsa.gov/data/report/2014-2016-nsduh-substate-region-shapefil
e>  divides USA into 395 substate regions.  For health and demographical
reasons it is important to include both Hawaii and Alaska in the
spatiotemporal analysis so I want to introduce these into the Southeast
coast of California and the Pacific northwest respectively.



This is just as Giovanni Millo added in spatial relationships for Sicily
across the Strait of Messina for splm on page 7 of the splm pdf.



I found edit.nb in spdep and operated it just as described in the
instructions and here
<https://github.com/r-spatial/spdep/blob/master/man/edit.nb.Rd> .  It
crashed RStudio many times but ran well in R3.6.1.  However even though I
assigned it to a new object it did not save well.  Although when I plotted
the dxxx file as the difference between the old and modified files it
plotted the changes beautifully in red and black respectively when plotted
by themselves it introduced many long distance extraneous relationships.  To
get the edited nb list file out of R 3.6.1 and into RStudio I saved it as an
RDS file.  However when opened in RStudio it was grossly erroneous and
included extraneous links from Hawaii to Boston and New York.  When I opened
the file in RStudio it again introduced these extraneous links.



Saving it as a further new object in R 3.6.1 did not remedy these
difficulties.



The other problem I have is that the spdep poly2nb function excludes
Richmond, an island off the southern tip of Long Island near New York as it
is an island.  Also one of the areas - Region 10 in Washington DC - is also
excluded for reasons of which I am unsure.



I found some code here to just patch single areas
<https://stat.ethz.ch/pipermail/r-sig-geo/2006-June/001073.html>  like this
but when I run it, it throws an integer error

"  INTEGER() can only be applied to a 'integer', not a 'double'  



No combination of bracketing around subscripts helps or works at all.  The
link mentioned has these statements in it



nb[[ij[1]]] <- sort(unique(c(nb[[ij[1]]], ij[2])))
nb[[ij[2]]] <- sort(unique(c(nb[[ij[2]]], ij[1])))



which makes me think that I should insert a vector "  c(i,j) "  where
indicated.  Even using "c(as.integer(i),as.integer(j)) "

or " as.integer(c(I,j)) "   doesn't work and still gives rise to the same
error.



I am sure I am not the only one to have encountered such difficulties but I
have really tried everything I can think of.



The other thing I would really like is some clear instructions as to the
true underlying structure of the nb list.  If I could clearly understand
this then I could just go into the affected lines of the list of lists and
edit them directly.



However I am quite unable to find any clear description of its structure on
line.



Similarly I cannot find the source code for drop.links online to try to
translate this code into add.links directly, as was also suggested.



But such a function would I think be enormously helpful and of invaluable
assistance for final editing.



Thankyou ever so much in advance for your kind and gracious assistance.



Yours sincerely,



Stuart Reece.


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email] <mailto:[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: Adding a Few Neighbour Relationships to a nb List

Wed, 08/07/2019 - 10:53
Maybe https://cran.r-project.org/web/packages/spdep/vignettes/nb_igraph.html
can help with all your questions.
https://cran.r-project.org/web/packages/spdep/vignettes/nb_sf.html contains
a little more detail about nb structure.  Finally, I encourage you to use
`str` to study the structure of R objects.
HTH,
Vijay.

On Wed, Aug 7, 2019 at 10:44 AM Stuart Reece <[hidden email]>
wrote:

> Dear R-Sig-Geo list,
>
>
>
> I was wondering if it might be possible please to request assistance with
> adding some nb relationships to a .nb.gal list composed either by GeoDa or
> poly2nb in R????
>
>
>
> The shapefile at this URL
> <
> https://www.samhsa.gov/data/report/2014-2016-nsduh-substate-region-shapefil
> e>  divides USA into 395 substate regions.  For health and demographical
> reasons it is important to include both Hawaii and Alaska in the
> spatiotemporal analysis so I want to introduce these into the Southeast
> coast of California and the Pacific northwest respectively.
>
>
>
> This is just as Giovanni Millo added in spatial relationships for Sicily
> across the Strait of Messina for splm on page 7 of the splm pdf.
>
>
>
> I found edit.nb in spdep and operated it just as described in the
> instructions and here
> <https://github.com/r-spatial/spdep/blob/master/man/edit.nb.Rd> .  It
> crashed RStudio many times but ran well in R3.6.1.  However even though I
> assigned it to a new object it did not save well.  Although when I plotted
> the dxxx file as the difference between the old and modified files it
> plotted the changes beautifully in red and black respectively when plotted
> by themselves it introduced many long distance extraneous relationships.
> To
> get the edited nb list file out of R 3.6.1 and into RStudio I saved it as
> an
> RDS file.  However when opened in RStudio it was grossly erroneous and
> included extraneous links from Hawaii to Boston and New York.  When I
> opened
> the file in RStudio it again introduced these extraneous links.
>
>
>
> Saving it as a further new object in R 3.6.1 did not remedy these
> difficulties.
>
>
>
> The other problem I have is that the spdep poly2nb function excludes
> Richmond, an island off the southern tip of Long Island near New York as it
> is an island.  Also one of the areas - Region 10 in Washington DC - is also
> excluded for reasons of which I am unsure.
>
>
>
> I found some code here to just patch single areas
> <https://stat.ethz.ch/pipermail/r-sig-geo/2006-June/001073.html>  like
> this
> but when I run it, it throws an integer error
>
> "  INTEGER() can only be applied to a 'integer', not a 'double'
>
>
>
> No combination of bracketing around subscripts helps or works at all.  The
> link mentioned has these statements in it
>
>
>
> nb[[ij[1]]] <- sort(unique(c(nb[[ij[1]]], ij[2])))
> nb[[ij[2]]] <- sort(unique(c(nb[[ij[2]]], ij[1])))
>
>
>
> which makes me think that I should insert a vector "  c(i,j) "  where
> indicated.  Even using "c(as.integer(i),as.integer(j)) "
>
> or " as.integer(c(I,j)) "   doesn't work and still gives rise to the same
> error.
>
>
>
> I am sure I am not the only one to have encountered such difficulties but I
> have really tried everything I can think of.
>
>
>
> The other thing I would really like is some clear instructions as to the
> true underlying structure of the nb list.  If I could clearly understand
> this then I could just go into the affected lines of the list of lists and
> edit them directly.
>
>
>
> However I am quite unable to find any clear description of its structure on
> line.
>
>
>
> Similarly I cannot find the source code for drop.links online to try to
> translate this code into add.links directly, as was also suggested.
>
>
>
> But such a function would I think be enormously helpful and of invaluable
> assistance for final editing.
>
>
>
> Thankyou ever so much in advance for your kind and gracious assistance.
>
>
>
> Yours sincerely,
>
>
>
> Stuart Reece.
>
>
>         [[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

Weight which work in Map Drawing but Fail in spml and spgm Regressions

Wed, 08/07/2019 - 10:07
Dear R-Sig-Geo List,

The other issue I have is weight which work beautifully in drawing maps, but
fail in regressions.

The usual error I get about this is that the subscript is out of bounds..
Which I can hardly understand...????

Thanks so much for any assistance,

Yours sincerely,

Stuart Reece.


        [[alternative HTML version deleted]]

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

Adding a Few Neighbour Relationships to a nb List

Wed, 08/07/2019 - 09:43
Dear R-Sig-Geo list,

 

I was wondering if it might be possible please to request assistance with
adding some nb relationships to a .nb.gal list composed either by GeoDa or
poly2nb in R????

 

The shapefile at this URL
<https://www.samhsa.gov/data/report/2014-2016-nsduh-substate-region-shapefil
e>  divides USA into 395 substate regions.  For health and demographical
reasons it is important to include both Hawaii and Alaska in the
spatiotemporal analysis so I want to introduce these into the Southeast
coast of California and the Pacific northwest respectively.

 

This is just as Giovanni Millo added in spatial relationships for Sicily
across the Strait of Messina for splm on page 7 of the splm pdf.

 

I found edit.nb in spdep and operated it just as described in the
instructions and here
<https://github.com/r-spatial/spdep/blob/master/man/edit.nb.Rd> .  It
crashed RStudio many times but ran well in R3.6.1.  However even though I
assigned it to a new object it did not save well.  Although when I plotted
the dxxx file as the difference between the old and modified files it
plotted the changes beautifully in red and black respectively when plotted
by themselves it introduced many long distance extraneous relationships.  To
get the edited nb list file out of R 3.6.1 and into RStudio I saved it as an
RDS file.  However when opened in RStudio it was grossly erroneous and
included extraneous links from Hawaii to Boston and New York.  When I opened
the file in RStudio it again introduced these extraneous links.

 

Saving it as a further new object in R 3.6.1 did not remedy these
difficulties.

 

The other problem I have is that the spdep poly2nb function excludes
Richmond, an island off the southern tip of Long Island near New York as it
is an island.  Also one of the areas - Region 10 in Washington DC - is also
excluded for reasons of which I am unsure.

 

I found some code here to just patch single areas
<https://stat.ethz.ch/pipermail/r-sig-geo/2006-June/001073.html>  like this
but when I run it, it throws an integer error

"  INTEGER() can only be applied to a 'integer', not a 'double'  

 

No combination of bracketing around subscripts helps or works at all.  The
link mentioned has these statements in it

 

nb[[ij[1]]] <- sort(unique(c(nb[[ij[1]]], ij[2])))
nb[[ij[2]]] <- sort(unique(c(nb[[ij[2]]], ij[1])))

 

which makes me think that I should insert a vector "  c(i,j) "  where
indicated.  Even using "c(as.integer(i),as.integer(j)) "

or " as.integer(c(I,j)) "   doesn't work and still gives rise to the same
error.

 

I am sure I am not the only one to have encountered such difficulties but I
have really tried everything I can think of.

 

The other thing I would really like is some clear instructions as to the
true underlying structure of the nb list.  If I could clearly understand
this then I could just go into the affected lines of the list of lists and
edit them directly.

 

However I am quite unable to find any clear description of its structure on
line.

 

Similarly I cannot find the source code for drop.links online to try to
translate this code into add.links directly, as was also suggested.

 

But such a function would I think be enormously helpful and of invaluable
assistance for final editing.

 

Thankyou ever so much in advance for your kind and gracious assistance.

 

Yours sincerely,

 

Stuart Reece.


        [[alternative HTML version deleted]]

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

ENFA in R: crossprod(Ze, DpZ) error

Tue, 08/06/2019 - 19:32
Dear Members,

I've like to use ENFA approach with my data set, but I have an *Error in
crossprod(Ze, DpZ) : non-conformable arguments* as output when I try to
use this function. I check my input objects and my coordinates are
in??'SpatialPointsDataFrame' and??the environmental attributes as
SpatialPixelsDataFrame" like in the help of the enfa() function. But, in
my example:

#Packages
library(dismo)
library(raster)
library(adehabitatMA)
library(adehabitatHS)

#Bradypus presence data set
file <- paste(system.file(package="dismo"), "/ex/bradypus.csv", sep="")
bradypus <- read.table(file, header=TRUE, sep="," )
# I don't need the first column
bradypus <- bradypus[,-1]
id=1:length(bradypus[,1])
bradypus.data <- as.data.frame(cbind(bradypus,id))
coordinates(bradypus.data) <- ~lon+lat
bradypus.data<-as(bradypus.data, 'SpatialPointsDataFrame')
proj4string(bradypus.data) <- CRS("+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0")


#Environmental variables
clim=getData('worldclim', var='bio', res=10)
clim.bradypus<-raster::crop(clim,c(-85.9333, -40.0667, -23.45, 13.95)) #
trim to a smaller region
plot(clim.bradypus[[1]])
plot(bradypus.data,add=T)
bradypus.map<-as(clim.bradypus, "SpatialPixelsDataFrame")
proj4string(bradypus.map) <- CRS("+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0")

## We prepare the data for the ENFA
tab <-bradypus.map
pr <- slot(count.points(bradypus.data, bradypus.map), "data")[,1]

# Run ENFA and make predictions of habitat suitability index
#First I create and removed NA (my original data set has NA)
idx = sample(1:nrow(tab),100)
tab[idx,] = NA
f1 <- function(vec) {
m <- mean(vec, na.rm = TRUE)
vec[is.na(vec)] <- m
return(vec)
}
tab2 = apply(tab,2,f1)
enfa.pink<- enfa(dudi.pca(tab2, scannf=FALSE), pr, scannf=FALSE, nf=2)##
Applied ENFA


*Error in crossprod(Ze, DpZ) : non-conformable arguments*

1: In x * pr : longer object length is not a multiple of shorter object
length ...


Any ideas about I can fix it?

Thanks in advanced

--
======================================================================
Alexandre dos Santos
Prote????o Florestal
IFMT - Instituto Federal de Educa????o, Ci??ncia e Tecnologia de Mato Grosso
Campus C??ceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
C??ceres - MT?????????????????????????????????????????? CEP: 78.200-000
Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)

[hidden email]
Lattes: http://lattes.cnpq.br/1360403201088680
OrcID: orcid.org/0000-0001-8232-6722
Researchgate: www.researchgate.net/profile/Alexandre_Santos10
LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
======================================================================





        [[alternative HTML version deleted]]

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

Re: raster::stack fails with "missing value where TRUE/FALSE needed"

Tue, 08/06/2019 - 15:36
For some reason, RStudio was linked to the R 32-bits version. I manually
defined the link to the 64-bits (in RStudio's global options) and the
problem went away.
Hugo

Hugo Costa <[hidden email]> escreveu no dia segunda, 5/08/2019 à(s)
22:49:

> Dear list,
>
> is there any reason for raster::stack retrieve and error (see below) for
> the raster provided in dropbox
> <https://www.dropbox.com/s/tdfxezlieka1ofd/gsw.tif?dl=0>?
> Thanks
>
> # this may read the raster from dropbox directly, but is slow
> # r<-raster("https://www.dropbox.com/s/tdfxezlieka1ofd/gsw.tif?dl=1")
>
> # download raster from
> https://www.dropbox.com/s/tdfxezlieka1ofd/gsw.tif?dl=0
> r<-raster("path/to/downloaded/raster")
> stack(r,r)
>
> Error in if (common.len == 1L) unlist(x, recursive = FALSE) else if
>> (common.len >  :
>>   missing value where TRUE/FALSE needed
>>
>
> # reduce extent (no error)
> r2<-crop(r, extent(r)/2)
> stack(r2,r2)
>
> R version 3.5.2 (2018-12-20)
>> Platform: i386-w64-mingw32/i386 (32-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
>> [4] LC_NUMERIC=C                            LC_TIME=English_United
>> Kingdom.1252
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> other attached packages:
>> [1] raster_2.9-23 sp_1.3-1
>>
>> loaded via a namespace (and not attached):
>> [1] compiler_3.5.2   rgdal_1.4-4      tools_3.5.2      Rcpp_1.0.0
>> codetools_0.2-15 grid_3.5.2       lattice_0.20-38
>>
>
        [[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

Tue, 08/06/2019 - 06:07
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

rgrass7 0.2-1 on CRAN

Tue, 08/06/2019 - 05:24
As use of sf and stars classes increases, the rgrass7 interface package
between R and GRASS has been changed to drop the assumption that sp (and
rgdal) will be loaded and attached. Users should now declare explicitly
with the new functions use_sp() and use_sf() which representation is
preferred for objects being moved between R and GRASS. At present,
use_sf() only supports vector objects.

The package is now on https://github.com/rsbivand/rgrass7 with
documentation on https://rsbivand.github.io/rgrass7/.

Please raise issues on https://github.com/rsbivand/rgrass7/issues, here,
or on [hidden email].

Roger

--
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

Pages