rr-sfr-spterra

divide shapefile into equal sized grids


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.


Solution

  • 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.