Dear list members,
I hope this message finds you well.
Here's a possible problem that I won't be able to solve on my own!
When trying to determine the perimeter of a plot of land based on the ellipsoid, I've just noticed that the sf::st_perimeter() function doesn't seem to return the right value.
Actually, the sf::st_perimeter() function returns a different perimeter value from that obtained using other approaches, such as :
- with R: using sf::st_boundary() |> sf::st_length()
or using sf::st_cast() |> sf::st_length()
- with POSTGIS: ROUND(ST_Perimeter(ST_Transform(geom,4326)::geography)::numeric,4)
- with QGIS: round($perimeter,4)
In my particular case, sf::st perimeter() returns the value 129.0982 [m], whereas all the other approaches mentioned above return the value 129.2854 [m].
Of course the difference is very small, but I guess the sf::st_perimeter() function is supposed to return exactly the same value as that return by the other approaches.
Please find below a small REPREX to quickly test the three approaches in R (i.e. sf::st_perimeter() vs. sf::st_boundary() vs. sf::st_cast()).
Thank you in advance for your help
Cheers,
Loïc
----------------------------
REPREX
# Generate the land plot as an sf object (french projection: EPSG:2154)
plot <- 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)))
# What the plot looks like
plot(plot)
# The geom of the plot seems O.K.:
sf::st_is_valid(plot)
sf::sf_use_s2(FALSE)
# Calculating ellipsoidal perimeter with st_perimeter() --> 129.0982 [m]
plot |>
sf::st_transform(4326) |>
sf::st_perimeter()
# Calculating ellipsoidal perimeter with st_boundary() --> 129.2854 [m]
plot |>
sf::st_transform(4326) |>
sf::st_boundary() |>
sf::st_length()
# Calculating ellipsoidal perimeter with st_cast() --> 129.2854 [m]
plot |>
sf::st_transform(4326) |>
sf::st_cast("MULTILINESTRING") |>
sf::st_length()