# Get approximately n points inside a polygon sf R

I need to get approximately n points in a pattern as regular as possible inside a polygon.

Current code:

library(tidyverse)
library(sf)
library(raster)
library(ggplot2)

n_points_wanted = 200

get_n_plants_grid <- function(shp, n_points) {

area_box <- st_bbox(shp) %>% st_as_sfc() %>% st_area()
area_shp <- shp %>% st_area() %>% sum()

n_grid <- ((area_box / (area_box - area_shp  ) )*n_points) %>% sqrt()

n_grid <- n_grid %>% as.numeric() %>% round()

return (  n_grid )

}

shp <- getData('GADM', country = 'aut', level = 0) %>% st_as_sf()

n_grid <- get_n_plants_grid(shp, n_points_wanted)

grid <- st_make_grid(shp, n=c(n_grid, n_grid),
what = "centers",
square=TRUE) %>% st_intersection(shp)

ggplot() +
geom_sf(data = shp) +
geom_sf(data = grid, col ="red")

print(length(grid))

In this example I got 393 points inside the polygon but wanted 200.

STOP USING SO MANY PIPES!!!

n_grid <- (area_box / (area_box - area_shp  ) )*n_points %>% sqrt()

is doing

n_grid <- (area_box/(area_box - area_shp)) * (n_points %>% sqrt())

in other words its square-rooting n_points. Fix that and I get exactly 200 points.

The pipe operator has high priority:

> 3 * 9 %>% sqrt()
[1] 9
> 3 * (9 %>% sqrt())
[1] 9
> (3 * 9) %>% sqrt()
[1] 5.196152

Pipes. You don't need them, you can always rewrite either with nested function calls or by saving intermediate results. They do not always add clarity or make for tidier code.

What about st_sample(shp, n_points_wanted, type="regular" )? It gives you approximately 200 points, as it has some degree of randomness:

library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 3.5.3
#> Warning: package 'ggplot2' was built under R version 3.5.3
#> Warning: package 'tibble' was built under R version 3.5.3
#> Warning: package 'tidyr' was built under R version 3.5.3
#> Warning: package 'readr' was built under R version 3.5.3
#> Warning: package 'purrr' was built under R version 3.5.3
#> Warning: package 'dplyr' was built under R version 3.5.3
#> Warning: package 'stringr' was built under R version 3.5.3
#> Warning: package 'forcats' was built under R version 3.5.3
library(sf)
#> Warning: package 'sf' was built under R version 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(raster)
#> Warning: package 'raster' was built under R version 3.5.3
#> Warning: package 'sp' was built under R version 3.5.3
#>
#> Attaching package: 'raster'
#> The following object is masked from 'package:dplyr':
#>
#>     select
#> The following object is masked from 'package:tidyr':
#>
#>     extract
library(ggplot2)

n_points_wanted = 200
shp <- getData('GADM', country = 'aut', level = 0) %>% st_as_sf()
shp=st_transform(shp,3857)
set.seed(1234)
grid <- st_sample(shp,n_points_wanted,type="regular")

ggplot() +
geom_sf(data = shp) +
geom_sf(data = grid, col ="red")

print(length(grid))
#> [1] 199

Created on 2020-02-15 by the reprex package (v0.3.0)

• The only thing here is that you need to project the shapefile in order to make st_sample works properly. Commented Feb 15, 2020 at 20:16