Finding the centroid of a polygon?
The formula is given here for vertices sorted by their occurance along the polygon's perimeter.
For those having difficulty understanding the sigma notation in those formulas, here is some C++ code showing how to do the computation:
#include <iostream>
struct Point2D
{
double x;
double y;
};
Point2D compute2DPolygonCentroid(const Point2D* vertices, int vertexCount)
{
Point2D centroid = {0, 0};
double signedArea = 0.0;
double x0 = 0.0; // Current vertex X
double y0 = 0.0; // Current vertex Y
double x1 = 0.0; // Next vertex X
double y1 = 0.0; // Next vertex Y
double a = 0.0; // Partial signed area
// For all vertices except last
int i=0;
for (i=0; i<vertexCount-1; ++i)
{
x0 = vertices[i].x;
y0 = vertices[i].y;
x1 = vertices[i+1].x;
y1 = vertices[i+1].y;
a = x0*y1 - x1*y0;
signedArea += a;
centroid.x += (x0 + x1)*a;
centroid.y += (y0 + y1)*a;
}
// Do last vertex separately to avoid performing an expensive
// modulus operation in each iteration.
x0 = vertices[i].x;
y0 = vertices[i].y;
x1 = vertices[0].x;
y1 = vertices[0].y;
a = x0*y1 - x1*y0;
signedArea += a;
centroid.x += (x0 + x1)*a;
centroid.y += (y0 + y1)*a;
signedArea *= 0.5;
centroid.x /= (6.0*signedArea);
centroid.y /= (6.0*signedArea);
return centroid;
}
int main()
{
Point2D polygon[] = {{0.0,0.0}, {0.0,10.0}, {10.0,10.0}, {10.0,0.0}};
size_t vertexCount = sizeof(polygon) / sizeof(polygon[0]);
Point2D centroid = compute2DPolygonCentroid(polygon, vertexCount);
std::cout << "Centroid is (" << centroid.x << ", " << centroid.y << ")\n";
}
I've only tested this for a square polygon in the upper-right x/y quadrant.
If you don't mind performing two (potentially expensive) extra modulus operations in each iteration, then you can simplify the previous compute2DPolygonCentroid function to the following:
Point2D compute2DPolygonCentroid(const Point2D* vertices, int vertexCount)
{
Point2D centroid = {0, 0};
double signedArea = 0.0;
double x0 = 0.0; // Current vertex X
double y0 = 0.0; // Current vertex Y
double x1 = 0.0; // Next vertex X
double y1 = 0.0; // Next vertex Y
double a = 0.0; // Partial signed area
// For all vertices
int i=0;
for (i=0; i<vertexCount; ++i)
{
x0 = vertices[i].x;
y0 = vertices[i].y;
x1 = vertices[(i+1) % vertexCount].x;
y1 = vertices[(i+1) % vertexCount].y;
a = x0*y1 - x1*y0;
signedArea += a;
centroid.x += (x0 + x1)*a;
centroid.y += (y0 + y1)*a;
}
signedArea *= 0.5;
centroid.x /= (6.0*signedArea);
centroid.y /= (6.0*signedArea);
return centroid;
}
How to calculate the centroid of a polygon shape file in r
It isn't clear between which centroids you want to measure, so there's a sample below of couple of different ways to find some or all distances.
The code below finds the centroid of each row, shows the distance between centroids of the first two rows, and then a small distance matrix of the first 10 rows.
The crs is the original wgs 84, distances are in meters. I didn't see a need to change the crs.
library(tidyverse)
library(sf)
x <- read_sf('shapefile_dir/') # Replace with local directory name
x <- x %>% mutate(centroids = st_centroid(st_geometry(.)))
#Distance between Goeree-Overflakkee (row 1) and Dordrecht (row 2):
st_distance(x$centroids[1], x$centroids[2])
#> Units: [m]
#> [,1]
#> [1,] 53317.83
#distance matrix of first 10 rows:
x[1:10,] %>%
st_set_geometry('centroids') %>% # Since there are 2 geometry cols, this sets the active one to 'centroids'
st_distance()
#> Units: [m]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.00 53317.826 49330.003 85844.745 76152.903 69198.288 89437.405
#> [2,] 53317.83 0.000 7355.978 38618.416 37934.810 31387.648 45356.081
#> [3,] 49330.00 7355.978 0.000 38510.970 34633.274 27580.688 44165.444
#> [4,] 85844.75 38618.416 38510.970 0.000 16873.684 19856.326 8741.436
#> [5,] 76152.90 37934.810 34633.274 16873.684 0.000 7298.181 14836.265
#> [6,] 69198.29 31387.648 27580.688 19856.326 7298.181 0.000 20611.458
#> [7,] 89437.40 45356.081 44165.444 8741.436 14836.265 20611.458 0.000
#> [8,] 84873.91 42274.715 40485.556 9634.087 9915.349 15823.963 4928.774
#> [9,] 88788.51 45795.090 44243.643 10476.055 13443.233 19690.611 2255.828
#> [10,] 94857.25 50872.407 49779.022 13205.312 19479.782 25802.514 5617.354
#> [,8] [,9] [,10]
#> [1,] 84873.911 88788.511 94857.246
#> [2,] 42274.715 45795.090 50872.407
#> [3,] 40485.556 44243.643 49779.022
#> [4,] 9634.087 10476.055 13205.312
#> [5,] 9915.349 13443.233 19479.782
#> [6,] 15823.963 19690.611 25802.514
#> [7,] 4928.774 2255.828 5617.354
#> [8,] 0.000 3922.534 9994.837
#> [9,] 3922.534 0.000 6115.383
#> [10,] 9994.837 6115.383 0.000
# A distance matrix of all pairs of points can be found by not
# subsetting the data:
# x %>% st_set_geometry('centroids') %>% st_distance()
Created on 2022-03-10 by the reprex package (v0.3.0)
Finding the centroid of a polygon
The problem, as it turns out, is that Point had unsigned members. There's also a potential secondary problem, regarding floating point rounding issues.
In this algorithm, if the wind order is the wrong way around, the computed area will be negative, and when you do something like
centroid.x += (x0 + x1) * a;
you'll get an incorrect value because a will be negative and centroid.x is unsigned.
You should store the intermediate centroid values as a floating point of some kind, so that a) you won't get these signed/unsigned issues and b) so that you won't get rounding errors on every vertex; I'm not sure if these errors would end up being large enough to cause you problems, but it seems silly to risk it when it can be so easily avoided. You should cast to unsigned int (or whatever) at the very end when you return.
How can you find the centroid of a concave irregular polygon in JavaScript?
For the centroid of the 2D surface (which is likely to be what you need),
The best is to start with a little bit of maths.
I adapted it here to your own notation :
function get_polygon_centroid(pts) {
var first = pts[0], last = pts[pts.length-1];
if (first.x != last.x || first.y != last.y) pts.push(first);
var twicearea=0,
x=0, y=0,
nPts = pts.length,
p1, p2, f;
for ( var i=0, j=nPts-1 ; i<nPts ; j=i++ ) {
p1 = pts[i]; p2 = pts[j];
f = p1.x*p2.y - p2.x*p1.y;
twicearea += f;
x += ( p1.x + p2.x ) * f;
y += ( p1.y + p2.y ) * f;
}
f = twicearea * 3;
return { x:x/f, y:y/f };
}
Center of gravity of a polygon
The center of gravity (also known as "center of mass" or "centroid" can be calculated with the following formula:
X = SUM[(Xi + Xi+1) * (Xi * Yi+1 - Xi+1 * Yi)] / 6 / A
Y = SUM[(Yi + Yi+1) * (Xi * Yi+1 - Xi+1 * Yi)] / 6 / A
Extracted from Wikipedia:
The centroid of a non-self-intersecting closed polygon defined by n vertices (x0,y0), (x1,y1), ..., (xn−1,yn−1) is the point (Cx, Cy), where
and where A is the polygon's signed area,
Example using VBasic:
' Find the polygon's centroid.
Public Sub FindCentroid(ByRef X As Single, ByRef Y As _
Single)
Dim pt As Integer
Dim second_factor As Single
Dim polygon_area As Single
' Add the first point at the end of the array.
ReDim Preserve m_Points(1 To m_NumPoints + 1)
m_Points(m_NumPoints + 1) = m_Points(1)
' Find the centroid.
X = 0
Y = 0
For pt = 1 To m_NumPoints
second_factor = _
m_Points(pt).X * m_Points(pt + 1).Y - _
m_Points(pt + 1).X * m_Points(pt).Y
X = X + (m_Points(pt).X + m_Points(pt + 1).X) * _
second_factor
Y = Y + (m_Points(pt).Y + m_Points(pt + 1).Y) * _
second_factor
Next pt
' Divide by 6 times the polygon's area.
polygon_area = PolygonArea
X = X / 6 / polygon_area
Y = Y / 6 / polygon_area
' If the values are negative, the polygon is
' oriented counterclockwise. Reverse the signs.
If X < 0 Then
X = -X
Y = -Y
End If
End Sub
For more info check this website or Wikipedia.
Hope it helps.
Regards!
How may I calculate the centroid of a multipolygon in shapely (not a polygon)
You can use MultiPolygon().centroid, it's just that you can't pass that coordinate_list directly to MultiPolygon constructor as it:
/../ takes a sequence of exterior ring and hole list tuples /../
/../ also accepts an unordered sequence of Polygon instances /../
https://shapely.readthedocs.io/en/stable/manual.html#collections-of-polygons
# Based on Multipolygon sample,
# https://shapely.readthedocs.io/en/stable/code/multipolygon.py
from matplotlib import pyplot
from shapely.geometry import Polygon, MultiPolygon
from descartes.patch import PolygonPatch
# from https://github.com/shapely/shapely/blob/main/docs/code/figures.py
from figures import BLUE, BLACK, SIZE, set_limits, plot_coords, color_isvalid
fig = pyplot.figure(1, figsize=SIZE, dpi=90)
ax = fig.add_subplot(121)
set_limits(ax, -1, 6, -1, 6)
coordinate_list = [[[1,2], [2,3], [5,5]], [[0,0], [0,1], [1,0]]]
# "constructor takes a sequence of exterior ring and hole list tuples" -
# https://shapely.readthedocs.io/en/stable/manual.html#collections-of-polygons
multi = MultiPolygon([(coordinate_list[0], []), (coordinate_list[1], [])])
# "the constructor also accepts an unordered sequence of Polygon instances"
#multi = MultiPolygon([Polygon(coordinate_list[0]),Polygon(coordinate_list[1])])
plot_coords(ax, multi.centroid, color=BLACK)
for polygon in multi.geoms:
plot_coords(ax, polygon.exterior)
patch = PolygonPatch(polygon, facecolor=BLUE, edgecolor=BLUE, alpha=0.5, zorder=2)
ax.add_patch(patch)
Related Topics
Extending the C++ Standard Library by Inheritance
Force Gcc to Notify About Undefined References in Shared Libraries
Difference Between *(Pointer + Index) and Pointer[]
Porting Clock_Gettime to Windows
How to Include Header Files in Gcc Search Path
Where Is the Lock for a Std::Atomic
Why Does C++ Allow Unnamed Function Parameters
How to Statically Link Standard Library to My C++ Program
Function Pointers Casting in C++
Why Is Including "Using Namespace" into a Header File a Bad Idea in C++
Change Floating Point Rounding Mode
Why Istream Object Can Be Used as a Bool Expression
How to Initialize a Const Field in Constructor
Differencebetween Pre-Increment and Post-Increment in the Cycle (For/While)