r - 在多边形中生成规则间隔的点

有没有办法使用 R 在多边形内生成规则间隔(例如，相距 500 米)的点？我一直在尝试使用 sp 包，但似乎无法定义一组彼此间隔一定距离的点。我的目标是生成点，然后将它们的经纬度坐标提取到新的数据帧中。任何帮助将非常感激!谢谢

最佳答案

非常直接，几乎是开箱即用的。

由于OP没有共享数据，请系好安全带，将您的座位垂直放置，让我们飞往巴黎。在那里，我们将调整一个 geosphere 函数，在它的帮助下，我们将巴黎的形状划分为 lon/lat 坐标，每个坐标相距 500 米(垂直和水平)。

# Load necessary libraries.
library(raster)
library(geosphere)
library(tidyverse)
library(sp)

# This is an adapted version of geosphere's destPoint() function that works with
# changing d (distance).
destPoint_v <- function (x, y, b, d, a = 6378137, f = 1/298.257223563, ...)
{
    r <- list(...)$r
    if (!is.null(r)) {
        return(.old_destPoint(x, y, b, d, r = r))
    }
    b <- as.vector(b)
    d <- as.vector(d)
    x <- as.vector(x)
    y <- as.vector(y)
    p <- cbind(x, y, b, d)
    r <- .Call("_geodesic", as.double(p[, 1]), as.double(p[, 2]),
               as.double(p[, 3]), as.double(p[, 4]),
               as.double(a), as.double(f),
               PACKAGE = "geosphere")
    r <- matrix(r, ncol = 3, byrow = TRUE)
    colnames(r) <- c("lon", "lat", "finalbearing")
    return(r[, 1:2, drop = FALSE])
}

# Data can be downloaded from
# http://osm13.openstreetmap.fr/~cquest/openfla/export/communes-20190101-shp.zip
# or
# https://www.data.gouv.fr/en/datasets/decoupage-administratif-communal-francais-issu-d-openstreetmap/
# ("Export simple de janvier 2019 (225Mo)")

# Load shapefile.
# shp <- raster::shapefile("Dropbox/work/crema/communes-20190101-shp/communes-20190101.shp")
# Extract Paris.
paris <- shp[shp$nom == "Paris", ]

# Set distance of points in meters.
dist <- 500

# Extract bounding box from Paris' SpatialPolygonDataFrame.
bbox <- raster::extent(paris)

# Calculate number of points on the vertical axis.
ny <- ceiling(geosphere::distGeo(p1 = c(bbox@xmin, bbox@ymin),
                               p2 = c(bbox@xmin, bbox@ymax)) / dist)
# Calculate maximum number of points on the horizontal axis.
# This needs to be calculated for the lowermost and uppermost horizontal lines
# as the distance between latitudinal lines varies when the longitude changes.
nx <- ceiling(max(geosphere::distGeo(p1 = c(bbox@xmin, bbox@ymin),
                                   p2 = c(bbox@xmax, bbox@ymin)) / dist,
                geosphere::distGeo(p1 = c(bbox@xmin, bbox@ymax),
                                   p2 = c(bbox@xmax, bbox@ymax)) / dist))

# Create result data frame with number of points on vertical axis.
df <- data.frame(ny = 1:ny)

# Calculate coordinates along the vertical axis.
pts <- geosphere::destPoint(p = c(bbox@xmin, bbox@ymin),
                            b = 0, d = dist * (1:ny - 1))
df$x <- pts[, 1]
df$y <- pts[, 2]

# Add points on horizontal axis.
df <- tidyr::crossing(nx = 1:nx, df)

# Calculate coordinates.
pts <- destPoint_v(df$x, df$y, b = 90, 500 * (df$nx - 1))

# Turn coordinates into SpatialPoints.
pts <- SpatialPoints(cbind(pts[, 1], pts[, 2]), proj4string = CRS(proj4string(paris)))

# Cut to boundaries of Paris.
result <- raster::intersect(pts, paris)

# Plot result.
plot(result)
title("Paris in Points")

有点像鱼，不是吗？

关于r - 在多边形中生成规则间隔的点，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/25547826/