[R-sig-Geo] How to find neighbor polygons and analysis of their attributes
Roger Bivand
Roger.Bivand at nhh.no
Wed Feb 21 12:38:02 CET 2007
On Tue, 20 Feb 2007, Debarchana Ghosh wrote:
> >
> > Please try writeOGR() in rgdal instead, with the appropriate
> > driver="ESRI Shapefile", probably:
> >
> > writeOGR(scot, dsn=".", layer="scot", driver="ESRI Shapefile")
>
>
>
> I tried this but were giving errors.
> > writeOGR(scot, dsn=".", layer="scot", driver="ESRI Shapefile")
> Error in writeOGR(scot, dsn = ".", layer = "scot", driver = "ESRI
> Shapefile") :
> unknown data type
>
> ## the scot file has a date field.
>
> 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.
>
>
> I've put the file in this website
> http://rguha.ath.cx/~debarchana/
OK, thanks. The 21 warnings come from the fact that 21 of the DBF columns
are just missing values (or are read by R as missing values). This is:
max(nchar(x[!is.na(x)]))
in write.dbf(), in trying to work out the field width.
These include two of the three dates, only SALE_DATE has values, but they
likely ought to be checked, the highest value is long in the future. Many
of the others also have many missing values, but I suppose they are as
they should be.
> scot$SALE_DATE <- as.character(scot$SALE_DATE)
> scot$AGPRE_ENRD <- as.character(scot$AGPRE_ENRD)
> scot$AGPRE_EXPD <- as.character(scot$AGPRE_EXPD)
> table(sapply(as(scot, "data.frame"), class))
character factor numeric
3 50 14
> table(sapply(as(scot, "data.frame"), typeof))
character double integer
3 14 50
> writeOGR(scot, dsn=".", layer="scot", driver="ESRI Shapefile")
now runs for me. So does:
> writePolyShape(scot, "scot_1")
There were 21 warnings (use warnings() to see them)
The writePolyShape() DBF was much smaller than the writeOGR() one, because
writePolyShape uses write.dbf() to try to guess the width of fields, while
writeOGR() - at the moment - just uses defaults. I read scot_1.shp into
ArcMAP 9.1 successfully, I didn't try scot.shp because of the very large
DBF - my Arc is on a laptop with 256M RAM. Your lagged values were
present.
I hope you said save(scotnbs, "scotnbs.RData", compress=TRUE), to save
repeating the generation of the neighbour lists. I feel that some of the
lists will be rather strange, because I think your plots are separated by
roads, etc., which will mean that many miss almost-contiguous neighbours
across a road. However, that is how the data are.
Best wishes,
Roger
>
> Thanks for the help.
> Debs.
>
>
> 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
> >
> >
>
>
>
--
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