1

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) 

Your Answer

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

Browse other questions tagged or ask your own question.