有没有办法确定 sf 中多边形的主要方向？

Question

我想知道是否有现成的工具，或者是否有人开发了一种方法来确定空间形状的人轴地理方向。一般来说，我希望能够确定一个形状是东西方向还是南北方向，但理想情况下，每个形状都会有一个角度或度数测量值。

ArcGIS 提供了“计算主角度工具”，但它是为正交形状设计的，我正在使用类似斑点或至少不是非常正交的野火周界。乍一看，Arc 工具提供了非常粗略的测量结果。

我想使用一个 sf 对象来做这件事，所以举个例子，也许可以使用 sf 包中的北卡罗来纳州数据。北卡罗来纳州 100 个县的地理定位是什么？

nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)

感谢您的帮助！

Answer 1

flightplanning-R 包具有计算最小边界矩形、方向角、高度和宽度的函数。 (https://github.com/caiohamamura/flightplanning-R)

我已经稍微调整了它，并在下面的另一个函数中使用它来返回一个带有方向角和一个 POLYGON 几何列的 sf 对象。角度是从0（东西）到180（也是东西），90是南北。

# Copied function getMinBBox()
# from https://github.com/caiohamamura/flightplanning-R/blob/master/R/utils.R
# credit there given to: Daniel Wollschlaeger <https://github.com/ramnathv>


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)
  
  ## rotating calipers algorithm using the convex hull
  H    <- grDevices::chull(x)      ## hull indices, vertices ordered clockwise
  n    <- length(H)      ## number of hull vertices
  hull <- x[H, ]        ## hull vertices
  
  ## unit basis vectors for all subspaces spanned by the hull edges
  hDir  <- diff(rbind(hull, hull[1, ])) ## hull vertices are circular
  hLens <- sqrt(rowSums(hDir^2))        ## length of basis vectors
  huDir <- diag(1/hLens) %*% hDir       ## scaled to unit length
  
  ## unit basis vectors for the orthogonal subspaces
  ## rotation by 90 deg -> y' = x, x' = -y
  ouDir <- cbind(-huDir[ , 2], huDir[ , 1])
  
  ## project hull vertices on the subspaces spanned by the hull edges, and on
  ## the subspaces spanned by their orthogonal complements - in subspace coords
  projMat <- rbind(huDir, ouDir) %*% t(hull)
  
  ## range of projections and corresponding width/height of bounding rectangle
  rangeH  <- matrix(numeric(n*2), ncol=2)  ## hull edge
  rangeO  <- matrix(numeric(n*2), ncol=2)  ## orthogonal subspace
  widths  <- numeric(n)
  heights <- numeric(n)
  
  for(i in seq(along=numeric(n))) {
    rangeH[i, ] <- range(projMat[  i, ])
    
    ## the orthogonal subspace is in the 2nd half of the matrix
    rangeO[i, ] <- range(projMat[n+i, ])
    widths[i]   <- abs(diff(rangeH[i, ]))
    heights[i]  <- abs(diff(rangeO[i, ]))
  }
  
  ## extreme projections for min-area rect in subspace coordinates
  ## hull edge leading to minimum-area
  eMin  <- which.min(widths*heights)
  hProj <- rbind(   rangeH[eMin, ], 0)
  oProj <- rbind(0, rangeO[eMin, ])
  
  ## move projections to rectangle corners
  hPts <- sweep(hProj, 1, oProj[ , 1], "+")
  oPts <- sweep(hProj, 1, oProj[ , 2], "+")
  
  ## corners in standard coordinates, rows = x,y, columns = corners
  ## in combined (4x2)-matrix: reverse point order to be usable in polygon()
  ## basis formed by hull edge and orthogonal subspace
  basis <- cbind(huDir[eMin, ], ouDir[eMin, ])
  hCorn <- basis %*% hPts
  oCorn <- basis %*% oPts
  pts   <- t(cbind(hCorn, oCorn[ , c(2, 1)]))
  
  ## angle of longer edge pointing up
  dPts <- diff(pts)
  e    <- dPts[which.max(rowSums(dPts^2)), ] ## one of the longer edges
  eUp  <- e * sign(e[2])       ## rotate upwards 180 deg if necessary
  deg  <- atan2(eUp[2], eUp[1])*180 / pi     ## angle in degrees
  
  return(list(pts=pts, width=heights[eMin], height=widths[eMin], angle=deg))
}

##############
## Use getMinBBox in a custom function to return an sf object
##############
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
}

# Testing on a county in the nc dataset with an unusual shape and orientation:

min_box_sf(nc[56,])
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -76.19819 ymin: 35.11926 xmax: -75.31058 ymax: 36.23016
#> Geodetic CRS:  NAD27
#>   id                       geometry    angle
#> 1  1 POLYGON ((-76.19819 36.0092... 117.4866

#Plotting county 56 & the associated minimum bounding box
ggplot() + 
  geom_sf(data = nc[56,], 
          fill = 'red', 
          alpha = .2) + 
  geom_sf(data = min_box_sf(nc[56,]), 
          fill = NA) 

形状不寻常的北卡罗来纳州戴尔县有一个“长”方向的最小边界框，其“长”方向约为 117 度，即北-北-西到南-南-东。

# Using the function on each row of an sf object.
#  note the crs is not retained.
pmap_dfr(nc, min_box_sf)
#> Simple feature collection with 100 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.86573 xmax: -75.31058 ymax: 36.87134
#> CRS:           NA
#> First 10 features:
#>    id       angle                       geometry
#> 1   1 177.0408464 POLYGON ((-81.74847 36.2486...
#> 2   1 179.0078231 POLYGON ((-81.3505 36.36728...
#> 3   1 178.4492784 POLYGON ((-80.97202 36.2365...
#> 4   1 136.8896308 POLYGON ((-75.59489 36.2906...
#> 5   1 149.5889916 POLYGON ((-77.71197 36.8713...
#> 6   1 179.5157854 POLYGON ((-77.21774 36.2322...
#> 7   1 147.1227419 POLYGON ((-75.90195 36.2792...
#> 8   1   0.1751954 POLYGON ((-76.95329 36.2937...
#> 9   1   0.1759289 POLYGON ((-78.32017 36.1949...
#> 10  1 179.0809855 POLYGON ((-80.02092 36.5467...

将所有县的最小边界框绘制在一起：

pmap_dfr(nc, min_box_sf) %>% 
  ggplot() +
  geom_sf(alpha = .2)

^{由reprex package (v2.0.1) 于 2021-08-20 创建}