[R-sig-Geo] [R] [ESRI-L] outline polygons of point clumps

[Roger Bivand]

On Fri, 12 May 2006, Xiaohua Dai wrote:

> My idea is to use your procedure to obtain food-tree patches of
> rhinos, then relate the patches to rhino movement pattern. For
> example, rhinos may move more slowly in food patches than outside the
> patches. I forgot one thing, could you also tell me how to write these
> hulls into a polygon shape file? write.polylistShape {maptools} or
> convert.to.shapefile {shapefiles}?

set.seed(1)
xy <- matrix(runif(500, 0, 10), ncol=2)
xy_clusts <- hclust(dist(xy), method="complete")
cl_10 <- cutree(xy_clusts, 10)
which_cl_10 <- tapply(1:nrow(xy), cl_10, function(i) xy[i,])
chulls_cl_10 <- lapply(which_cl_10, function(x) x[chull(x),])
# so far as before - now:

n <- length(chulls_cl_10)
library(sp)
polygons <- lapply(1:n, function(i) {
  chulls_cl_10[[i]] <- rbind(chulls_cl_10[[i]], chulls_cl_10[[i]][1,])
# the convex hulls do not join first and last points, so we copy here
  Polygons(list(Polygon(coords=chulls_cl_10[[i]])), ID=i) })
out <- SpatialPolygonsDataFrame(SpatialPolygons(polygons),
  data=data.frame(ID=1:n))
plot(out)
# note standard-violating intersecting polygons!
tempfile <- tempfile()
library(maptools)
writePolyShape(out, tempfile)
in_again <- readShapePoly(tempfile)
plot(in_again, border="blue", add=TRUE)

The trick is to construct a SpatialPolygons object defined in the sp 
package, and write it and a single data field out to a shapefile. The 
shapefile is strictly wrong, because the polygons can intersect, but at 
least it's something to start with.

Roger


