R для вычисления площади и центроида многоугольников

Question

Создайте функции для вычисления области и центра тяжести списка полигонов (как в формате georgia.polys). Формула для площади многоугольника

, где А представляет собой площадь многоугольника, XI является я й х-координата границы многоугольника (х [г] в R), уг это я й y-координатой границы многоугольника (y [i] в ​​R), а n - количество точек, используемых для задания границы многоугольника. Предполагается, что многоугольник находится в замкнутом виде, так что xi и yi принимают то же значение, что и xn и yn. Центроид имеет координаты (Cx, Cy), где:

Вот код, который уже создан, но им не уверен, что центр тяжести координат правильно

library(GISTools) 
data("georgia") 


polyn<-function(x){ 

    poly.df<-data.frame() 

    for(d in 1:159){ 
    poly.d<-x[[d]] 
    n<-length(poly.d[,1]) 

    i<-1 
    A.sum<-0 
    C.xsum<-0 
    C.ysum<-0 

    while(i<n){ 

     A.area<-0.5*(poly.d[i,2]*poly.d[i+1,1]-poly.d[i+1,2]*poly.d[i,1]) 
     A.sum<-A.sum+A.area 

     C.x<-(1/(6*A.sum))*(poly.d[i,2]+poly.d[i+1,2])*(poly.d[i,2]*poly.d[i+1,1]-poly.d[i+1,2]*poly.d[i,1]) 
     C.xsum<-C.xsum+C.x 

     C.y<-(1/(6*A.sum))*(poly.d[i,1]+poly.d[i+1,1])*(poly.d[i,2]*poly.d[i+1,1]-poly.d[i+1,2]*poly.d[i,1]) 
     C.ysum<-C.ysum+C.y 

     i<-i+1 
    } 

    poly.df<-rbind(poly.df, c(A.sum,C.xsum,C.ysum)) 
    colnames(poly.df) <- c("Area", "Cx", "Cy") 
    } 

    poly.df 

} 

polyn(georgia.polys) 

Это некоторый результат эта функция,

  Area   Cx   Cy 
1 1326077000 4044403.4 4855396.03 
2 891511462 -2237689.5 -2962558.41 
3 740601936 10709355.7 12996988.27 

Есть ли кто-нибудь, кто может мне помочь с кодом?

Answer 1

Площадь A.sum в C.ysum и C.xsum должна быть общая площадь, но не та область, которая зависит от итератора i. Самым простым способом является разделение после вычисления площади.

Кроме того, уравнения должны проходить по индексам 1,2,...,n+1, причем последняя вершина совпадает с первой вершиной. Поэтому вы также должны изменить свой код, чтобы он зацикливался на последнем случае в суммировании уравнений.

.... 

while(i<n+1){ 

    j <- ifelse(i+1==n+1,1,i+1) # j=i+1 and j=1 for the last iteration 

    A.area<-0.5*(poly.d[i,2]*poly.d[j,1]-poly.d[j,2]*poly.d[i,1]) 
    A.sum<-A.sum+A.area 

    C.x<-(poly.d[i,2]+poly.d[j,2])*(poly.d[i,2]*poly.d[j,1]-poly.d[j,2]*poly.d[i,1]) 
    C.xsum<-C.xsum+C.x 

    C.y<-(poly.d[i,1]+poly.d[j,1])*(poly.d[i,2]*poly.d[j,1]-poly.d[j,2]*poly.d[i,1]) 
    C.ysum<-C.ysum+C.y 

    i<-i+1 
} 

C.ysum<-C.ysum/(6*A.sum) 
C.xsum<-C.xsum/(6*A.sum) 
....