Test whether coordinate reference system has units of meters

Question

Is there a simple way to test whether a spatial object's coordinate reference system expresses coordinates in units of meters? The snippet below shows my initial attempt at a function that provides such a test, but I don't have a great sense for how broadly reliable it is.

Is a better (more reliable and/or more widely tested) approach already available? Failing that, are there commonly used coordinate reference systems for which my approach is likely to fail?

library(sf)

## Function to test whether an sf object's projection expresses
## coordinates in units of meters by checking that the projection:
## (a) is not longlat; and (b) comes with no multiplier factor for
## converting units to meters.
is_proj_in_meters <- function(obj) {
    !is.na(st_crs(obj)) &
        !st_is_longlat(obj) &
         is.null(st_crs(obj)$to_meter)
}


## Create four sample polygons, each with a different CRS

g <- st_as_sfc("POLYGON ((-61.66957 10.69214, -61.565 10.75728, -61.37453 10.77654, -61.40721 10.60681, -61.66957 10.69214))")

## (1) No CRS supplied (units are unspecified)
a <- st_as_sf(data.frame(id = 1, geometry = g))
st_crs(a)$to_meter ## NA
st_is_longlat(a)   ## NA

## (2) WGS 84 projection (units are degrees)
b <- st_as_sf(data.frame(id = 1, geometry = g), crs = 4326)
st_crs(b)$to_meter ## NULL
st_is_longlat(b)   ## [1] TRUE

## (3) 1903 Trinidad grid (unit is Clarke's foot)
c <- st_transform(b, 2314)
st_crs(c)$to_meter ## [1] 0.30479726540000002
st_is_longlat(c)   ## [1] FALSE

## (4) UTM zones applied to Trinidad and Tobago (units are meters)
d <- st_transform(b, 32621)
st_crs(d)$to_meter ## NULL
st_is_longlat(d)   ## [1] FALSE


## Finally, test `is_proj_in_meters()`
sapply(list(a,b,c,d), is_proj_in_meters)
## [1] FALSE FALSE FALSE  TRUE

Answer 1

I cannot comment on your method, but out of interest I implemented a linearUnits method in terra that calls the GDAL method GetLinearUnits. Perhaps this is useful for you for comparison (although I assume that sf calls the same GDAL function). The units are also expressed in meters. I set the value for empty crs to NA and for lon/lat to zero. Note that GDAL assumes 1 when there is no data.

library(terra)
#terra version 0.9.14
x <- rast()
 
crs(x) <- ""
linearUnits(x)
#[1] NaN
 
crs(x) <- "+proj=longlat +datum=WGS84"
linearUnits(x)
#[1] 0
 
crs(x) <- "+proj=utm +zone=1 +units=cm"
linearUnits(x)
#[1] 0.01
 

With your example data

g <- vect("POLYGON ((-61.66957 10.69214, -61.565 10.75728, -61.37453 10.77654, -61.40721 10.60681, -61.66957 10.69214))")
 
linearUnits(g)
#[1] NaN
 
crs(g) <- "epsg:4326"
linearUnits(g)
#[1] 0
 
x <- project(g, "epsg:2314")
linearUnits(x)
#[1] 0.3047973
 
linearUnits(project(g, "epsg:32621"))
#[1] 1

I now use this in terra to compute areas and distances in meters, in stead of "in the units of the crs".

Answer 2

I suppose your method is reliable, as well as the GDAL getLinearUnits method. At the end they do the same thing by checking the units which are defined for the coordinate system. Compare

https://epsg.org/crs/wkt/id/2314

LENGTHUNIT["Clarke's foot",0.3047972654,ID["EPSG",9005]],

https://epsg.org/crs/wkt/id/32621

LENGTHUNIT["metre",1,ID["EPSG",9001]],

https://epsg.org/crs/wkt/id/4326

ANGLEUNIT["degree (supplier to define representation)",0.0174532925199433,ID["EPSG",9122]],

Because the last one is not projected it does not have lengthunit.

Answer 3

Try projection_units() from the metacapa package.