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

Roger Bivand Roger.Bivand at nhh.no
Tue Feb 20 18:51:51 CET 2007 

On Tue, 20 Feb 2007, Debarchana Ghosh wrote:

> Hi,
> 
> The following codes worked except step6. I thought I'll write all the steps
> once agian to remind you of this particular problem.
> 
> # step1
> library(maptools)
> scot<-readShapePoly("parcels_scot.shp")
> 
> # step2
> scotnbs<-poly2nb(as(scot, "SpatialPolygons"))
> #(It kept on grindng overnight, but worked)
> 
> # step 3
> scotnbs
> Neighbour list object:
> Number of regions: 49958
> Number of nonzero links: 186210
> Percentage nonzero weights: 0.007460929
> Average number of links: 3.727331
> 189 regions with no links:
> 
> # step 4
> scotwts<-nb2listw(scotnbs, style="W", zero.policy=TRUE)
> 
> # step 5
> scot$W_EMV_TOTAL<-lag(scotwts, scot$EMV_TOTAL, zero.policy=TRUE)
> 
> # step 6
> writePolyShape(scot, fn="scot.shp")

I don't think we can get to the cause of the warnings by remote debugging. 
Please try writeOGR() in rgdal instead, with the appropriate 
driver="ESRI Shapefile", probably:

writeOGR(scot, dsn=".", layer="scot", driver="ESRI Shapefile")

We may need to convert some data frame columns to other classes, but 
readOGR() is more restrictive than writeOGR().

If you'd like to get to the bottom of the warnings, please do save(scot, 
file="scot.RData", compress=TRUE) and put it on a temporary website. I 
can't see any obvious calls to max() that might not have suitable 
arguments, so it needs running under debug. 

Can readShapePoly() read the output shapefile?

Roger

> There were 21 warnings (use warnings() to see them)
> warnings()
> Warning messages:
> 1: no non-missing arguments to max; returning -Inf
> 2: no non-missing arguments to max; returning -Inf
> 3: no non-missing arguments to max; returning -Inf
> 4: no non-missing arguments to max; returning -Inf
> 5: no non-missing arguments to max; returning -Inf
> 6: no non-missing arguments to max; returning -Inf
> 7: no non-missing arguments to max; returning -Inf
> 8: no non-missing arguments to max; returning -Inf
> 9: no non-missing arguments to max; returning -Inf
> 10: no non-missing arguments to max; returning -Inf
> 11: no non-missing arguments to max; returning -Inf
> 12: no non-missing arguments to max; returning -Inf
> 13: no non-missing arguments to max; returning -Inf
> 14: no non-missing arguments to max; returning -Inf
> 15: no non-missing arguments to max; returning -Inf
> 16: no non-missing arguments to max; returning -Inf
> 17: no non-missing arguments to max; returning -Inf
> 18: no non-missing arguments to max; returning -Inf
> 19: no non-missing arguments to max; returning -Inf
> 20: no non-missing arguments to max; returning -Inf
> 21: no non-missing arguments to max; returning -Inf
> 
> I can't understand what these warnings mean. writePolyShape() writes 3
> files, "scot.shp", "scot.shp.dbf", "scot.shx". I can view the dbf file but
> not the shapefile in ESRI arccatalog.
> 
> Thanks,
> Debs.
> 
> On 2/19/07, Debarchana Ghosh <debarchana.ghosh at gmail.com> wrote:
> >
> >
> >
> > library(foreign)
> > > x <- read.dbf("parcels_scot.dbf", as.is = TRUE)
> >
> >
> > This  worked with no problem. I also did
> > summary(x) to see whether it was able to read in all the attributes. All
> > the attributes are there.
> >
> > to see if it can read the file. If it can, readShapePoly() will work too.
> > > Are there for example any Date fields?
> >
> >
> > Yes there is one Date field.  But  I'm  trying  to see  whether
> > readShapePoly() will  work  as well.
> >
> > Thanks for all the codes.
> > Debs.
> >
> > 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
> > >
> > >
> >
> >
> > --
> > Debarchana Ghosh
> > Research Assistant,
> > Department of Geography
> > University of Minnesota.
> > PH: 8143607580
> > www.tc.umn.edu/~ghos0033/ <http://www.tc.umn.edu/%7Eghos0033/>
> >
> 
> 
> 
> 

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

More information about the R-sig-Geo mailing list