Convert points layer into a circular polygon on R

Question

As the title says, I am supposed to convert a point layer similar to this one (it represents the perimeter of a tree at a certain height, and it's a sf that I called "points" in R):

I would like to convert this layer into a circular polygon in order to calculate its circumference just like in this example:

I tried to use the "st_cast" function in my script:

Unfortunately, this script doesn't generated what I expected... Indeed, the polygon generated looks like a weird shape:

I don't really know what to do to correct this dysfunction because I am not very familiar with R. Do you have any idea about what I am supposed to do in order to create a shape just like in the example?

Answer 1

Here's another method which orders the points by angle to get a first approximation to the "circumference" polygon. Use this function to get a polygon from the points, connecting points in increasing angle from the centroid (mean location) of the points:

centrit <- function(pts){
    centre = st_coordinates(st_centroid(st_union(pts)))
    pts = st_coordinates(pts)

    theta = atan2(pts[,2]-centre[,2], pts[,1]-centre[,1])
    poly = cbind(
        pts[order(theta),1],
        pts[order(theta),2])
    poly = st_polygon(list(rbind(poly, poly[1,])))
    poly
    
}

That gets you this from my test data:

and you can see how it joins the dots. The length of this line is:

> pp = centrit(pts)
> st_length(st_cast(pp,"LINESTRING"))
[1] 2.71652

which is a bit longer than the other method because it is not very smooth.

Smoothing it with st_simplify gives something smoother, but requires you to get a parameter in the sweet spot:

pps = st_simplify(pp,dTolerance=.02)

With a length now:

> st_length(st_cast(pps,"LINESTRING"))
[1] 2.158476

which is smaller than the smallest buffer in the other answer but perhaps is over-smoothed. Experimenting with the smoothing parameter can give answers within the buffer-based answer:

> pps = st_simplify(pp,dTolerance=.01)
> st_length(st_cast(pps,"LINESTRING"))
[1] 2.442864

Answer 2

Another method - fit a smoothed line based on polar coordinates of the points from the centroid.

Ingredients: a function to convert x,y to r,theta from the centroid of a set of points:

polar <- function(pts,centre){
    pts = st_coordinates(pts)
    theta = atan2(pts[,2]-centre[,2], pts[,1]-centre[,1])
    r = sqrt((pts[,1]-centre[,1])^2 + (pts[,2]-centre[,2])^2)
    data.frame(theta=theta, r=r)
}

A function to take the points, call the polar conversion, then fit a smoothed line to the polar coordinates, then transform and shift back to cartesian coordinates and make a polygon:

fitline <- function(pts){
    centre = st_coordinates(st_centroid(st_union(pts)))
    rt = polar(pts,centre)
    rtsmooth  = supsmu(rt$theta, rt$r)
    xysmooth = cbind(centre[,1]+rtsmooth$y*cos(rtsmooth$x),
                     centre[,2]+rtsmooth$y*sin(rtsmooth$x))
    xysmooth = rbind(xysmooth, xysmooth[1,])
    st_polygon(list(xysmooth))
}

To test on my sample data:

> tree = fitline(pts)
> plot(tree)
> plot(pts, add=TRUE)

This line has a length of:

> st_length(st_cast(tree,"LINESTRING"))
[1] 1.958053

The only tuning this method needs is to control the degree of smoothing. Lots of other smoothing functions could be used, I've used supsmu since its in the base R packages and is pretty much guaranteed to do something reasonable.

Ideally you'd fit a periodic function so that your line joins up at the 360-0 degree gap but with dense enough data this probably isn't a problem.

Answer 3

Here's some sample points I created:

> plot(pts)

Now using a carefully selected buffer size wb (not too big, not too small) I can create a single thinnish polygon that contains the points:

> wb = 0.015
> buf = st_buffer(st_union(pts), wb)
> plot(buf,add=TRUE)

Experiment with wb sizes to see the effect.

That polygon is an outer ring and an inner ring, which you can get by casting the polygon to LINESTRING:

> rings = st_cast(buf, "LINESTRING")

Then you can get the length of the rings:

> st_length(rings)
[1] 2.535989 2.195873

and then a plausible estimate for your "circumference" is:

> mean(st_length(rings))
[1] 2.365931

Note this doesn't actually compute the perimeter line, it just has to have a length between the inner and outer one. Again experiment with wb to see what happens if you have it smaller or bigger. If you have a lot of these to do with different size trees you'll have to think of a way to automate choice of wb. Hard to give advice on that without knowing exactly what your point-generating process is, but if you choose it such that you get exactly two rings with similar lengths then that might be a good start, but it might get in trouble if you have any really outlying points.