Dear Micha,
First of all, many thanks for exploring my request so quickly.
However, the results you obtained don't remove the doubt I have when trying to
determine the perimeter on the ellipsoid with the sf::st_perimeter() function.
There area two main reasons for this:
1 - When I focus on the last part of your e-mail (i.e. from where you
transform the crs from EPSG:2154 to EPSG:4326), I get the same results as you
do (namely, 129.0982 [m] whether I use the st_perimeter approach or the other
two approaches, i.e., sf::st_boundary() and sf::st_cast()), but if, and only
if, I leave the following default setting: sf::sf_use_s2(TRUE).
The problem is that if this setting is left as default, I understand that
the perimeter is calculated based on a sphere and not an ellipsoid (at least,
if I've correctly understood the following explanations here:
https://r-spatial.github.io/sf/reference/geos_measures.html).
To compute the perimeter on the ellipsoid, it seems that we need to use
sf::sf_use_s2(FALSE) before any calculation. And in doing so, I always end up
with the same values as those indicated in my initial e-mail, namely 129.0982
[m] with the sf::st_perimeter() function and 129.2854 [m] for the other two
approaches (i.e. sf::st_boundary() and sf::st_cast()).
2 - This value of 129.2854 [m] is, moreover, precisely the one found when the
calculation of the perimeter based on the ellipsoid is performed in QGIS
[$perimeter] or POSTGIS [ST_Perimeter(ST_Transform(geom,4326)::geography)].
This leads me to guess that it must be the right value and that,
consequently, there may be a problem with the sf::st_perimeter() function when
sf::sf_use_s2() is turned to FALSE.
Thank you very much again for your tests... and for your comment about the use
of "plot" in my REPREX (sorry for that)
For members who would also like to do the tests, please find below the REPREX
modified accordingly (using "plt" instead of "plot")
I have also added the session information at the bottom of this message.
Many thanks in advance.
Cheers,
Loïc
----------------------------
REPREX
# Generate the land plot as an sf object (french projection: EPSG:2154)
plt <- 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(plt)
# The geom of the plot seems O.K.:
sf::st_is_valid(plt)
# To compute ellipsoidal metrics
sf::sf_use_s2(FALSE)
# Calculating ellipsoidal perimeter with st_perimeter() --> 129.0982 [m]
plt |>
sf::st_transform(4326) |>
sf::st_perimeter()
# Calculating ellipsoidal perimeter with st_boundary() --> 129.2854 [m]
plt |>
sf::st_transform(4326) |>
sf::st_boundary() |>
sf::st_length()
# Calculating ellipsoidal perimeter with st_cast() --> 129.2854 [m]
plt |>
sf::st_transform(4326) |>
sf::st_cast("MULTILINESTRING") |>
sf::st_length()
--------------------------------------------------------------------------------------
> sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=French_France.utf8 LC_CTYPE=French_France.utf8
LC_MONETARY=French_France.utf8 LC_NUMERIC=C
LC_TIME=French_France.utf8
time zone: Europe/Paris
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sf_1.0-21
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 wk_0.9.4 sys_3.4.3
rstudioapi_0.17.1 jsonlite_2.0.0 happign_0.3.3
magrittr_2.0.3 rmarkdown_2.29 farver_2.1.2
[10] fs_1.6.6 vctrs_0.6.5 memoise_2.0.1
base64enc_0.1-3 terra_1.8-50 htmltools_0.5.8.1 usethis_3.1.0
s2_1.1.8 raster_3.6-32
[19] cellranger_1.1.0 pins_1.4.1 sass_0.4.10
shinyFiles_0.9.3 KernSmooth_2.23-21 bslib_0.9.0
htmlwidgets_1.6.4 keyring_1.3.2 httr2_1.1.2
[28] stars_0.6-8 cachem_1.1.0 igraph_2.1.4
mime_0.13 lifecycle_1.0.4 pkgconfig_2.0.3 R6_2.6.1
fastmap_1.2.0 shiny_1.10.0
[37] connections_0.2.0 digest_0.6.37 colorspace_2.1-1
mapview_2.11.2 lobstr_1.1.2.9000 shinycssloaders_1.1.0 leafem_0.2.4
pkgload_1.4.0 nngeo_0.4.8
[46] crosstalk_1.2.1 lwgeom_0.2-14 abind_1.4-8
RPostgreSQL_0.7-8 compiler_4.3.1 proxy_0.4-27 remotes_2.5.0
bit64_4.6.0-1 backports_1.5.0
[55] viridis_0.6.5 DBI_1.2.3 pkgbuild_1.4.7
rappdirs_0.3.3 sessioninfo_1.2.3 leaflet_2.2.2 waiter_0.2.5
classInt_0.4-11 tools_4.3.1
[64] formattable_0.2.1 units_0.8-7 rpostgis_1.4.4
odbc_1.6.1 zip_2.3.3 httpuv_1.6.16 glue_1.8.0
satellite_1.0.5 promises_1.3.2
[73] grid_4.3.1 generics_0.1.4 gtable_0.3.6
tzdb_0.5.0 class_7.3-22 RPostgres_1.4.8 tidyr_1.3.1
data.table_1.17.2 hms_1.1.3
[82] sp_2.2-0 xml2_1.3.8 pillar_1.10.2
stringr_1.5.1 later_1.4.2 dplyr_1.1.4
lattice_0.21-8 bit_4.6.0 tidyselect_1.2.1
[91] miniUI_0.1.2 knitr_1.50 gridExtra_2.3
xfun_0.52 stats4_4.3.1 devtools_2.4.5 dm_1.0.11
stringi_1.8.7 pacman_0.5.1
[100] geojsonsf_2.0.3 evaluate_1.0.3 shinyWidgets_0.9.0
codetools_0.2-19 archive_1.1.12 tibble_3.2.1 cli_3.6.5
xtable_1.8-4 rscontract_0.1.2
[109] jquerylib_0.1.4 dichromat_2.0-0.1 Rcpp_1.0.14
readxl_1.4.5 dbplyr_2.5.0 png_0.1-8
parallel_4.3.1 yyjsonr_0.1.20 ellipsis_0.3.2
[118] ggplot2_3.5.2 readr_2.1.5 assertthat_0.2.1
blob_1.2.4 profvis_0.4.0 urlchecker_1.0.1
viridisLite_0.4.2 scales_1.4.0 e1071_1.7-16
[127] purrr_1.0.4 rlang_1.1.6 shinyjs_2.1.0
rgeos_0.6-4
________________________________________
De : Micha Silver <micha_silver+...@proton.me>
Envoyé : vendredi 16 mai 2025 18:14
À : Loïc Valéry <lval...@outlook.fr>
Cc : r-sig-geo <r-sig-geo@r-project.org>
Objet : Re: [R-sig-Geo] The sf::st_perimeter() function doesn't seem to work
properly when trying to determine the perimeter of an sf object on the ellipsoid
Interesting. Thanks for the MWE!
Obviously, when you transform to WGS84 you'll get different length/perimeter
values than the original projection.
Here's what I see (I used plt as variable name to avoid conflicts with the
'plot` func):
st_crs(plt)[1]
$input
[1] "EPSG:2154"
# Calculating ellipsoidal perimeter with st_perimeter()
plt |> sf::st_perimeter()
129.0982 [m]
# Calculating ellipsoidal perimeter with st_boundary()
plt |>
sf::st_boundary() |>
sf::st_length()
129.0982 [m]
# Calculating ellipsoidal perimeter with st_cast()
plt |>
sf::st_cast("MULTILINESTRING") |>
sf::st_length()
129.0982 [m]
# Looks OK
# Now transformed
plt <- st_transform(plt, 4326)
# Calculating ellipsoidal perimeter with st_perimeter()
plt |> sf::st_perimeter()
129.0982 [m]
# Calculating ellipsoidal perimeter with st_boundary()
plt |>
sf::st_boundary() |>
sf::st_length()
129.0982 [m]
# Calculating ellipsoidal perimeter with st_cast()
plt |>
sf::st_cast("MULTILINESTRING") |>
sf::st_length()
129.0982 [m]
# Also OK
sessionInfo()
R version 4.2.2 Patched (2022-11-10 r83330)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 12 (bookworm)
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.21.so
locale:
[1] LC_CTYPE=en_IL.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_IL.UTF-8 LC_COLLATE=en_IL.UTF-8
[5] LC_MONETARY=en_IL.UTF-8 LC_MESSAGES=en_IL.UTF-8
[7] LC_PAPER=en_IL.UTF-8 LC_NAME=en_IL.UTF-8
[9] LC_ADDRESS=en_IL.UTF-8 LC_TELEPHONE=en_IL.UTF-8
[11] LC_MEASUREMENT=en_IL.UTF-8 LC_IDENTIFICATION=en_IL.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sf_1.0-20 rkward_0.7.5 devtools_2.4.5 usethis_3.1.0
HTH,
Micha
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
