Rasterizing more than 11000 polygons in QGIS

Question

I'm using QGIS 3.28.3. I have more than 11000 species ranges, which I want to transform into binary rasters (1 for species occurence, 0 for species non occurence), with one raster for each species.

I tried the rasterize tool in batch mode, but QGIS can't handle that much information, it freezes. Nor can it work with the data split in smaller batches.

Does anyone have a solution? I'm not very familiar with Python, but can work with R (although I'm not sure how to make a batch proccess in R)

What is the aim of converting polygons to rasters in your case? — Babel, Jan 12 '24 at 10:00
@Babel The goal is to be able to perform calculations of occurence "density". A sort of heatmap — Spectron, Jan 12 '24 at 10:07
The species ranges are all in one vector layer? Can you create a grid, intersect it with your vector layer, then build a model with a vector features input to create one raster per grid cell. Then merge all raster outputs? — BERA, Jan 12 '24 at 10:19
I asked because you might want to consider Union geoprocessing tool to identify areas where polygons overlap and to create a separate polygon of each area with a different number of overlaps. Something like this: https://gis.stackexchange.com/a/419732/88814 or https://gis.stackexchange.com/a/392905/88814 And: did you test if your polygons don't have errors and are valid? That might be a problem — Babel, Jan 12 '24 at 10:59
@user30184 X: 4500 Y: 2126 (whole world divided in 10x10km approximately — Spectron, Jan 14 '24 at 22:11
If your data is fairly sparse - which I'm guessing it is - you will be likely better off sticking with vector data, as others have suggested. What's the total number of occurrences across the 11,000 species? Point data (one point per 10kmx10km square) may be a good option depending on your data size. — Tom Brennan, Jan 15 '24 at 05:37

score 0 · Answer 1 · answered Jan 21 '24 at 14:39

My approach to map species occurrences is: to download GBIF data for species and their synonyms, merge them and run a cleanup. Within R there are rgbif and CoordinateCleaner packages as a helpers. Then the process looks like below (sticking to vectorized data/maps)

As first step we have to prepare the grid, which be used as a occurrence counter:

sf::sf_use_s2(FALSE)
data(World, package = "tmap")
all <- World |>
  dplyr::summarize() |>
  terra::vect() |>
  terra::fillHoles() |>
  sf::st_as_sf()
grid <- sf::st_make_grid(all, cellsize = 3, what = "polygons", square = FALSE) |>
  sf::st_intersection(all) |>
  sf::st_as_sf()
grid |>
  plot(lwd = 0.3)

Having the grid ready we can import species data:


p <- sf::st_read("/home/sapi/prunus_serotina_ehrh.gpkg",
                 layer = "prunus_serotina_ehrh")
#> Simple feature collection with 183634 features and 21 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -123.3403 ymin: -43.66203 xmax: 177.7184 ymax: 60.67848
#> Geodetic CRS:  WGS 84

Let's show them on the map:

tmap::tm_shape(grid) +
  tmap::tm_polygons() +
  tmap::tm_shape(p) +
  tmap::tm_dots()

Now, we are counting how many occurrences are in particular grid. A small explanation: st_intersects() returns FALSE/TRUE, and lengths() of it counts how many TRUE values was there.

grid$num_points <- lengths(sf::st_intersects(grid, p))

And time to draw them on the map. As you are interested only to show presence/absence without counts it's enough to show them as gray polygons without any scale:

g <- grid |>
  subset(num_points > 0)
tmap::tm_shape(all) +
  tmap::tm_borders() +
  tmap::tm_shape(g) +
  tmap::tm_polygons()

If your data (species range) is in form of polygons, then the similar approach will work, as long as the polygons will intersects with grid. Just make sure both have the same CRS.

If your desired output is a raster, then create one, like:

r <- terra::rast(nrows=180, ncols=360, xmin=-180, xmax=180, 
     ymin=-90, ymax=90, crs = "EPSG:4326", vals = NA)

and rasterize your data:

s <- terra::rasterize(terra::vect(p), r)
terra::plot(s)

^{Created on 2024-01-21 with reprex v2.0.2}

Rasterizing of polygons works the same way, just have a look on cover option in terra::rasterize() function. And maybe touches as well.

If there is 11000 species, create a list of input files with list.files() function and either create a function for lapply(), either a for() {} loop. Plot the resulting graphics with png() function + dev.off(). Or with tmap::save() or whatever you will use.

Rasterizing more than 11000 polygons in QGIS

1 Answers1