R – Creating a Regular Sampling Grid with Specific Point Distances

rspatial-analysisterra

I wish to create a regular sampling grid within m spatial polygon. Similar question is answered as Creating regular soil sampling grids using QGIS but I wish to use R.

In terra package, I can create spatial sample using spatSample(vector,numberOfPoints, method="regular") within a SpatVector. But, this creates the maximum number of points, but ignores the exact distance between the two points.

I was wondering to use two ways ho to derive my 10×10 m point grid:

  • get the extent of my polygon (eg using SpatExt) and estimate the number of points scaling its size into 10 m chunks: ((round(b1$radius)/10)*2+1)^2 – this however leads to not exactly 10 m, but a bit different number (distance of 8.9m, ..)
  • create a centroid of the buffer and add 10 m in all cardinal directions (and on the sides)? this seems like creating nicely set distance, but I am not here sure where to start.

What would be your suggestions?

Dummy example:

   ## SpatVector
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
#sample geometries 
i <- sample(nrow(v), 1)
vv <- v[i,]

# convert to projected system:
vv_3035 <- project(vv,"EPSG:3035") 

# calculate the assumed radius to estimate the number of points at 10 m distance:
vv_3035$radius = sqrt(vv_3035$AREA/pi) 


# Create a spatial sample in regular pattern: 10 indicates total number of points
# not the distance between them!
samp10 <- spatSample(vv_3035, 10, method="regular")
samp_area <- spatSample(vv_3035, ((round(vv_3035$radius)/10)*2+1)^2, method="regular")

Regular sampling grid with 10 points: but how to set up the exact distance between points?

Red: grid for 10 points, Blue: distance by approximation with radius

Best Answer

I would create a SpatRaster using the SpatVector as template, and set the resolution to 10 meters. Then you got a Raster with the desired gap that you can transform to points if needed.

On this example I used a distance of 1000 meters for demo purposes, just set it to 10:

## SpatVector
library(terra)
#> terra 1.5.21

f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)

#sample geometries 
i <- sample(nrow(v), 1)
vv <- v[i,]

# convert to projected system:
vv_3035 <- project(vv,"EPSG:3035") 

# Create a raster with the desired res
# In this example, for demo, I used 1000 m intead of 10 m
template <- rast(vv_3035, resolution = c(1000,1000))

# Fake values, for demo only
values(template) <- 1:ncell(template)

# Mask to polygon
template <- mask(template, vv_3035)


plot(vv_3035)
plot(template, add=TRUE, alpha=0.5)



# To points: This is now a vector
points <- as.points(template, na.rm = TRUE)

# Intersect to keep only points on the shape
points <- points[vv_3035]

# Check
distance(points[20], points[21])
#>      [,1]
#> [1,] 1000


plot(vv_3035)
plot(points, add=TRUE)

Created on 2022-05-23 by the reprex package (v2.0.1)

Related Question