通过给定角度的点创建空间线

Question

我有一个空间多边形，例如格陵兰的形状：

library("rnaturalearth")
Greenland <- ne_countries(scale = "medium", country="Greenland")
plot(Greenland)

我知道如何计算它的几何中心，即质心：

library("rgeos")
Greenland_cent <- gCentroid(Greenland)
plot(Greenland_cent,add=T,col="red")

现在我想创建各种角度的线，这些线都穿过质心，以便将格陵兰岛切成两半（大约，编辑：正如@Henry 所提到的，这可能不完全是一半）。但是，我还没有找到任何允许创建由角度指定的线的方法...... ...有什么想法吗？

或者还有其他方法可以将一个多边形切割成大小差不多甚至完全一样的两部分吗？

Answer 1

我用三角函数实现了这条线。我假设你总是希望线条在格陵兰岛或格陵兰岛之外开始和结束，所以我找到了最小的封闭圆圈并使线条开始或结束于此。将一个形状切成两半需要更多的工作，但我想你可以使用坐标创建两个半圆，并使用它们来掩盖原始格陵兰形状以制作两个新的。从那里质心和区域很容易。

library(rnaturalearth)
library(sp)
library(rgeos)
library(raster)
projection <- '+proj=stere +lat_0=90 +lat_ts=90 +lon_0=-33 +k=0.994 +x_0=2000000 +y_0=2000000 +datum=WGS84 +units=m +no_defs '
Greenland.raw <- ne_countries(scale = "medium", country="Greenland")
# use a different projection to avoid some warnings
Greenland = spTransform(Greenland.raw, CRSobj = projection)
Greenland_cent <- gCentroid(Greenland)
# Make a simpler Greenland shape - expand and contract to smooth out the crinkly edges 
SimplerGreenland <- gBuffer(gBuffer(Greenland, width = 1), width=-1)
# Make a line from the shape
SimplerGreenlandLine <- as(SimplerGreenland, 'SpatialLines')
# sample all around the simpler line
PointsAroundGreenland <- spsample(x = SimplerGreenlandLine, n = 1000, type = 'regular')
# find the furthest point from the centroid
distances <- gDistance(PointsAroundGreenland, Greenland_cent, byid = T)
furthestpointfromcentroid <- distances[which.max(distances)]
# make the smallest circle that encloses Greenland to ensure that any line we draw will be in Greenland
enclosingcircle <- gBuffer(spgeom = Greenland_cent, width = furthestpointfromcentroid, quadsegs = 100)
# choose an angle and make a line that goes through the centroid and starts and ends on the enclosing circle
angleindegrees <- 5
angle <- angleindegrees * pi / 180
dx <- furthestpointfromcentroid * cos(angle)
dy <- furthestpointfromcentroid * sin(angle)
centroid <- as.numeric(coordinates(Greenland_cent))
x1 <- centroid[1] + dx
y1 <- centroid[2] + dy
x2 <- centroid[1] - dx
y2 <- centroid[2] - dy
l1 <- Line(cbind(c(x1,x2), c(y1,y2)))
l2 <- Lines(list(l1), ID = 'line')
lineatangle <- SpatialLines(list(l2))
crs(lineatangle) <- crs(Greenland)

# plot everything
plot(Greenland)
plot(Greenland_cent, add = T, col = 'red')
plot(lineatangle, add = T, col = 'green')

Answer 2

这是另一个问题的答案 - 即找到将形状分成 2 个相等部分的角度。它不使用直线计算所以我决定将其作为第二个答案以避免丢失第一个答案。

library(rnaturalearth)
library(sp)
library(rgeos)
library(raster)
library(dplyr)
library(swfscMisc)
library(ggplot2)
projection <- '+proj=stere +lat_0=90 +lat_ts=90 +lon_0=-33 +k=0.994 +x_0=2000000 +y_0=2000000 +datum=WGS84 +units=m +no_defs '
Greenland.raw <- ne_countries(scale = "medium", country="Greenland")
# use a different projection to avoid some warnings
Greenland = spTransform(Greenland.raw, CRSobj = projection)
Greenland_cent <- gCentroid(Greenland)
# Make a simpler Greenland shape - expand and contract to smooth out the crinkly edges 
SimplerGreenland <- gBuffer(gBuffer(Greenland, width = 1), width=-1)
# Make a line from the shape
SimplerGreenlandLine <- as(SimplerGreenland, 'SpatialLines')
# sample all around the simpler line
PointsAroundGreenland <- spsample(x = SimplerGreenlandLine, n = 1000, type = 'regular')
# find the furthest point from the centroid
distances <- gDistance(PointsAroundGreenland, Greenland_cent, byid = T)
furthestpointfromcentroid <- distances[which.max(distances)]
findArc <- function(angleindegrees, country_centroid, radius) {
    startandend = c(90-angleindegrees, 90-angleindegrees+180)
    centroid <- as.numeric(coordinates(country_centroid))
    arc = circle.polygon(
        x = centroid[1], 
        y = centroid[2], 
        radius = radius, units = 'km', 
        poly.type = 'cartesian', sides = 100, by.length = F,
        brng.limits = startandend)
    arc = arc[!duplicated(arc),]
    segment = centroid %>% rbind(arc) %>% rbind(centroid)
    p = Polygon(segment)
    ps = Polygons(list(p), 1)
    sps = SpatialPolygons(list(ps))
    crs(sps) = crs(country_centroid)
    sps
}
findHalf <- function(angleindegrees, country_centroid, radius, country_shape) {
    areaCountry = gArea(country_shape)
    sps <- findArc(angleindegrees, country_centroid, radius)
    halfCountry <- crop(country_shape, sps)
    areaHalfCountry <- gArea(halfCountry)
    ratio <- areaHalfCountry / areaCountry
    ratio
}

angle <- seq(from = 1, to = 70, by = 1)
lratios <- lapply(angle, function(a) findHalf(a, Greenland_cent, furthestpointfromcentroid, Greenland))
ratios.df <- do.call(rbind.data.frame, lratios)
names(ratios.df) <- 'ratio'
results.df <- data.frame(angle = angle) %>% dplyr::bind_cols(ratios.df)

ggplot(results.df, aes(x=angle, y=ratio)) + geom_line() + ggtitle('Ratio of half country to whole country')

结果图如下

10 度角和大约 69 度角将格陵兰岛分成 2 个相等的区域（假设是一个 2d 国家）。

通过给定角度的点创建空间线

Create spatial line through point at a given angle

r

polygon

spatial

angle

centroid