R für die Berechnung der Fläche und des Schwerpunkts von Polygonen

Question

Erstellen Sie Funktionen zum Berechnen der Fläche und des Schwerpunkts einer Polygonliste (wie im Format von georgia.polys). Die Formel für die Fläche eines Polygons ist

wobei A die Polygonfläche, xi der i-te x-Koordinate der Polygongrenze (x [i] in R), yi der i Die y-Koordinate der Polygongrenze (y [i] in R) - und n ist die Anzahl der Punkte, die zur Angabe der Polygongrenze verwendet werden. Es wird angenommen, dass das Polygon in geschlossener Form ist, so dass xi und yi den gleichen Wert wie xn und yn annehmen. Der Schwerpunkt hat die Koordinaten (Cx, Cy) Wo:

Hier den Code ist, die bereits erstellt, aber im nicht sicher, dass der Schwerpunkt Koordinate korrekt ist

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) 

Dies ist einige Ergebnis diese Funktion,

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

Gibt es jemanden kann mir mit dem Code helfen?

Answer 1

Das Gebiet A.sum in C.ysum und C.xsum sollten die Gesamtfläche, aber nicht ein Gebiet, das i auf Ihrem Iterator abhängt. Der einfachste Weg ist die Division nach der Flächenberechnung.

Außerdem sollten die Gleichungen die Indizes 1,2,...,n+1 durchlaufen, wobei der letzte Eckpunkt derselbe wie der erste Eckpunkt ist. Daher sollten Sie auch Ihren Code so ändern, dass der letzte Fall in den Summierungen der Gleichungen durchlaufen wird.

.... 

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