总结:我认为在 R 中使用 gCentroid 会返回一组点的质心，但是我意识到由于某种原因它实际上返回几何平均值而不是质心
我想复制我在 R 中所做的质心计算:
gCentroid {rgeos}
这些点的质心:
34.7573, -86.678606
38.30088, -76.520266
38.712147, -77.158616
39.704905, -84.126463
...使用 r 脚本...
require(rgdal)
require(rgeos)
no_am_eq_co <- "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
wgs84 <- "+proj=longlat +datum=WGS84"
df <- as.data.frame(list(c(34.7573,
38.30088,
38.712147,
39.704905),
c(-86.678606,
-76.520266,
-77.158616,
-84.126463)))
df$Name <- "points_A"
colnames(df) <- c("lat", "lon", "Name")
# FROM: Coordinates are geographic latitude/longitudes
coordinates(df) <- c("lon", "lat")
proj4string(df) <- CRS(wgs84)
# TO: Project into North America Equidistant Conic
df <- spTransform(df, CRS(no_am_eq_co))
# Get centroids
ctrs <- lapply(unique(df$Name),
function(x) gCentroid(SpatialPoints(df[df$Name==x,])))
ctrsout <- setNames( ctrs , unique(df$Name ) )
# Create data frame
df <- do.call(rbind, lapply(ctrsout, data.frame, stringsAsFactors=FALSE))
coordinates(df) <- c("x", "y")
proj4string(df) <- CRS(no_am_eq_co)
df <- as.data.frame(spTransform(df, CRS(wgs84)))
names(df) <- c("longitude", "latitude")
print(df$latitude)
print(df$longitude)
来到:
37.94873834, -81.18378815
我在 python 中构建了以下示例 - 我想使用以下方法复制计算:
import numpy as np
from pyproj import Proj, transform
# Using: http://www.spatialreference.org/ref/esri/102010/ we get the Proj4js format
na_eq_co = "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
wgs84 = "+proj=longlat +datum=WGS84"
def proj_arr(points,proj_from,proj_to):
inproj = Proj(proj_from)
outproj = Proj(proj_to)
func = lambda x: transform(inproj,outproj,x[0],x[1])
return np.array(list(map(func, points)))
def get_polygon_centroid(polygon):
#https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
pol = np.array(polygon)
if np.any(pol[-1] != pol[0]):
pol = np.append(pol,[pol[0]], axis=0)
pol_area = get_polygon_area(pol)
x = pol[:,0]
y = pol[:,1]
Cx = np.sum((x[:-1] + x[1:]) * ((x[:-1] * y[1:]) - (y[:-1] * x[1:]))) / (6. * pol_area)
Cy = np.sum((y[:-1] + y[1:]) * ((x[:-1] * y[1:]) - (y[:-1] * x[1:]))) / (6. * pol_area)
return np.array([Cx, Cy])
def get_polygon_area(polygon):
pol = np.array(polygon)
x = pol[:,0]
y = pol[:,1]
return np.sum( (x[:-1] * y[1:]) - (y[:-1] * x[1:]) ) / 2
def get_polygon_mean(polygon):
pol = np.array(polygon)
x = pol[:,0]
y = pol[:,1]
return np.array([np.mean(x),np.mean(y)])
def run_test(points):
points = points[:,::-1] #Flip-axis (so that longitude x-axis, latitude y-axis)
points_proj = proj_arr(points,wgs84,na_eq_co)
centroid_proj = get_polygon_centroid(points_proj)
mean_proj = get_polygon_mean(points_proj)
centroid = proj_arr([centroid_proj],na_eq_co,wgs84)
mean = proj_arr([mean_proj],na_eq_co,wgs84)
return (centroid[:,::-1][0], mean[:,::-1][0])
if __name__ == '__main__':
my_points = np.array([[34.7573,-86.678606],
[38.30088,-76.520266],
[38.712147,-77.158616],
[39.704905,-84.126463]])
test = run_test(my_points)
print("Centroid calculation: {0}\nMean calculation {1}".format(test[0],test[1]))
由此我得到:
37.72876321 -82.35113685
不是:
37.94873834,-81.18378815
经过进一步挖掘，我添加了一个函数，给出几何平均值:
Centroid calculation: [ 37.72876321 -82.35113685]
Mean calculation [ 37.94873834 -81.18378815]
我意识到，由于某种原因，gCentroid 似乎正在计算几何平均值而不是特征质心(我添加了一个均值函数，您可以看到它与 R 结果匹配)
编辑:
我认为原因可能是:因为我有一组点，而不是通过它们拟合随机多边形(就像示例中的我一样)，甚至是凸包，然后取其质心，该命令将如果数据类型为“点”，则默认为平均值计算。所以我明确地向它传递了一个多边形:
x = readWKT(paste("POLYGON((-6424797.94257892 7164920.56353916,
-5582828.69570672 6739129.64644454,
-5583459.32266293 6808624.95123077,
-5855637.16642608 7316808.01148585,
-5941009.53089084 7067939.71641507,
-6424797.94257892 7164920.56353916))"))
python_cent = readWKT(paste("POINT(-5941009.53089084 7067939.71641507)"))
r_cent = gCentroid(x)
plot(x)
plot(r_cent,add=T,col='red')
plot(python_cent, add=T,col='green')
Python 质心在哪里:
centroid = get_polygon_centroid(np.array([[-6424797.94257892, 7164920.56353916],
[-5582828.69570672, 6739129.64644454],
[-5583459.32266293, 6808624.95123077],
[-5855637.16642608, 7316808.01148585],
[-6424797.94257892, 7164920.56353916]]))
然后将其质心绘制为红色(-5875318 7010915)，然后将同一多边形(使用Python)上的质心绘制为绿色(-5941009 7067939)以及蓝色的简单平均值 (-5974304 7038880):
最佳答案
事实证明:如果提供了一组“点”，那么该命令将自动为您提供投影坐标的平均值，而不是通过这些点猜测多边形或生成凸包。
但是，如果你提供一个多边形，那么你会得到一个质心(与 python 脚本相同) - 在我的 python 示例中，我缺少一个坐标:
centroid = get_polygon_centroid(np.array([[-6424797.94257892, 7164920.56353916],
[-5582828.69570672, 6739129.64644454],
[-5583459.32266293, 6808624.95123077],
[-5855637.16642608, 7316808.01148585],
[-5941009.53089084, 7067939.71641507],
[-6424797.94257892, 7164920.56353916]]))
#polygon closed
#[-5875317.84402261 7010915.37286505]
所以运行这个 R 脚本:
x = readWKT(paste("POLYGON((-6424797.94257892 7164920.56353916,
-5582828.69570672 6739129.64644454,
-5583459.32266293 6808624.95123077,
-5855637.16642608 7316808.01148585,
-5941009.53089084 7067939.71641507,
-6424797.94257892 7164920.56353916))"))
python_cent = readWKT(paste("POINT(-5875317.84402261 7010915.37286505)"))
r_cent = gCentroid(x)
plot(x)
plot(r_cent,add=T,col='red', pch = 0)
plot(python_cent, add=T,col='green', pch = 1)
一切都很匹配:
我在我的 blog 上添加了更多信息如果有兴趣的话。
关于python - gCentroid (rgeos) R 与实际质心(Python 中)，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/35720614/