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

Extracting the measurement unit (degrees, metres, etc.) from spatial data in R

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

我尝试使用 sf 包中的 st_crs() 函数,它 return 是众所周知的文本格式的坐标参考系统。但是,我正在努力确定从该知名文本中提取测量单位的正则表达式是否能够在广泛的坐标系中可靠地运行。

是否存在 return 空间对象测量单位的现有函数?

例如,以下代码生成一个使用 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) 并且那个函数 return 单位 "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 创建

是否存在 return 空间对象测量单位的现有 R 函数?

st_crs() 有一个 parameters 参数,当 TRUE 时 returns 有用的 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