主页 / coding / 问题 / 79074934
Accepted
M. Beausoleil
Asked: 2024-10-10 22:33:40 +0800 CST

如何从 sf 对象的打印头中提取“大地测量 CRS”或“投影 CRS”？

  • 772

将 sf 对象打印到控制台时，有一个带有Geodetic CRS或 的标题Porjected CRS，其中包含 CRS 的良好格式名称。有没有办法提取此信息，而不是

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:
[...]

看看如何st_crs(sfc)$input给予"EPSG:4326"而不是WGS 84。然而，st_crs(nc)$input给予"NAD27"，但st_crs(st_transform(sfc, 3857))$input给予 "EPSG:3857"而不是Projected CRS: WGS 84 / Pseudo-Mercator

地球上是否也存在与栅格数据相当的东西？

2 个回答

  1. 您可以使用wkt元素st_crs()。您将获得一个包含众所周知的文本的字符串；然后是使用正则表达式进行文本提取的问题：

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

    要得到

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

    检查来源，例如

    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"))
#>     }
#> /.../
#> }

    虽然你基本上可以使用自己，但它也可以通过返回的对象sf:::crs_parameters()$Name公开：crsst_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"

    它类似于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)"

    创建于 2024-10-10，使用reprex v2.1.1

    • 2

相关问题