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

Sat May 13 11:07:52 CEST 2006

On Sat, 13 May 2006, Xiaohua Dai wrote:

> Thanks. How about merging the intersected hulls to make a new hull? Is
> it easy to do that? It seems that I have to decide how many hulls I
> needed at first?

Yes, deciding the cluster criteria and numbers comes first. Then there are
tricks in the gpclib and (off-CRAN) r-spatial sourceforge repository spgpc
packages for doing things to Polygons objects and lists of Polygons
objects (see CRAN -> "Task Views" -> Spatial for the link). I don't think
we have an intersection checker yet, but it wouldn't be difficult to
write. Alternately, we could cut out the intersection from both, because
we know that no points lie inside it - the polygons would still touch as
the two points on each side of the intersection (if the should be one
cluster under the chosen cluster criteria, the cluster algorithm would
have joined them).

Roger

> 
> Regards
> Xiaohua
> 
> On 5/13/06, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> > 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