0

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.


Code corrected after @Spacedman’s advice.

8

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.

3

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
#> Loading required package: sp
#> 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)

1
  • The only thing here is that you need to project the shapefile in order to make st_sample works properly. – dieghernan Feb 15 '20 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.