Generate regularly spaced points in polygon

Question

Is there a way to generate regularly spaced (e.g., 500 meters apart) points within a polygon using R? I have been trying to use the sp package but can't seem to define a set of points that are spaced a certain distance apart from one another. My aim is to generate the points, then extract their lat/long coordinates into a new dataframe. Any help would be much appreciated! Thanks

You seem to be looking for [Circle Packing](https://en.wikipedia.org/wiki/Circle_packing), specifically, packing circles of radius 250 meters (then any two centers of adjoining circles will be 500 meters apart). — Stephan Kolassa, Aug 28 '14 at 11:43
Maybe `SpatialGrid` or `GridTopology` from package `sp` in combination with `over()` — rcs, Aug 28 '14 at 11:52
Can you just create a regular array of points on 500-m spacings, and then use one of the "within polygon" functions such as `point.in.polygon` to reject unwanted elements of the array ? — Carl Witthoft, Aug 28 '14 at 12:43
Thanks for the comments. Carl, how would I go about generating a regular array of points? — Ross, Aug 28 '14 at 13:00
Maybe `cbind(rep(seq(0,5000,500),each=10),rep(seq(0,5000,500),times=10))` ? That is, a matrix w/ a coordinate pair in each row. — Carl Witthoft, Aug 28 '14 at 14:24
On a plane, can't you just create a rectangular lattice or triangular lattice, then mask with the polygon you want? On a sphere (you mention lat/long), any solution will be approximate. Over what scale is your polygon: 100's of meters? 100's of kilometers? — Dan H, Apr 07 '18 at 03:15

score 2 · Answer 1 · answered May 29 '19 at 23:10

Quite straight forward and almost out-of-the-box.

As OP did not share data, buckle up, put your seats in a vertical position and let us fly to Paris. There, we will adapt a geosphere function, and with its help we will divide up Paris' shape into lon / lat coordinates that are 500 meters apart each (vertically and horizontally).

# 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")

Kind of looks like a fish, doesn't it?

geosphere functions are vectorized. In stead of creating `destPoint_v` , you can directly use `destPoint`, as in `geosphere::pts <- destPoint(df[, c("x", "y")], b = 90, 500 * (df$nx - 1))` — Robert Hijmans, May 30 '19 at 00:53

Robert Hijmans · Answer 2 · 2019-05-30T03:29:16.977

Here is a way to do assuming you have a lonlat polygon by first transforming it to a planar crs (not as nifty as Roman's solution with destPoint).

Packages and example data

library(raster)
library(rgdal)
p <- shapefile(system.file("external/lux.shp", package="raster"))[1,]

Transform to planar crs (pick one that matches your data!)

putm <- spTransform(p, "+proj=utm +zone=32 +datum=WGS84")

Create a raster with 500 m resolution, rasterize the polygon and transform to points

r <- raster(putm, res=500)
r <- rasterize(putm, r)
pts <- rasterToPoints(r, spatial=TRUE)

Transform the points to lon/lat and plot the results

pts_lonlat <- spTransform(pts, "+proj=longlat +datum=WGS84")
result <- coordinates(pts_lonlat)

plot(p)
points(result, pch="+", cex=.5)

(looks like an elephant)

Generate regularly spaced points in polygon

2 Answers2

Linked