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

[Roger Bivand]

Debs:

Please diregard my previous message, the problem is that one of the 
columns in the DBF is of an unsupported type, so not the shapes, the 
problem is the data. readOGR() only supports Real, Integer, and String. 
orgInfo() on the file might help to show what is going on.

readShapePoly() uses read.dbf() from foreign, which is probably more 
robust. See if readShapePoly() works better fot the types of fields you 
have in your dbf file, or try:

library(foreign)
x <- read.dbf("parcels_scot.dbf", as.is = TRUE)

to see if it can read the file. If it can, readShapePoly() will work too. 
Are there for example any Date fields?

Roger


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