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:
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)