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

Question

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

我正在从激光雷达点云中导出各种栅格层，我想将复合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组件。下面的函数是我尝试编写的。我已经用旧的和新的激光雷达数据对它进行了测试，旧的数据有一个简单的(只是水平的)CRS，新数据有一个复合CRS。

我相信一定有更好的方法。

更新2021-11-01:调整该函数接受任何与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)
}
}