0

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)

8
  • 1
    What is the aim of converting polygons to rasters in your case?
    – Babel
    Jan 12 at 10:00
  • @Babel The goal is to be able to perform calculations of occurence "density". A sort of heatmap
    – Spectron
    Jan 12 at 10:07
  • @BERA Yeah, it freezes
    – Spectron
    Jan 12 at 10:08
  • 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 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: gis.stackexchange.com/a/419732/88814 or 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 at 10:59

1 Answer 1

0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.