Extract biogeographical genions based on coodinates-CodePudding

I would like to extract the biogeographical for each sample sites based on the coordinates (i.e. latitude and longitude) Searching through different resources I found a code that might work, unfortunately, I am getting an error. I obtained the shapefile of the biogeographical regions from: https://www.eea.europa.eu/data-and-maps/data/biogeographical-regions-europe-3

pg <- rgdal::readOGR(dsn = "BiogeoRegions2016.shp",
                     layer = "BiogeoRegions2016", p4s = " init=epsg:3035") #load data
proj4string(pg) #Check projection
#plot(pg) #Plot

To not overload the code I just used some example point:

# example points
pt <- data.frame(x = c(14,28), y = c(41, 59), id = 1:2) #example points
coordinates(pt) <- ~x y #Transform to spatial coordinates
proj4string(pt) <- CRS(" init=epsg:4326") #projection
proj4string(pt) #projection

#  project to CRS of pg
pt <- spTransform(pt, CRS(" init=epsg:3035")) #Transform the point to the same projection
over(pt, pg) #This function should report the biogeographical region but I got the next error:
#Error message:
Error in .local(x, y, returnList, fn, ...) : 
  identicalCRS(x, y) is not TRUE

More info at: https://stat.ethz.ch/pipermail/r-sig-geo/2015-February/022329.html

Thanks so much in advance for your help

CodePudding user response：

This error Error in .local(x, y, returnList, fn, ...) : identicalCRS(x, y) is not TRUE suggests the CRS assignment and/or reprojection didn't work as required.

If you don't mind using package {sf}, here's an approach:

library(sf)
library(dplyr)

pt <- 
  data.frame(x = c(14,28), y = c(41, 59), id = 1:2) |>
  ## add geometry column:
  st_as_sf(coords = c('x', 'y')) |>
  ## set CRS:
  st_set_crs(4326) |>
  ## transform to EEA projection:
  st_transform(3035)

## read in EEA shapefile:
bioregions <- sf::read_sf('path/to/BiogeoRegions2016.shp')

pt |>
rowwise() |>
mutate(bioregion_index = unlist(st_within(geometry, bioregions)), ## see note
       bioregion_name = bioregions[bioregion_index,]$short_name
       )

note: we have to unlist, because st_within returns a list of row matches, even if there's only a single match.

output:

Simple feature collection with 2 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 4658848 ymin: 1997901 xmax: 5344678 ymax: 4119297
Projected CRS: ETRS89-extended / LAEA Europe
# A tibble: 2 x 4
# Rowwise: 
     id          geometry bioregion_index bioregion_name
* <int>       <POINT [m]>           <int> <chr>         
1     1 (4658848 1997901)               9 mediterranean 
2     2 (5344678 4119297)               6 boreal