gCentroid(rgeos)R与实际质心(在python中) [英] gCentroid (rgeos) R vs. Actual Centroid (in python)

问题描述

摘要:我认为在R中使用gCentroid会返回一组点的质心，但是我意识到由于某种原因，它实际上返回的是几何平均值而不是质心

我想复制我在R中所做的质心计算:

I wanted to replicate a centroid calculation I did in R:

gCentroid {rgeos}

gCentroid {rgeos}

这些点的重心:

34.7573,    -86.678606  
38.30088,   -76.520266  
38.712147,  -77.158616  
39.704905,  -84.126463  

...使用R脚本...

... using the r-script ...

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中构建了以下示例-我想使用以下方法来复制计算结果:

I constructed the following example in python - I wanted to replicate the calculation, using:

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 

通过更多的挖掘，我添加了一个给我几何均值的函数:

With a bit more digging I added a function give me the geometric mean:

Centroid calculation: [ 37.72876321 -82.35113685]
Mean calculation [ 37.94873834 -81.18378815]

我意识到，由于某种原因，gCentroid似乎是在计算几何均值，而不是特征质心(我添加了一个均值函数，您可以看到它与R结果匹配)

I realised that for some reason the gCentroid seems to be calculating the geometric mean not the feature centroid (I have added a mean function, which you can see matches the R-result)

我认为原因可能是:由于我有一组点，而不是像它们一样，在它们中拟合随机多边形(例如本例)，甚至是凸包，然后取其质心，因此命令将如果数据类型为点"，则默认为均值计算.因此，我明确为其传递了一个多边形:

I thought that perhaps the reason was: since I had a grouping of points, instead of fitting a random polygon through them - like me in the example - or even a convex hull and then taking the centroid of that, the command would default to a mean calculation if the data-type was 'point'. So I explicitly passed it a polygon:

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质心在哪里:

Where the python centroid is:

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 )绘制其质心，然后以绿色( -5941009 7067939 )在同一多边形上绘制质心(使用python)以及蓝色的简单均值( -5974304 7038880 ):

And then plotted the centroid of this in red (-5875318 7010915) and then the centroid on the same polygon (using python) in green (-5941009 7067939) and the simple mean (-5974304 7038880) in blue:

推荐答案

事实证明:如果提供了一组点"，则无需通过点猜测多边形或生成凸包-该命令将自动执行给您投影坐标的平均值.

It turns out that: if a group of 'Points' are supplied, then instead of guessing a polygon through the points or producing a convex hull - the command automatically gives you the mean of the projected co-ordinates.

但是，如果提供多边形，则会得到一个质心(与python脚本相同)-在我的python示例中，我缺少一个坐标:

However, if you supply a polygon then you get a centroid (the same as the python script) - in my python example I was missing one co-ordinate:

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)

一切都很好:

我在博客如果有兴趣.

这篇关于gCentroid(rgeos)R与实际质心(在python中)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！