I am updating code from rgdal to use the sf package.

I have a SpatialPolygonsDataFrame that I need to project. In rgdal I used proj4string(). With sf, I tried to use st_crs():

library(sf)
ext <- extent(c(0, 20, 0, 20))
r <- raster(ext, res=1)  
r[] = 1:ncell(r)

# convert the raster to polygon:
Output_Shapefile <- rasterToPolygons(r) 

prj <- "+proj=eqdc +lat_0=17.8333333333333 +lon_0=-66.4333333333333 +lat_1=18.4333333333333 +lat_2=18.0333333333333 +x_0=200000 +y_0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

#project:
#  proj4string(Output_Shapefile) <- prj


st_crs(Output_Shapefile) <- prj

Here is the error I am getting:

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

I am updating code from rgdal to use the sf package.

I have a SpatialPolygonsDataFrame that I need to project. In rgdal I used proj4string(). With sf, I tried to use st_crs():

library(sf)
ext <- extent(c(0, 20, 0, 20))
r <- raster(ext, res=1)  
r[] = 1:ncell(r)

# convert the raster to polygon:
Output_Shapefile <- rasterToPolygons(r) 

prj <- "+proj=eqdc +lat_0=17.8333333333333 +lon_0=-66.4333333333333 +lat_1=18.4333333333333 +lat_2=18.0333333333333 +x_0=200000 +y_0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

#project:
#  proj4string(Output_Shapefile) <- prj


st_crs(Output_Shapefile) <- prj

Here is the error I am getting:

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

• r
• spatial
• r-sf

• sp versions starting with 2.0 are using sf for crs. Did you have any issues with proj4string(Output_Shapefile) <- prj ? – margusl Commented Mar 30 at 15:57
• this is what I am getting: 'Error: unable to find an inherited method for function ‘proj4string<-’ for signature ‘obj = "SpatialPolygonsDataFrame", value = "crs"’' – Ilik Commented Mar 30 at 16:01
• 1 I'm not able to reproduce this with sp-2.2-0 & R 4.4.3 on Windows - gist.github/marguslt/… . Please consider including a complete example with all library() calls. And also your system details, e.g. from sessionInfo(). – margusl Commented Mar 30 at 16:36
• Reproduced the error using R 4.4.3, sf_1.0-19, raster_3.6-31, and sp_2.2-0 using Windows and glad you got a working answer. Two points 1) maybe you're just using raster for reprex purposes (?), but if updating to sf, consider also switching to raster 's replacement terra. The transition is relatively seamless as function naming has for the most part been 'intuitively' migrated to terra e.g. raster::raster() == terra::rast(), raster::extent() == terra::ext() etc. 2) need to project != set projection, so consider editing your question remove that ambiguity. – L Tyrone Commented Mar 30 at 20:16

2 Answers 2

First of all, consider working completely with sf and avoid the use of sp as sp is in maintenance mode right now.

Solution is to convert the SpatialPolygonsDataFrame to sf, assign the crs and back to SpatialPolygonsDataFrame:

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(raster)
#> Cargando paquete requerido: sp
ext <- extent(c(0, 20, 0, 20))
r <- raster(ext, res = 1)
r[] <- 1:ncell(r)

# convert the raster to polygon:
Output_Shapefile <- rasterToPolygons(r)

prj <- "+proj=eqdc +lat_0=17.8333333333333 +lon_0=-66.4333333333333 +lat_1=18.4333333333333 +lat_2=18.0333333333333 +x_0=200000 +y_0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

# project:
#  proj4string(Output_Shapefile) <- prj


st_crs(Output_Shapefile) <- prj
#> Error in UseMethod("st_crs<-"): no applicable method for 'st_crs<-' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialVector')"

And here the solution:

## Solution

# Object with no sp
Output_Shapefile
#> class       : SpatialPolygonsDataFrame 
#> features    : 400 
#> extent      : 0, 20, 0, 20  (xmin, xmax, ymin, ymax)
#> crs         : NA 
#> variables   : 1
#> names       : layer 
#> min values  :     1 
#> max values  :   400

# Instead do
Output_Shapefile_sf <- st_as_sf(Output_Shapefile)
Output_Shapefile_sf
#> Simple feature collection with 400 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: 0 xmax: 20 ymax: 20
#> CRS:           NA
#> First 10 features:
#>    layer                       geometry
#> 1      1 POLYGON ((0 20, 1 20, 1 19,...
#> 2      2 POLYGON ((1 20, 2 20, 2 19,...
#> 3      3 POLYGON ((2 20, 3 20, 3 19,...
#> 4      4 POLYGON ((3 20, 4 20, 4 19,...
#> 5      5 POLYGON ((4 20, 5 20, 5 19,...
#> 6      6 POLYGON ((5 20, 6 20, 6 19,...
#> 7      7 POLYGON ((6 20, 7 20, 7 19,...
#> 8      8 POLYGON ((7 20, 8 20, 8 19,...
#> 9      9 POLYGON ((8 20, 9 20, 9 19,...
#> 10    10 POLYGON ((9 20, 10 20, 10 1...
st_crs(Output_Shapefile_sf) <- prj
Output_Shapefile_sf
#> Simple feature collection with 400 features and 1 field
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: 0 xmax: 20 ymax: 20
#> Projected CRS: +proj=eqdc +lat_0=17.8333333333333 +lon_0=-66.4333333333333 +lat_1=18.4333333333333 +lat_2=18.0333333333333 +x_0=200000 +y_0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
#> First 10 features:
#>    layer                       geometry
#> 1      1 POLYGON ((0 20, 1 20, 1 19,...
#> 2      2 POLYGON ((1 20, 2 20, 2 19,...
#> 3      3 POLYGON ((2 20, 3 20, 3 19,...
#> 4      4 POLYGON ((3 20, 4 20, 4 19,...
#> 5      5 POLYGON ((4 20, 5 20, 5 19,...
#> 6      6 POLYGON ((5 20, 6 20, 6 19,...
#> 7      7 POLYGON ((6 20, 7 20, 7 19,...
#> 8      8 POLYGON ((7 20, 8 20, 8 19,...
#> 9      9 POLYGON ((8 20, 9 20, 9 19,...
#> 10    10 POLYGON ((9 20, 10 20, 10 1...

# Back to sp
Output_Shapefile <- as(Output_Shapefile_sf, "Spatial")

Output_Shapefile
#> class       : SpatialPolygonsDataFrame 
#> features    : 400 
#> extent      : 0, 20, 0, 20  (xmin, xmax, ymin, ymax)
#> crs         : +proj=eqdc +lat_0=17.8333333333333 +lon_0=-66.4333333333333 +lat_1=18.4333333333333 +lat_2=18.0333333333333 +x_0=200000 +y_0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
#> variables   : 1
#> names       : layer 
#> min values  :     1 
#> max values  :   400

^{Created on 2025-03-30 with reprex v2.1.1}

As you are moving to "sf", it makes sense to also leave "raster" behind and use "terra" instead.

library(terra)
ext <- ext(c(0, 20, 0, 20))
prj <- "+proj=eqdc +lat_0=17.8333333333333 +lon_0=-66.4333333333333 +lat_1=18.4333333333333 +lat_2=18.0333333333333 +x_0=200000 +y_0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
r <- rast(ext, res=1, crs=prj)  
values(r) = 1:ncell(r)

# you can do 
pols <- as.polygons(r) 
sf1 <- sf::st_as_sf(pols)
sf1::st_transform(sf1, "+proj=longlat")

#or 
p <- project(pols, "+proj=longlat")
sf2 <- sf::st_as_sf(pols)