2

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here enter image description here

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?

4
  • So it has created a polygon, just not the one you were expecting.
    – nmtoken
    Dec 26, 2020 at 13:19
  • Yes you are right, this is not a bad start I guess!
    – Alexis
    Dec 26, 2020 at 13:53
  • Do the points have any information about their ordering round the "circle"? Or are they purely the X and Y coords? If there's no info, then you'll have to infer it by joining nearest points, or finding the centre and joining on increasing angle...
    – Spacedman
    Dec 26, 2020 at 14:02
  • Another idea might be to "buffer" the points (see st_buffer) and merge the buffers until you have a polygon of two rings, then take the average perimeter of the two rings. Your second image looks like you have a thick yellow polygon make from two rings...
    – Spacedman
    Dec 26, 2020 at 14:06

3 Answers 3

4

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:

enter image description here

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)

enter image description here

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
1
  • Thank you so much for all these tools and functions you just shared! I was not aware of their existence but they will help me a lot! Thank you again and have a nice week!
    – Alexis
    Dec 27, 2020 at 20:59
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)

enter image description here

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.

1

Here's some sample points I created:

> plot(pts)

enter image description here

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)

enter image description here

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.