r - Set projection of SpatialPolygonsDataFrame using the sf package - Stack Overflow

admin管理员组

文章数量:1376572

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

Share Improve this question Follow edited Mar 30 at 19:51 L Tyrone 7,5182323 gold badges2828 silver badges4141 bronze badges asked Mar 30 at 15:08 IlikIlik 18711 silver badge1111 bronze badges 4

• 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

Add a comment  | 

2 Answers 2

Sorted by: Reset to default 1

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)

本文标签： rSet projection of SpatialPolygonsDataFrame using the sf packageStack Overflow