I want to divide a shapefile into equal size grids (polygons). This is how I am doing it:
# load shapefile
library(geodata)
library(terra)
can_shp <- geodata::gadm('CAN', level = 1, path = file.path(dir_ls$input, 'shapefile'))
can_shp <- can_shp[can_shp$NAME_1 == 'British Columbia']
# using the extent of the shapefile, create a raster
temp_raster <- terra::rast(xmin = ext(can_shp)[1],
xmax = ext(can_shp)[2],
ymin = ext(can_shp)[3],
ymax = ext(can_shp)[4],
resolution = 0.008333333,
crs = crs(can_shp),
val = 0)
# retain only those raster cells that are within the boundary of the shapefile
temp_raster <- terra::mask(temp_raster, can_shp)
temp_raster <- terra::crop(temp_raster, can_shp)
# convert the resulting raster into grid cells (polygons)
temp_grid <- as.polygons(temp_raster, trunc = FALSE, dissolve = FALSE)
temp_grid
is what I am after which will give me equal grids within the countary boundary but the last step is taking hours and wondered if there's anything I am doing wrong here or if there's any faster way to do this.
Here is your script in more concise form, and for a smaller country that does not need to be downloaded. It shows that your approach works, in principle.
library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))
r <- terra::rast(v, resolution = 0.008333333, vals=1)
r <- terra::mask(r, v)
g <- as.polygons(r, dissolve=FALSE)
Doing the same for British Columbia creates a very large set of polygons that would be hard to manipulate. Keeping your data as a raster is much more efficient. So, to me, the question is why you want these polygons. You would almost certainly be better off without them.