Terra function intersect() and crop() on large SpatVector returns Large list in R

user_19883393 0
Viewed 60

I have what I think is a fairly simple workflow on ArcGIS that I am trying to automatise using r/terra to run different scenarios for the same flow. I am quite an advanced ArcMap user and fluent in R and data manipulation, but 100% new to Terra and this has left me stuck for days!

Context: I have a large dataset of (multipart) polygons: ~2,000 species distributions in Australia including marine ones, with a spatial resolution of 1 or 10km each. There are 150,000 single polygons in total when fully disaggregated.

Goal: I need to be able to calculate metrics related to species for each cell of a given grid. Metrics will include, but will not limited to, species number or area covered by each species. A data.frame containing all the information in the species distribution shapefile within each cell would be the ideal product.

Issue: I tried rasterize() on aggregated data but it did not return the correct species count (I went through this https://github.com/rspatial/terra/issues/553 and was not able to fix, possibly because of all the tiny polygons involved). I chose a less straightforward solution (but more practical for my needs) using intersect(). I end up with a Large list of thousands of elements (which I have no idea how to deal with) instead of a SpatVector (which I ultimately need for spatial data processing). It worked before but I am unable to pinpoint the issue. The same work flow works fine in ArcGIS on subset data.

library(terra)
library(dplyr)

# Load dummy data 
p <- vect(system.file("ex/lux.shp", package="terra"))
v <- p
cell_size <- 0.1 # in Decimal degrees

# Create grid cells
r <- rast(v, res=cell_size)
# Give a value (ID) to each cell and name this value
values(r) <- 1:ncell(r)
names(r) <- "CELLID"
# Transform raster grid to polygons grid
z <- as.polygons(r)
# Intersect with species
u <- intersect(z,v)

This is all I needed for subsequent data analysis (which will look like the line below, only including other metrics than simple species count).

# Create vector with unique species per cell
pa <- aggregate(u, by=c("CELLID","NAME_1"))

Actual data: The data can be downloaded here: http://www.environment.gov.au/fed/catalog/search/resource/details.page?uuid=%7B337B05B6-254E-47AD-A701-C55D9A0435EA%7D

Tried: I've tried many fixes on the base data, including subsetting the data spatially (crop() to a smaller area also returns a Large list), removing a bunch of those "issue" species, disaggregating, aggregating, fixing geometries, etc. I always end up with this large list either at the crop() step or the intersect() step.

Am I missing something obvious?

Sorry for the long post, I tried to include as much as I could.

Thanks in advance.

terra large-data gis polygon r
1 Answers

Turns out I was, indeed, missing something obvious. The code snippet runs perfectly fine, except there was a conflict between library(dplyr) and library(terra) preventing terra functions such as intersect() and crop() from working properly. Duh! Thanks a lot Robert for pointing out the packages!

There are different ways to fix this but I will use the suggestions in this link:

This is the new header of my script:

library(dplyr)
requireNamespace("terra")

Then I use terra:: each time I will use a terra function. This option better suits the analyses that I will perform after that bit of code below.

# Load dummy data 
p <- terra::vect(system.file("ex/lux.shp", package="terra"))
v <- p
cell_size <- 0.1 # in Decimal degrees

# Create grid cells
r <- terra::rast(v, res=cell_size)
# Give a value (ID) to each cell and name this value
terra::values(r) <- 1:terra::ncell(r)
names(r) <- "CELLID"
# Transform raster grid to polygons grid
z <- terra::as.polygons(r)
# Intersect with species
u <- terra::intersect(z,v)
user_198833930
Related