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@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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
