https://stat.ethz.ch/pipermail/r-sig-geo/2025-May/029541.html
Below we calculate individual distances along polygon edges, and sum up. This works here only because this is a single-ring polygon.
The numbers in play are
129.248935: the planar distance in the crs
129.285384: the "true" distance, to some insanely accurate level
129.09815: the spherical distance according to S2 (reproduced below by setting radius and flattening)
WIP, I'd love to write this up and make it generally easy to explore in detail.
p <- structure(list(GEOM = structure(list(structure(list(structure(c(320102.67,
320111.85, 320131.82, 320134.33, 320138.5, 320141.93, 320135.33,
320118.58, 320113.45, 320105.18, 320102.67, 6845954.56, 6845962.15,
6845971.55, 6845971.9, 6845953.21, 6845936.12, 6845935.07, 6845931.54,
6845929.91, 6845950.57, 6845954.56), dim = c(11L, 2L))), class = c("XY",
"POLYGON", "sfg"))), n_empty = 0L, class = c("sfc_POLYGON", "sfc"
), precision = 0, bbox = structure(c(xmin = 320102.67, ymin = 6845929.91,
xmax = 320141.93, ymax = 6845971.9), class = "bbox"), crs = structure(list(
input = "EPSG:2154", wkt = "PROJCRS[\"RGF93 v1 / Lambert-93\",\n BASEGEOGCRS[\"RGF93 v1\",\n DATUM[\"Reseau Geodesique Francais 1993 v1\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4171]],\n CONVERSION[\"Lambert-93\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",46.5,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",3,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",44,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",700000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",6600000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting (X)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"France - onshore and offshore, mainland and Corsica (France métropolitaine including Corsica).\"],\n BBOX[41.15,-9.86,51.56,10.38]],\n ID[\"EPSG\",2154]]"), class = "crs"))), row.names = 1L, class = c("sf",
"data.frame"), sf_column = "GEOM", agr = structure(integer(0), levels = c("constant",
"aggregate", "identity"), class = "factor", names = character(0)))
options(digits = 9)
sf::sf_use_s2(TRUE)
pll <- sf::st_transform(p, "EPSG:4326")
sf::st_perimeter(pll)
#> 129.09815 [m]
sf::st_length(sf::st_boundary(pll))
#> 129.09815 [m]
sf::st_length(sf::st_cast(pll, "MULTILINESTRING"))
#> 129.09815 [m]
## that number is "s2 spherical", which we can reproduce with generic methods, get the coordinates for manual usage
xy <- as.matrix(wk::wk_coords(p)[, c("x", "y")])
ll <- terra::project(xy, to = "EPSG:4326", from = terra::crs(terra::vect(p)))
sum(geosphere::distMeeus(ll[-nrow(ll), ], ll[-1, ], a = 6371010, f = 0)) #https://github.com/r-spatial/s2/blob/342bbc23ff39d1be1c23537e8288f47b5ffe7c9d/src/s2/s2earth.h#L117
#> [1] 129.09815
sf::sf_use_s2(F)
#> Spherical geometry (s2) switched off
sf::st_perimeter(p)
#> Linking to GEOS 3.12.2, GDAL 3.10.3, PROJ 9.4.1; sf_use_s2() is FALSE
#> 129.248935 [m]
sf::st_length(sf::st_boundary(p))
#> 129.248935 [m]
sf::st_length(sf::st_cast(p, "MULTILINESTRING"))
#> 129.248935 [m]
sf::st_perimeter(pll) ## this appears to be the issue, as per the OP
#> 129.09815 [m]
sf::st_length(sf::st_boundary(pll))
#> 129.285384 [m]
sf::st_length(sf::st_cast(pll, "MULTILINESTRING"))
#> 129.285384 [m]
p0 <- sf::st_set_crs(p, NA)
sf::st_perimeter(p0)
#> [1] 129.248935
sf::st_length(sf::st_boundary(p0))
#> [1] 129.248935
sf::st_length(sf::st_cast(p0, "MULTILINESTRING"))
#> [1] 129.248935
## terra is Karney-enabled, for all crs (maybe not everywhere, but mostly, I think)
terra::perim(terra::vect(pll))
#> [1] 129.285384
terra::perim(terra::vect(p))
#> [1] 129.248935
terra::perim(terra::vect(sf::st_set_crs(p, NA)))
#> Warning: [perim] unknown CRS. Results can be wrong
#> [1] 129.248935
## decompose and sum up manually
idx <- seq(1, nrow(xy) - 1)
vlen <- function(x) sqrt(x[,1] ^2 + x[,2] ^ 2)
sum(vlen(xy[idx, ] - xy[idx + 1, ]))
#> [1] 129.248935
sum(geodist::geodist(terra::project(xy, to = "EPSG:4326", from = terra::crs(terra::vect(p))), sequential = TRUE, measure = "geodesic"))
#> object has no named columns; assuming order is lon then lat
#> [1] 129.285384
ll <- terra::project(xy, to = "EPSG:4326", from = terra::crs(terra::vect(p)))
geosphere::distGeo(ll[1, ], ll[2, ])
#> [1] 11.9147186
## what you get from Karney is what you would get from an individual equal distant projection on each pair
## function to local project to aeqd with one of the points on the centre
fun <- function(x, source) {
x <- terra::project(x, to = "EPSG:4326", from = source)
pr <- sprintf("+proj=aeqd +lon_0=%f +lat_0=%f", x[1, 1], x[1, 2])
x <- terra::project(x, to = pr, from = "EPSG:4326")
sum(vlen(x[1, , drop = FALSE] - x[2, ]))
}
d <- numeric(0)
for (i in 1:(nrow(xy)-1)) {
d <- c(d, fun(xy[i:(i+1), ], terra::crs(terra::vect(p))))
}
sum(d)
#> [1] 129.285384Created on 2025-05-17 with reprex v2.1.1
Maybe it's a lwgeom problem.
Do it!