I need to create a grid inside a shapefile with a fixed the number of cells n (inside the shapefile). I tried a few stuff see bellow but nothing worked so far. Any ideas?
The data shapefile can be found here:
https://geodacenter.github.io/data-and-lab/scotlip/
library(maptools)
library(sp)
library(raster)
library(rgeos)
scotlip <- readShapeSpatial("~/Desktop/scotlip/scotdistricts.shp")
plot(scotlip)
n <- 100 # the number of cells I want include the shapefile.
# Attempt 1
set.seed(12)
sppoint <- spsample(scotlip, n = n, type = "regular")
nrow(sppoint@coords)
# first problem: 92 instead of 100 points
points(sppoint, col = "red", pch = 19, cex = .5)
grid <- SpatialPixels(sppoint)
grid <- as(grid, "SpatialPoygon")
plot(grid)
plot(scotlip, border = "red", add = T)
# second problem: it doesnt cover the shapefile, and thus the areas differ
# Attempt 2
xmin <- scotlip@bbox[1,1] - 10000
xmax <- scotlip@bbox[1,2] + 10000
ymin <- scotlip@bbox[2,1] - 10000
ymax <- scotlip@bbox[2,2] + 10000
ext <- extent(c(xmin, xmax, ymin, ymax))
# try fixing the resolution based on the number of grid cells required to cover the area, but
# this doesn't work too, because of the grid cells that intersect with the land.
tot.area <- sum(area(scotlip))
res(r) <- sqrt(tot.area/n)
grid <- SpatialPixelsDataFrame(SpatialPoints(coordinates(r)),
data = data.frame(layer = 1:nrow(coordinates(r))))
grid <- as(grid, "SpatialPolygonsDataFrame")
plot(grid)
plot(scotlip, border = "red", add = T)
int <- intersect(grid, sctolip)
sum(area(int)) - sum(area(scotlip))
plot(int)
nrow(int)