> 
> Thanks
> Xiaohua
> 
> On 5/12/06, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > On Fri, 12 May 2006, Xiaohua Dai wrote:
> >
> > > Hi Roger,
> > >
> > > Thank you for your kindly help. It is exactly what I want. I will try
> > > to understand your programming.  I think the program should be useful
> > > in many fields as you indicated.
> > >
> > > I did not expect the problem can be solved so easily. Again, R is very
> > > powerful with so many good functions. I am quite interested in lapply
> > > and tapply. I will find some references/websites on their
> > > applications.
> >
> > Note that using a different method= argument to hclust, or a different
> > metric to dist(), you will get different clusters, as well as getting
> > different numbers for different cutree() values. If you have specific
> > range distances (from an off-list message, the points are rhino), you
> > could also use other methods to create your own clusters where the first
> > nearest neighbour distance is less than your empirically observed range -
> > there should be a literature on this among ecologists, which others on the
> > list may be able to help with. Is there anything in the adehabitat
> > package, or others?
> >
> > Roger
> >
> > >
> > > Also thanks to remind me of the post guide.
> > >
> > > Best wishes
> > > Xiaohua
> > >
> > > On 5/12/06, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > > > On Fri, 12 May 2006, Xiaohua Dai wrote:
> > > >
> > > > > Sorry, I did not make my question clear. Since I have a point theme
> > > > > with many points, some of them may clump together. the problems here
> > > > > are:
> > > >
> > > > A more specific follow-up after my reply directly to R-help.
> > > >
> > > > I'm assuming you have a running R, and that convex hulls are what you
> > > > want.
> > > >
> > > > Take a random data set:
> > > >
> > > > set.seed(1)
> > > > xy <- matrix(runif(500, 0, 10), ncol=2)
> > > > xy_clusts <- hclust(dist(xy), method="complete")
> > > > # complete linkage hierarchical clustering
> > > > plot(xy_clusts)
> > > > # shows the clustering tree
> > > > cl_10 <- cutree(xy_clusts, 10)
> > > > cl_20 <- cutree(xy_clusts, 20)
> > > > cl_30 <- cutree(xy_clusts, 30)
> > > > # cut the tree - the objects contain the memberships
> > > > which_cl_10 <- tapply(1:nrow(xy), cl_10, function(i) xy[i,])
> > > > chulls_cl_10 <- lapply(which_cl_10, function(x) x[chull(x),])
> > > > # construct convex hull polygons for each cluster
> > > > plot(xy)
> > > > res <- lapply(chulls_cl_10, polygon)
> > > > # and repeat for cl_20 and cl_30
> > > > which_cl_20 <- tapply(1:nrow(xy), cl_20, function(i) xy[i,])
> > > > chulls_cl_20 <- lapply(which_cl_20, function(x) x[chull(x),])
> > > > plot(xy)
> > > > res <- lapply(chulls_cl_20, polygon)
> > > > which_cl_30 <- tapply(1:nrow(xy), cl_30, function(i) xy[i,])
> > > > chulls_cl_30 <- lapply(which_cl_30, function(x) x[chull(x),])
> > > > plot(xy)
> > > > res <- lapply(chulls_cl_30, polygon)
> > > >
> > > > If you need the list of convex hulls out as a shapefile, we can do that,
> > > > if you need a raster, I'd suggest using kernel density with different
> > > > bandwidths for a start, and NA out the densities below a chosen threshold.
> > > >
> > > > Hope this helps,
> > > >
> > > > Roger Bivand
> > > >
> > > > PS. Use package maptools, function readShapePoints() to read your
> > > > shapefile, and coordinates() of the input object to extract the
> > > > coordinates. If need be, project from geographical coordinates using
> > > > transform methods in package rgdal.
> > > >
> > > >
> > > > > 1.  how to find clumps in a point theme?
> > > > > 2.  the convex-hull extension I found only deal with all the points in
> > > > > a theme at each time?  how to make each convex hull around each point
> > > > > clump automatically?
> > > > >
> > > > > Thanks.
> > > > >
> > > > > Xiaohua
> > > > >
> > > > >
> > > > >
> > > > > On 5/12/06, Bob Booth <bbooth at esri.com> wrote:
> > > > > > Xiaohua,
> > > > > >
> > > > > > That would be one way to do it. There are others.
> > > > > > Try searching ArcScripts for "convex hull"
> > > > > > http://arcscripts.esri.com/
> > > > > >
> > > > > > For example:
> > > > > > http://arcscripts.esri.com/details.asp?dbid=14535
> > > > > > or
> > > > > > http://arcscripts.esri.com/details.asp?dbid=12084
> > > > > >
> > > > > > Bob
> > > > > >
> > > > > >
> > > > > > -----Original Message-----
> > > > > > From: ESRI-L [mailto:ESRI-L at esri.com] On Behalf Of Xiaohua Dai
> > > > > > Sent: Friday, May 12, 2006 2:33 AM
> > > > > > To: ESRI-L at esri.com
> > > > > > Subject: [ESRI-L] outline polygons of point clumps
> > > > > >
> > > > > > Dear all,
> > > > > >
> > > > > > How to generate one outline polygon for each point clump? Are there
> > > > > > any present functions in ArcView, ArcGIS, R or some freewares? I just
> > > > > > had a quick look at the package adehabitat and did not find the
> > > > > > function.
> > > > > >
> > > > > > To my knowledge, I could do it as follows: 1) make a grid map of my
> > > > > > study area with cell values = 0; 2) assign 1 to the cells containing
> > > > > > at least one point; 3) convert the 1-value cells into polygons.
> > > > > >
> > > > > > Am I right?
> > > > > >
> > > > > > Thanks
> > > > > > Xiaohua
> > > > > >
> > > > > >
> > > > >
> > > > >
> > > > >
> > > >
> > > > --
> > > > 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