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

I would like to extract the unit of measurement (decimal degrees, metres, feet, etc.) from a spatial object in R. For example, if I have an SF data frame that uses the WGS84 co-ordinate reference system (EPSG:4326), I would like to be able to determine that the co-ordinates are specified in decimal degrees. Similarly, I'd like to be able to determine that UTM co-ordinates (e.g. EPSG:32615) are specified in metres.

I have tried using the st_crs() function from the sf package, which returns the co-ordinate reference system in well-known text format. However, I'm struggling to be certain that a regex that extracts the unit of measurement from that well-known text will operate reliably for a wide range of co-ordinate systems.

Is there an existing function that returns the measurement unit for a spatial object?

For example, the following code produces an SF data frame that uses the WGS84 co-ordinate system:

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)

^{Created on 2021-12-21 by the reprex package (v2.0.1)}

I am ideally looking for a function that allows me to determine that the spatial unit for this dataset is decimal degrees, e.g. if the function was called st_crs_unit() I would like to call st_crs_unit(cities) and that function return the unit "degrees" or similar.

st_crs() produces information about the CRS in well-known text format, including that the co-ordinate system (CS[]) uses the ANGLEUNIT "degree" for both axes, but the structure of this text varies considerably across co-ordinate systems, so I cannot be sure a regex trained on some systems will work for all.

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]]

^{Created on 2021-12-21 by the reprex package (v2.0.1)}

For example, if we transform the same data to use the UTM zone 30N co-ordinate system, the output from st_crs() changes substantially.

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]]

^{Created on 2021-12-21 by the reprex package (v2.0.1)}

Is there an existing R function that returns the measurement unit for a spatial object?

CodePudding user response：

st_crs() has a parameters argument that returns a list of useful CRS parameters when TRUE, including the units of the CRS. Here's an example with the built-in nc data:

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"

Note that for some purposes st_is_longlat() might be sufficient:

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