da <- data.frame(Tree=rep(c("CB6", "CB7", "MB10"), each=8), ID=rep(1:8, 3),
AZ=c(240.7, 276.6, 301.1, 51.9, 110.8, 130, NA, NA, 211.3, 182.1,
178, 111.6, 29.4, 13.8, 307.7, 274.2, 325.5, 228.8, 192.9, 138.9,
92.4, 35.9, 330.8, NA),
DH=c(3.6, 5.6, 3.3, 3.3, 2.5, 6.4, NA, NA, 4.5, 4.3, 3.4, 4.6, 2.4,
4.5, 7.6, 3.2, 6.1, 10, 3.2, 3.9, 6.9, 4.2, 3.4, NA))
### remove NA's e calcula coordenadas x e y
da <- da[complete.cases(da),]
da$x = da$DH*cos(da$AZ*pi/180)
da$y = da$DH*sin(da$AZ*pi/180)
da
### "fecha" polígonos
da2 <- NULL
for (tree in unique(da$Tree)) {
sel <- da[da$Tree==tree,]
sel <- rbind(sel, sel[1,])
da2 <- rbind(da2, sel)
}
### visualiza formas das copas
par(mfrow=c(2,2))
for (tree in unique(da$Tree)) {
sel <- da[da$Tree==tree,]
plot(sel$y~sel$x, type="n", asp=T, xlab=NA, ylab=NA, xlim=c(-10,10), ylim=c(-10,10))
polygon(sel$y~sel$x, col="light green", bor=3)
area <- paste(round(splancs::areapl(cbind(sel$x, sel$y)),2 ), "m²")
text(6,6, area)
}
### Exemplos de operações com um polígono
p1 <- as.matrix(da2[1:7, 5:6])
splancs::areapl(p1) # [1] 31.77187 m²
p1.d <- as.matrix(dist(p1))
p1.l <- sapply(1:(nrow(p1.d)-1), function(x) p1.d[x, x+1])
sum(p1.l) # perimeter - 27.21928 m
require(sp)
p1.sp <- SpatialPolygons(list(Polygons(list(Polygon(p1)), "pol1")))
rgeos::gArea(p1.sp) # [1] 31.77187 m²
require(geosphere)
areaPolygon(p1.sp) # 390888183993 ???
perimeter(p1.sp) # 3014297 ???
### geosphere considera que seriam coordenadas em graus
### é possível calcular um fator de conversão
conv <- sqrt(areaPolygon(p1.sp)/rgeos::gArea(p1.sp)); conv # ~ 110919m
areaPolygon(p1.sp)/conv^2 # 31.77187 m²
perimeter(p1.sp)/conv # 27.17573 m
### Cálculo para os três polígonos
for (tree in unique(da2$Tree)) {
sel <- da2[da2$Tree==tree, 5:6]
sel.sp <- SpatialPolygons(list(Polygons(list(Polygon(sel)), "pol1")))
conv <- sqrt(areaPolygon(sel.sp)/rgeos::gArea(sel.sp))
print(paste(tree, round(areaPolygon(sel.sp)/conv^2, 1),
round(perimeter(sel.sp)/conv, 1)))
}
# Área e perímetros
# [1] "CB6 31.8 27.2"
# [1] "CB7 48.1 31.3"
# [1] "MB10 74 40.8"
### </code>