[R-sig-Geo] How to find neighbor polygons and analysis of their attributes

[Roger Bivand]

On Thu, 15 Feb 2007, Debarchana Ghosh wrote:

> Hi All,
> 
> I am trying to select a polygon (parcel), then select all its neighboring
> polygons that share a border with the original polygon so that all the
> neighbor polygons and the original one are selected. Then I want to
> calculate the mean value of one of the attributes from all the selected
> polys and assign that value to a new filed in the original poly. I would
> like to loop this so that it will do it for all the polygons in my layer. I
> have approx 50,000 polygons.

What we don't know is how your 50k polygons are stored, nor which platform 
and version of R you are using.

Assuming current R on a modern OS[1], and reading from a file format read
by OGR, then:

# step 1
library(rgdal)
my_50k_polys <- readOGR(dsn=".", layer="my_layer")
# see ?readOGR for format-dependent variants on dsn= and layer=
# my_50k_polys is an "SpatialPolygonsDataFrame" object

# step 2
library(spdep)
nbs <- poly2nb(as(my_50k_polys, "SpatialPolygons"))
# nbs is an object of class "nb"
# the need to coerce using as() will go away in the next release of spdep;
# this step will take a lunch break but the output object can be saved
# and reused

# step 3
wts <- nb2listw(nbs, style="W")
# wts is an object of class "listw"
# convert neighbours to spatial weights

#step 4
my_50k_polys$W_var <- lag(wts, my_50k_polys$var)
# to assign to the original polygons

You'll need to see whether the neighbour object has no-neighbour polygons, 
which will be reported if you just say:

nbs

at the R prompt (this invokes the print method for this class of object). 
If there are no-neighbour polygons, you can choose to add zero.policy=TRUE 
to steps 3 and 4.

Hope this helps,

Roger

[1] Windows, including XP, allocates memory differently from Unix/Linux or
OSX, so if each of the 50k polygons has thousands of coordinates in its
border, there may be problems with timings and memory allocation, and
these problems will be harder to handle on Windows. If the underlying data
are "odd" in some sense, like having rivers between polygons that should
be made neighbours, more details will be needed (see the snap= argument to
poly2nb()).

> 
> Any help with will be greatly appreciated.
> 
> Thanks,
> Debs.
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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