Set projection of SpatialPolygonsDataFrame using the sf package

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; sfuses2() is TRUE
library(raster)
#> Cargando paquete requerido: sp
ext  names       : layer 
#> min values  :     1 
#> max values  :   400
Instead do
OutputShapefilesf  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...
stcrs(OutputShapefilesf)  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 +lat0=17.8333333333333 +lon0=-66.4333333333333 +lat1=18.4333333333333 +lat2=18.0333333333333 +x0=200000 +y0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +nodefs
#> 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
OutputShapefile  class       : SpatialPolygonsDataFrame 
#> features    : 400 
#> extent      : 0, 20, 0, 20  (xmin, xmax, ymin, ymax)
#> crs         : +proj=eqdc +lat0=17.8333333333333 +lon0=-66.4333333333333 +lat1=18.4333333333333 +lat2=18.0333333333333 +x0=200000 +y0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +nodefs 
#> variables   : 1
#> names       : layer 
#> min values  :     1 
#> max values  :   400

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