CODEWITHSUNDEEP
rr-sfr-spterra
89_Simple
89_Simple

Reputation: 3805

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.

Upvotes: 2

Views: 1056

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47071

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.

Upvotes: 3

caldwellst
caldwellst

Reputation: 5956

You can do it easily with the sf package, which allows easy grid creation and reduction of polygons to those within boundaries. Here's an example using the sample North Carolina data from the package itself. ggplot2 just used for clearer visuals.

library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.1, PROJ 9.0.1; sf_use_s2() is TRUE
library(ggplot2)

# sample North Carolina data
nc <- st_read(system.file("gpkg/nc.gpkg", package="sf"))
#> Reading layer `nc.gpkg' from data source 
#>   `/usr/local/lib/R/site-library/sf/gpkg/nc.gpkg' using driver `GPKG'
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27

plot(nc, max.plot = 1)

enter image description here


# make grid from the extents
nc_grid <- st_make_grid(nc)

ggplot() +
  geom_sf(data = nc, fill = "grey", color = NA) +
  geom_sf(data = nc_grid, fill = NA) +
  theme_void()

enter image description here


# filter the grid to cells within NC
nc_filter <- st_within(nc_grid, st_union(nc)) %>% lengths > 0

ggplot() +
  geom_sf(data = nc, fill = "grey", color = NA) +
  geom_sf(data = nc_grid[which(nc_filter)], fill = NA) +
  theme_void()

enter image description here

Upvotes: 1

Related Questions