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

[Xiaohua Dai]

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?

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