Início / coding / Perguntas / 79544779
Accepted
Ilik
Asked: 2025-03-30 23:08:28 +0800 CST

Defina a projeção de SpatialPolygonsDataFrame usando o pacote sf

  • 772

Estou atualizando o código rgdalpara usar o sfpacote.

Tenho um SpatialPolygonsDataFrame que preciso projetar. Em rgdaleu usei proj4string(). Com sf, tentei usar 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

Aqui está o erro que estou recebendo:

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

2 respostas

  1. Best Answer

    Primeiro, considere trabalhar completamente com sf e evite usar sp, pois o sp está em modo de manutenção no momento.

    A solução é converter SpatialPolygonsDataFramepara sf, atribuir o crs e voltar para 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')"

    E aqui está a solução:

    
## 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

    Criado em 2025-03-30 com reprex v2.1.1

    • 1

  2. Como você está migrando para "sf", faz sentido também deixar "raster" para trás e usar "terra".

    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)
    • 1

relate perguntas