`spatialBlock` fails for data in EPSG:3035

[AMBarbosa]

Hi, I have species occurrences in EPSG 3035 (ETRS89-extended / LAEA Europe), which for some reason causes an error in spatialBlock. Here's some code based on the help file examples:

awt <- raster::brick(system.file("extdata", "awt.grd", package = "blockCV"))
PA <- read.csv(system.file("extdata", "PA.csv", package = "blockCV"))
pa_data <- sf::st_as_sf(PA, coords = c("x", "y"), crs = raster::crs(awt))

pa_data <- sf::st_transform(pa_data, 3035)  # to show what happens if data are in this CRS

sb1 <- blockCV::spatialBlock(speciesData = pa_data, theRange = 70000)
# Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared,  : 
#  st_crs(x) == st_crs(y) is not TRUE

I'm not sure if it's because this is a "custom" European CRS. Indeed, if I use the "plain" LAEA projection, I don't get the error:

pa_data <- sf::st_transform(pa_data, "+proj=laea")  # to show what happens if data are in this CRS
sb1 <- blockCV::spatialBlock(speciesData = pa_data, theRange = 70000)
# ok

For now I'm circumventing this by projecting my data before using spatialBlock, but a fix would be appreciated. Cheers!

[rvalavi]

Hi @AMBarbosa , thanks for pointing to the issue. The data stored in the blockCV package is located in eastern Australia and changing projection to an EU system will completely distort them. Do have the same error on data from EU with EPSG 3035? I know where the issue is coming from; just want to make sure it's not because what I mentioned.

[AMBarbosa]

Hi, yes I found the problem with my EU species data, which came originally in this CRS.

[rvalavi]

@AMBarbosa I added a minor update to the package. Re-download the package from GitHub and check if the problem is solved. You'll probably get some warnings that are because of the new update of projection systems in R. This is a temporary fix and I'll go into more depth when I want to update blockCV with the terra package. Please let me know if you still have the issue.

[AMBarbosa]

Hi, I've just had a chance to install the latest GitHub version of 'blockCV' and try again with my European data. If I don't use the 'rasterLayer' argument, the error no longer occurs, although there's a "_Warning message: st_crs<- : replacing crs does not reproject data; use sttransform for that". However, if I do use a 'rasterLayer' for visualisation (in the same CRS as the species data, i.e. EPSG:3035), I still get "_Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, : st_crs(x) == stcrs(y) is not TRUE". Here's some code with a European species:

library(terra)
library(geodata)
library(blockCV)

# get species occurrence data in EPSG:3035
download.file("https://easin.jrc.ec.europa.eu/easin/docs/baseline/gis/Baseline2019.zip", destfile = "Baseline2019.zip")
unzip("Baseline2019.zip", exdir = "Baseline2019")
list.files("Baseline2019")
bl <- vect("Baseline2019.shp")

# choose a species:
unique(bl$SpeciesNam)
occdata <- bl[bl$SpeciesNam == "Lepomis gibbosus", ]
occdata <- centroids(occdata)
plot(occdata)

# download a raster map, crop to species area and project to EPSG:3035
elev <- elevation_global(res = 10, path = tempdir())
elev <- crop(elev, project(occdata, crs(elev)))
elev <- project(elev, crs(occdata))
plot(elev)
plot(occdata, add = TRUE)

# compute spatial blocks:
set.seed(2)
blocks <- spatialBlock(speciesData = as(occdata, "Spatial"), theRange = 200000, k = 5)
# OK except "Warning message: st_crs<- : replacing crs does not reproject data; use st_transform for that"

# however, if a 'rasterLayer' is used:
set.seed(2)
blocks <- spatialBlock(speciesData = as(occdata, "Spatial"), rasterLayer = raster::raster(elev), theRange = 200000, k = 5)
# Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared,  : st_crs(x) == st_crs(y) is not TRUE

# check that maps do overlap correctly:
plot(elev)
plot(blocks$blocks, add = TRUE)
plot(occdata, cex = 0.1, add = TRUE)

(NOTE: with set.seed(1), 'spatialBlock' gives me a different error which I've posted under a separate issue here)

[benikm91]

I had the same error "st_crs(x) == st_crs(y) is not TRUE" when trying to make an older project run. With blockCV 2.1.4 it happened even without rasterLayer. With blockCV 2.1.5 it is only happening with parameter rasterLayer (I think it is because my value for rasterLayer is not a "sf" object and fix from 2.1.5 only sets crs if parameter is a sf object).

I found out that the problem is the raster package. This is where the actual exception occures. If I down grade "raster" to the version around a year ago (where my project was originally developed), everything works fine.

The version I use is "raster" version "3.4-5".

I installed it with devtools:

devtools::install_version("raster", version="3.4-5")

[rvalavi]

Thanks for reporting this @benikm91 I try to check this tonight and fix the bogs. This is all about the newer updates of the sf package and the updates of the PROJ 6 library for spatial reference systems.

[Lukpawlik]

Having the same problem today when using sf object with points from France in crs=3035. I've got this error: Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, : zmienna st_crs(x) == st_crs(y) nie ma wartości TRUE blockCV version 2.1.4

Have you got any upate on this issue? Thanks, Lukasz

[Lukpawlik]

After updating to blockCV 2.1.5 it works now for function blockCV::spatialBlock(speciesData = df09, species = "damage", theRange = 50000, k = 10, selection = "random")

[rvalavi]

Please try blockCV v3.0, the issue must be resolved. I'm closing this issue.