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

[Roger Bivand]

On Fri, 16 Feb 2007, Debarchana Ghosh wrote:

> Hi,
> 
> I followed the steps as provided by Roger in the previous emails on a test
> ESRI shapefile with 1300 polygons. Everything went fine and I got the
> intended output. But when I tried to use the same functions on a ESRI
> shapefile with 50,000 polygons It gave the following error.
> 
> > scot<-readOGR("parcels_scot.shp", layer="parcels_scot")
> OGR data source with driver: ESRI Shapefile
> Source: "parcels_scot.shp", layer: "parcels_scot"
> with  49958  rows and  66  columns
> Error in readOGR("parcels_scot.shp", layer = "parcels_scot") :
>         unsupported type
> 
> The parcel_scot.shp has lakes and some other water features in it. If this
> is a problem, do I have to use readShapePoly() in maptools?

No, the underlying code for both of these is the same, shapelib, (but in
maptools has some extra tweaks for non-conforming shapefiles).

You can use getinfo.shape() in maptools to show what shapefile type is
being returned. The ones permitted in readShapePoly() are 1 (POINT), 8
(MULTIPOINT), 11 (POINTZ), 3 (ARC), 13 (ARCZ - 3rd dimension discarded), 5
(POLYGON), and 15 (POLYGONZ - 3rd dimension discarded). The readOGR()  
function using the current OGR shapefile driver appears to handle POLYGON
5, POLYGONM 25, and POLYGONZ 15, ARC 3, ARCM 23, and ARCZ 13, and
MULTIPATCH 31, as well as POINT 1, POINTM 21, POINTZ 11, and MULTIPOINT 8,
MULTIPOINTM 28, and MULTIPOINTZ 18.

The third dimensions are discarded from non-points.

The extra polygons for lakes and water bodies could be the problem, 
because there are no data for them, but it appears that the shape type is 
not among those supported. So the first question is what does 
getinfo.shape() say? Could the file (or a subset with the same problems) 
be made available (say on a website for a short time) for debugging? And 
how was the file created - were there options that could be turned off to 
make it simpler?

Hope this helps,

Roger

> 
> Thanks,
> Debs.
> 
> On 2/15/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> >
> > On Thu, 15 Feb 2007, Debarchana Ghosh wrote:
> >
> > > Hi,
> > >
> > > I'm extremely sorry for not giving the entire information about the
> > problem.
> >
> > OK, it makes a difference for some things.
> >
> > > I'm using current R version on Windows on a 2 GB Ram PC. The polygon is
> > in
> > > the ESRI shapefile format. I tried to do this in Arc but of course
> > failed.
> > > So can I do this in R following the steps you've mentioned in your
> > reply?
> >
> > In principle, yes. Reading the shapefile with readOGR() ought to work, and
> > will use less memory than readShapePoly() in maptools. The poly2nb()
> > function will grind on, overnight if need be, but does depend on the
> > numbers of coordinates on polygon boundaries, the more there are, the
> > longer it takes. It uses C functions internally, and was made faster
> > thanks to helpful comments by Hisaji Ono, but is not designed for speed
> > (rather to build an "nb" object once and save, reuse, or export as GAL).
> > Please keep us posted on your progress, there are other tricks that can be
> > tried if you get stuck!
> >
> > writeOGR() should also be able to get the data back out again.
> >
> > <ad>There will be a workshop on using R with spatial data at the
> > Association of American Geographers Conference in San Francisco, Tuesday
> > 17 April</ad>
> >
> > Roger
> >
> >
> > > Will there be any changes?
> > >
> > > Thanks,
> > > Debs.
> > >
> > > On 2/15/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > > >
> > > > 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
> > > >
> > > >
> > >
> > >
> > >
> >
> > --
> > 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
> >
> >
> >
> 
> 
> 

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