是否有 R 函数来解析复合坐标参考系？

Question

我最近开始使用机载 LiDAR 数据，该数据具有由水平（投影）和垂直分量组成的复合坐标参考系统。下面显示了一个示例，其中包含从 WKT 描述创建复合 CRS 对象的代码。

我正在从 LiDAR 点云中导出各种栅格图层，我只想将复合 CRS 的水平分量分配给这些图层（示例中为 EPSG:7856）。有谁知道现有的包函数可以可靠地提取水平 PROJCRS 组件，即允许各种新旧 CRS 定义？

2021-11-01 更新：调整了 WKT 字符串的原始示例，以提供在 R 中创建复合 CRS 对象的代码。

# Create a compound CRS object of the type used for
# publicly available LiDAR point cloud data in Australia.
# Requires the glue and sf packages.
#
wkt <- glue::glue('COMPOUNDCRS["GDA2020 / MGA zone 56 + AHD height - AUSGeoid2020 (Meters) (with axis order normalized for visualization) (with axis order normalized for visualization)",
    PROJCRS["GDA2020 / MGA zone 56",
        BASEGEOGCRS["GDA2020",
            DATUM["Geocentric Datum of Australia 2020",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",7844]],
        CONVERSION["UTM zone 56S",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",153,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",0.9996,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",500000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",10000000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]]],
        CS[Cartesian,2],
            AXIS["easting",east,
                ORDER[1],
                LENGTHUNIT["metre",1]],
            AXIS["northing",north,
                ORDER[2],
                LENGTHUNIT["metre",1]],
        ID["EPSG",7856]],
    VERTCRS["AHD height - AUSGeoid2020 (Meters)",
        VDATUM["Australian Height Datum"],
        CS[vertical,1],
            AXIS["gravity-related height",up,
                LENGTHUNIT["metre",1]],
        ID["EPSG",5711]]]')

# Create the compound CRS object
compound_crs <- st_crs(wkt)

Answer 1

在 R 空间丛林中进行了大量拖网之后，我还没有找到提取水平 CRS 分量的现有函数。下面的函数是我尝试写的一个。我已经使用新旧 LiDAR 数据对其进行了测试，其中旧数据具有简单（仅水平）CRS，而新数据具有复合 CRS。

我确信一定有一个现有的更好的方法。

2021 年 11 月 1 日更新：调整了函数以接受与 sf::st_crs() 函数兼容的任何空间对象。

# Retrieve the horizontal component of a compound CRS.
# The object x can be an 'sf' package 'crs' object or any
# spatial object from which a CRS can be queried using the
# sf::st_crs function.
#
get_horizontal_crs <- function(x) {
  xcrs <- sf::st_crs(x)
  if (is.na(xcrs)) stop("No CRS defined")

  wkt <- sf::st_as_text(xcrs)
  
  if (!grepl("COMPD_CS", wkt)) {
    # Should just be a horizontal CRS - simply return it
    xcrs
  } else {
    # Extract the horizontal component
    i <- regexpr("PROJCS\\[", wkt)
    wkt <- substring(wkt, i)
    
    # Match square brackets to discard any trailing
    # component (e.g. the vertical CRS)
    wkt_chars <- strsplit(wkt, "")[[1]]
    level <- 1
    k <- match("[", wkt_chars)
    while (level > 0) {
      k <- k + 1
      if (wkt_chars[k] == '[') {
        level <- level + 1
      } else if (wkt_chars[k] == ']') {
        level <- level - 1
      }
    }
    
    wkt <- substring(wkt, 1, k)
    
    sf::st_crs(wkt)
  }
}