How to find neighbor polygons and analysis of their attributes

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

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