Início / coding / Perguntas / 79074934
Accepted
M. Beausoleil
Asked: 2024-10-10 22:33:40 +0800 CST

Como extrair o `Geodetic CRS` ou `Projected CRS` do cabeçalho de impressão do objeto sf?

  • 772

Ao imprimir objetos sf no console, há um cabeçalho com o Geodetic CRSou Porjected CRSque tem um nome de formato agradável para o CRS. Existe uma maneira de extrair essas informações em vez de

p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
sfc = st_sfc(p1, p2, crs = 4326)
sfc
Geometry set for 2 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -30 ymin: 20 xmax: 7 ymax: 52
Geodetic CRS:  WGS 84                  <- this line
POINT (7 52)
POINT (-30 20)

st_transform(sfc, 3857)
Geometry set for 2 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
Projected CRS: WGS 84 / Pseudo-Mercator                  <- this line
POINT (779236.4 6800125)
POINT (-3339585 2273031)

nc = st_read(system.file("shape/nc.shp", package="sf"))
nc
Simple feature collection with 100 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27                  <- this line
First 10 features:
[...]

Veja como st_crs(sfc)$inputdá "EPSG:4326"e não WGS 84. No entanto, st_crs(nc)$inputdá "NAD27", mas st_crs(st_transform(sfc, 3857))$inputdá "EPSG:3857"e nãoProjected CRS: WGS 84 / Pseudo-Mercator

Existe também o equivalente em terra para dados raster?

2 respostas

  1. Você pode usar o wktelemento from st_crs(). Você obtém uma string com o texto bem conhecido ; então é uma questão de extração de texto com uma expressão regular :

    library(stringr)
str_extract(st_crs(st_transform(sfc, 3857))$wkt, "(?<=PROJCRS\\[\\\").*(?=\\\")")

    obter

    [1] "WGS 84 / Pseudo-Mercator"
    • 3
  2. Best Answer

    Verifique a fonte, por exemplo

    library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

sloop::s3_dispatch(print(st_sfc()))
#>    print.sfc_GEOMETRY
#> => print.sfc
#>  * print.default

sloop::s3_get_method(print.sfc)
#> function (x, ..., n = 5L, what = "Geometry set for", append = "") 
#> {
#> /.../
#>     crs = st_crs(x)
#>     if (is.na(crs)) 
#>         cat(paste0("CRS:           NA\n"))
#>     else {
#>         p = crs_parameters(crs)
#>         if (p$Name == "unknown") {
#>             if (is.character(crs$input) && !is.na(crs$input) && 
#>                 crs$input != "unknown") 
#>                 p$Name = crs$input
#>             else p$Name = crs$proj4string
#>         }
#>         if (p$IsGeographic) 
#>             cat(paste0("Geodetic CRS:  ", p$Name, "\n"))
#>         else cat(paste0("Projected CRS: ", p$Name, "\n"))
#>     }
#> /.../
#> }

    Embora você possa basicamente usar sf:::crs_parameters()$Namevocê mesmo, ele também é exposto por meio crsdo objeto retornado por st_crs():

    st_crs(st_sfc(crs = 4326))$Name
#> [1] "WGS 84"
st_crs(st_sfc(crs = 3857))$Name
#> [1] "WGS 84 / Pseudo-Mercator"
st_crs(st_sfc(crs = 4267))$Name
#> [1] "NAD27"

    É semelhante com terra:

    x <- terra::rast(nrows=108, ncols=21, xmin=0, xmax=10)
x
#> class       : SpatRaster 
#> dimensions  : 108, 21, 1  (nrow, ncol, nlyr)
#> resolution  : 0.4761905, 1.666667  (x, y)
#> extent      : 0, 10, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
terra::crs(x, describe = TRUE)$name
#> [1] "WGS 84 (CRS84)"

# sf::st_crs() would also work:
sf::st_crs(x)$Name
#> [1] "WGS 84 (CRS84)"

    Criado em 2024-10-10 com reprex v2.1.1

    • 2

relate perguntas