从 R 中的空间数据中提取测量单位（度、米等）

Posted 2023-02-14

【中文标题】从 R 中的空间数据中提取测量单位（度、米等）

【英文标题】：Extracting the measurement unit (degrees, metres, etc.) from spatial data in R

【发布时间】：2022-01-22 21:48:25

【问题描述】：

我想从 R 中的空间对象中提取测量单位（十进制度、米、英尺等）。例如，如果我有一个使用 WGS84 坐标参考系统的 SF 数据框（ EPSG：4326），我希望能够确定坐标以十进制度数指定。同样，我希望能够确定 UTM 坐标（例如 EPSG:32615）以米为单位。

我尝试使用 sf 包中的 st_crs() 函数，它以众所周知的文本格式返回坐标参考系统。但是，我很难确定从该众所周知的文本中提取测量单位的正则表达式是否能在各种坐标系中可靠地运行。

是否有返回空间对象测量单位的现有函数？

例如，以下代码生成使用 WGS84 坐标系的 SF 数据框：

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1

cities <- st_sf(city = "London", geometry = st_sfc(st_point(c(-0.1276, 51.5072))), crs = 4326)

cities
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -0.1276 ymin: 51.5072 xmax: -0.1276 ymax: 51.5072
#> Geodetic CRS:  WGS 84
#>     city                geometry
#> 1 London POINT (-0.1276 51.5072)

^{由reprex package (v2.0.1) 于 2021 年 12 月 21 日创建}

理想情况下，我正在寻找一个函数，该函数允许我确定此数据集的空间单位是十进制度，例如如果该函数被称为 st_crs_unit() 我想调用 st_crs_unit(cities) 并且该函数返回单位 "degrees" 或类似的。

st_crs() 以众所周知的文本格式生成有关 CRS 的信息，包括坐标系 (CS[]) 对两个轴使用 ANGLEUNIT "degree"，但此文本的结构差异很大跨坐标系统，所以我不能确定在某些系统上训练的正则表达式是否适用于所有人。

st_crs(cities)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
#>   wkt:
#> GEOGCRS["WGS 84",
#>     DATUM["World Geodetic System 1984",
#>         ELLIPSOID["WGS 84",6378137,298.257223563,
#>             LENGTHUNIT["metre",1]]],
#>     PRIMEM["Greenwich",0,
#>         ANGLEUNIT["degree",0.0174532925199433]],
#>     CS[ellipsoidal,2],
#>         AXIS["geodetic latitude (Lat)",north,
#>             ORDER[1],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         AXIS["geodetic longitude (Lon)",east,
#>             ORDER[2],
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>     USAGE[
#>         SCOPE["Horizontal component of 3D system."],
#>         AREA["World."],
#>         BBOX[-90,-180,90,180]],
#>     ID["EPSG",4326]]

^{由reprex package (v2.0.1) 于 2021 年 12 月 21 日创建}

例如，如果我们将相同的数据转换为使用 UTM 区域 30N 坐标系，st_crs() 的输出会发生很大变化。

st_crs(st_transform(cities, crs = 32630))
#> Coordinate Reference System:
#>   User input: EPSG:32630 
#>   wkt:
#> PROJCRS["WGS 84 / UTM zone 30N",
#>     BASEGEOGCRS["WGS 84",
#>         DATUM["World Geodetic System 1984",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",4326]],
#>     CONVERSION["UTM zone 30N",
#>         METHOD["Transverse Mercator",
#>             ID["EPSG",9807]],
#>         PARAMETER["Latitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8801]],
#>         PARAMETER["Longitude of natural origin",-3,
#>             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",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]],
#>     USAGE[
#>         SCOPE["Engineering survey, topographic mapping."],
#>         AREA["Between 6°W and 0°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Burkina Faso. Côte' Ivoire (Ivory Coast). Faroe Islands - offshore. France. Ghana. Gibraltar. Ireland - offshore Irish Sea. Mali. Mauritania. Morocco. Spain. United Kingdom (UK)."],
#>         BBOX[0,-6,84,0]],
#>     ID["EPSG",32630]]

^{由reprex package (v2.0.1) 于 2021 年 12 月 21 日创建}

是否存在返回空间对象测量单位的现有 R 函数？

【问题讨论】：

发布一些可重复的数据以及预期的输出，以便 ppl 可以帮助您

@ChrisRuehlemann 通常我会，但在这种情况下，我看不出它有什么帮助（存在从 SF 对象中提取测量单位的函数，或者不存在），它也有风险通过延长问题使问题变得不那么清楚。也就是说，我现在添加了一些虚拟代码，以防它对任何人有所帮助。

【参考方案1】：

st_crs() 有一个parameters 参数，当TRUE 时返回有用的CRS 参数列表，包括CRS 的单位。这是一个带有内置 nc 数据的示例：

library(sf)

nc_4267 <- read_sf(system.file("shape/nc.shp", package="sf"))
nc_3857 <- st_transform(nc_4267, 3857)

st_crs(nc_4267, parameters = TRUE)$units_gdal
#> [1] "degree"
st_crs(nc_3857, parameters = TRUE)$units_gdal
#> [1] "metre"

请注意，对于某些用途，st_is_longlat() 可能就足够了：

st_is_longlat(nc_4267)
#> [1] TRUE
st_is_longlat(nc_3857)
#> [1] FALSE

【讨论】：

太好了，谢谢！