I have a shapefile ‘shp’ with boundaries, and I have a data .csv file with coordinates and some values. I want to join the df to shp, so that each coordinate/observation would be in the correct zone/boundary.

What I’ve done, but I get an error:

shape_df <- st_join(shp, df, join = st_intersects)

(st_contains gives the same error)

Data types (same for both):

> class(shp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

 > class(df)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

> st_crs(shp)
Coordinate Reference System:
  User input: OSGB36 / British National Grid 
  wkt:
PROJCRS["OSGB36 / British National Grid",

>st_crs(pollutants2)
Coordinate Reference System:
  User input: OSGB36 / British National Grid 
  wkt:
PROJCRS["OSGB36 / British National Grid",
    BASEGEOGCRS["OSGB36",

Error:

Error in UseMethod("st_join") : 
  no applicable method for 'st_join' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial')"

What am I doing wrong? How do I solve this?