[R-sig-Geo] how to code for neighbour list editing

[James Rooney]

Many thanks Roger,

I'll study that and see how I get on!

James
________________________________________
From: Roger Bivand [Roger.Bivand at nhh.no]
Sent: 15 July 2013 19:51
To: James Rooney
Cc: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] how to code for neighbour list editing

On Mon, 15 Jul 2013, James Rooney wrote:

> Hi all,
>
> Ok so I have a spatial data set with over 3000 polygons. I used
> poly2nb() to identify neighbour relationships and make an nb object. Of
> the polygons, about 30 or 40 are islands and have no neighbours
> according to polynb. Prior to doing Bayesian smoothing, I wish to assign
> neighbours to those islands (and probably weight them with a weaker
> weighting than true neighbours). Anyhow I've tried the edit tool for
> manually editing neighbours. Its basically not useful however since
> there are so many polygons its too hard to see the islands.
>
> So I'd like to write some code to do something like this by manipulating
> the nblist:

Please provide example code and data to set up the problem, say US
counties or similar:

http://www2.census.gov/geo/tiger/GENZ2010/gz_2010_us_050_00_20m.zip

library(rgdal)
usc <- readOGR(".", "gz_2010_us_050_00_20m")
library(spdep)
nb <- poly2nb(usc)

>
> 1) Identify islands

no_neighs <- which(card(nb) == 0) # returns set cardinality.

To check that you get the same result using counts of graph components:

nc <- n.comp.nb(nb)
table(nc$comp.id)
no_neighs_comps <- which(unname(table(nc$comp.id)==1))
match(no_neighs_comps, nc$comp.id)

> 2) For each island find nearest neighbour

k1nb <- knn2nb(knearneigh(coordinates(usc), k=1,
   longlat=!is.projected(usc)))

It isn't just a matter of finding the one-way neighbor, as the reverse
should also be assigned (if i neighbours j, then j neighbours i, as in
poly2nb()).

> 3) Assign a neighbour relationship between the island and its nearest
> neighbour

is.symmetric.nb(nb, force=TRUE)
nb1 <- nb

res <- k1nb[no_neighs]
nb1[no_neighs] <- res
attr(nb1, "sym") <- is.symmetric.nb(nb1, force=TRUE)
# re-assign the now incorrect symmetry attribute
nb1

nb2 <- make.sym.nb(nb1)
is.symmetric.nb(nb2, force=TRUE)
nb2

I think that some new links here have contributed to changes in the
smaller subgraphs:

table(n.comp.nb(nb2)$comp.id)

> 4) Weight it appropriately (half normal weighting for example)
>

This is harder to arrange, because of the symmetric linkages, and although
it could be done, is unlikely to affect the result. Indeed, judicious use
of set.ZeroPolicyOption(TRUE) to avoid having to set zero.policy=TRUE
might be quite acceptable if there are few no-neighbour observations.

> My problem is I don't really have any understanding of the nb object or
> how it stores relationships. I've tried summary() attributes() print()
> etc etc and can't really glean an insight into how to access or edit the
> inner structure of the nblist.

Use a smaller example to get on top of it and read the vignette in spdep
as well as the help pages and their references. An nb object is a list of
integer vectors, but uses an integer vector with a single value of 0 to
signal no neighbour.

Hope this clarifies,

Roger

>
> I was hoping someone might have dealt with this before. Or maybe there
> are functions for this I have yet to discover ?
>
> Many thanks,
> James
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no