[R-sig-Geo] Extracting regions with no neighbors from Map2poly object

Fri Jun 24 15:23:19 CEST 2005

On Fri, 24 Jun 2005 Alexander.Herr at csiro.au wrote:

> Thank Roger,
> 
> I can replicate most of you example using my Map2poly object (ie
> plotting and subsetting works on the centroids and the poly2nb object).
> However when subsetting the Map2poly object I get an "Error in
> order(na.last, decreasing,...): argument 1 is not a vector.

Not on the centroids, they were used in that case only to generate an "nb" 
object (by the way your Map2poly object is a "polylist" object, Map2poly() 
is just one function which returns such an object; "nb" objects can be 
returned by many more functions than poly2nb()).

Please try debug(subset.polylist) and re-run to see where this is 
happening, most likely in:

    area <- sapply(res, function(x) attr(x, "area"))
    if (any(sapply(area, is.null))) 
        pO <- order(subset.default(attr(x, "plotOrder"), subset))
    else pO <- order(area, decreasing = TRUE)

It looks as though any(sapply(area, is.null)) is TRUE as well as the 
return value from attr(x, "plotOrder") is not a vector. This suggests that 
your "polylist" object is pretty old (?) and/or that the maptools package 
is not very fresh(?) - which version of maptools are you using, and when 
was the "polylist" object generated? What are the components of 
attributes(<"polylist">)? Is there an area attribute and/or a plotOrder 
attribute?

A warning, I'll be offline from Sunday (GMT) for some weeks, so if 
problems continue, please join in, other maptools users!

Roger

> 
> The shapefile I am using is messy and has multiple polygons in one
> region (ie islands off the coast assocaited with coastal polygons). Is
> there a way of matching region.id in the Map2poly object with the
> isolated/non-isolated polygons from the neigborlist so I can
> delete/retain polygons?

PS. I don't understand. A shape can have multiple polygons, and in 
maptools, you can only subset shapes easily, removing one or more 
sub-polygons from a shape can be done, but you have to manipulate the 
object directly. This is one of the reasons why the sp package has been 
developed, and it is where development will continue.

> 
> Thanks Herry
> 
> --------------------------------------------
> Alexander Herr - Herry
> Northern Futures
> Davies Laboratory, CSIRO
> PMB, Aitkenvale, QLD 4814
> Phone (07) 4753 8510
> Fax   (07) 4753 8650
> Home: http://herry.ausbats.org.au
> Webadmin ABS: http://ausbats.org.au
> Sustainable Ecosystems: http://www.cse.csiro.au/
> --------------------------------------------
>   
> 
> 
> 
> -----Original Message-----
> From: Roger Bivand [mailto:Roger.Bivand at nhh.no] 
> Sent: Thursday, 23 June 2005 7:44 PM
> To: Herr, Alexander Herr - Herry (CSE, Townsville)
> Cc: r-sig-geo at stat.math.ethz.ch
> Subject: Re: [R-sig-Geo] Extracting regions with no neighbors from Map2poly object
> 
> 
> On Thu, 23 Jun 2005 Alexander.Herr at csiro.au wrote:
> 
> > Hi List,
> > 
> > I am trying to delete polygons with no neighbors in a Map2poly object.
> > 
> > Region.id exists, however, I cannot find a way to identify/delete the 
> > list of polygons (from poly2nb()) with no neighbors. I have trouble 
> > subsetting the Map2poly object, so I can create a new object without 
> > these polygons.
> > 
> 
> I'll try to replicate:
> 
> > library(spdep)
> > try3 <- read.shape(system.file("shapes/sids.shp", 
> > package="maptools")[1]) mappolys <- Map2poly(try3, 
> > as.character(try3$att.data$FIPSNO))
> > cents <- do.call("rbind", lapply(lapply(mappolys, function(x) attr(x, 
> +   "centroid")), unlist))
> > str(cents)
>  num [1:100, 1:2] -81.5 -81.1 -80.7 -76.1 -77.4 ...
>  - attr(*, "dimnames")=List of 2
>   ..$ : NULL
>   ..$ : chr [1:2] "x" "y"
> > nb <- dnearneigh(cents, 0, 0.3)
> > nb
> Neighbour list object:
> Number of regions: 100 
> Number of nonzero links: 78 
> Percentage nonzero weights: 0.78 
> Average number of links: 0.78 
> 44 regions with no links:
> 6 8 10 11 12 14 24 33 37 38 39 42 44 46 51 54 55 56 57 58 62 63 64 66 67 
> 68 70 71 75 78 79 80 82 84 85 86 88 89 92 94 95 96 98 100
> > table(card(nb))
> 
>  0  1  2  3 
> 44 36 18  2 
> > plot(mappolys)
> > plot(nb, cents, add=TRUE, lwd=3, col="blue")
> 
> This generates a neighbour list with lots of no-neighbour counties, so 
> subsetting is just a matter of dropping the counties with zero neighbour 
> cardinality (subset= are the ones to retain):
> 
> > sub_mappolys <- subset(mappolys, subset=card(nb) > 0)
> > sub_nb <- subset(nb, subset=card(nb) > 0)
> > sub_cents <- subset(cents, subset=card(nb) > 0)
> > sub_nb
> Neighbour list object:
> Number of regions: 56 
> Number of nonzero links: 78 
> Percentage nonzero weights: 2.487245 
> Average number of links: 1.392857 
> > plot(mappolys, col="lightgreen")
> > plot(sub_mappolys, col="white", add=TRUE)
> > plot(sub_nb, sub_cents, add=TRUE, lwd=3, col="blue")
> 
> Of course, this is just for this data set, and maybe there is something 
> different about yours.
> 
> Hope this helps,
> 
> Roger
> 
> 
> > Any help appreciated
> > 
> > Thanks HErry
> > 
> > --------------------------------------------
> > Alexander Herr - Herry
> > Northern Futures
> > Davies Laboratory, CSIRO
> > PMB, Aitkenvale, QLD 4814
> > Phone (07) 4753 8510
> > Fax   (07) 4753 8650
> > Home: http://herry.ausbats.org.au
> > Webadmin ABS: http://ausbats.org.au
> > Sustainable Ecosystems: http://www.cse.csiro.au/
> > --------------------------------------------
> >  
> > 
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch 
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > 
> 
> 

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