6

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?

2
  • 1
    Post some reproducible data along with an expected output so ppl can help you Dec 21, 2021 at 13:20
  • @ChrisRuehlemann usually I would, but in this case I don't see how it would help (either a function exists to extract the unit of measurement from an SF object or it doesn't) and it also risks making the question less clear by making it longer. That said, I've added some dummy code now in case it does help anyone.
    – Matt Ashby
    Dec 21, 2021 at 14:33

1 Answer 1

8

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
1
  • That's perfect, thank you!
    – Matt Ashby
    Dec 21, 2021 at 19:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.