flightplanning-R
包有一个函数可以计算最小边界矩形、方向角度、高度和宽度。(https://github.com/caiohamamura/flightplanning-R)
我稍作调整并在另一个函数中使用它,返回一个带有方向角度和POLYGON几何列的sf对象。方向角度从0(东西)到180(也是东西),90为南北。
library(tidyverse)
library(sf)
library(sfheaders)
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE) %>%
st_geometry() %>% st_as_sf()
getMinBBox <- function(x) {
stopifnot(is.matrix(x), is.numeric(x), nrow(x) >= 2, ncol(x) == 2)
H <- grDevices::chull(x)
n <- length(H)
hull <- x[H, ]
hDir <- diff(rbind(hull, hull[1, ]))
hLens <- sqrt(rowSums(hDir^2))
huDir <- diag(1/hLens) %*% hDir
ouDir <- cbind(-huDir[ , 2], huDir[ , 1])
projMat <- rbind(huDir, ouDir) %*% t(hull)
rangeH <- matrix(numeric(n*2), ncol=2)
rangeO <- matrix(numeric(n*2), ncol=2)
widths <- numeric(n)
heights <- numeric(n)
for(i in seq(along=numeric(n))) {
rangeH[i, ] <- range(projMat[ i, ])
rangeO[i, ] <- range(projMat[n+i, ])
widths[i] <- abs(diff(rangeH[i, ]))
heights[i] <- abs(diff(rangeO[i, ]))
}
eMin <- which.min(widths*heights)
hProj <- rbind( rangeH[eMin, ], 0)
oProj <- rbind(0, rangeO[eMin, ])
hPts <- sweep(hProj, 1, oProj[ , 1], "+")
oPts <- sweep(hProj, 1, oProj[ , 2], "+")
basis <- cbind(huDir[eMin, ], ouDir[eMin, ])
hCorn <- basis %*% hPts
oCorn <- basis %*% oPts
pts <- t(cbind(hCorn, oCorn[ , c(2, 1)]))
dPts <- diff(pts)
e <- dPts[which.max(rowSums(dPts^2)), ]
eUp <- e * sign(e[2])
deg <- atan2(eUp[2], eUp[1])*180 / pi
return(list(pts=pts, width=heights[eMin], height=widths[eMin], angle=deg))
}
min_box_sf <- function(x){
crs <- st_crs(x)
x_as_matrix <- st_coordinates(x)[,1:2]
min_box <- getMinBBox(x_as_matrix)
box <- sfheaders::sf_polygon(min_box$pts) %>%
st_set_crs(crs)
box$angle <- min_box$angle
box
}
min_box_sf(nc[56,])
ggplot() +
geom_sf(data = nc[56,],
fill = 'red',
alpha = .2) +
geom_sf(data = min_box_sf(nc[56,]),
fill = NA)
不寻常的 Dare 县,NC 具有最小外接矩形,其“长”方向约为 117 度,即北北西至南南东。
pmap_dfr(nc, min_box_sf)
将所有县的最小边界框绘制在一起:
pmap_dfr(nc, min_box_sf) %>%
ggplot() +
geom_sf(alpha = .2)
![这里输入图片描述](https://istack.dev59.com/rcI0o.webp)
此内容由 reprex package (v2.0.1) 于2021-08-20创建。