[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

Sat May 17 12:59:00 CEST 2025

After another check, I see you are correct: the st_perimeter() function returns a different value **when S2 is turned off**

sf::st_is_valid(plt)
[1] TRUE

plt <- st_transform(plt, 4326)
sf::sf_use_s2(FALSE)
Spherical geometry (s2) switched off
plt |>  sf::st_perimeter()
129.0982 [m]
# Calculating ellipsoidal perimeter with st_boundary()
plt |>
  sf::st_boundary() |>
  sf::st_length()
129.2854 [m]
# Calculating ellipsoidal perimeter with st_cast() 
plt |>
  sf::st_cast("MULTILINESTRING") |>
  sf::st_length()
129.2854 [m]

I don't know for sure, but I would guess that st_boundary() and st_cast() break the polygon into individual line strings, and then calculate length of each edge, whereas st_perimeter takes the polygon area and calculates perimeter, without using S2 so assuming a cartesian CRS, (which would be distorted).

In the help for distance functions it says:
"if sf_use_s2() is FALSE, ellipsoidal distances are computed using st_geod_distance which uses function geod_inverse from GeographicLib (part of PROJ); see Karney, Charles FF, 2013, Algorithms for geodesics, Journal of Geodesy 87(1), 43–55"

I take that to mean that *length* calculations use the st_geod_distance(). But perimeter must work differently.

Best, Micha

--
Micha Silver
cell: +972-523-665918

On Saturday, 17 May 2025 at 01:29, Josiah Parry <josiah.parry using gmail.com> wrote:

> I believe that using s2 is calculating a more accurate* perimeter because
> it is used only with spherical geometries. It uses a spherical
> calculation because you’re using a spherical projection. Otherwise you’re
> calculating the perimeter of a geography in 2D ignoring the fact that the
> data are not themselves projected onto a 2D plane.
> 
> On Fri, May 16, 2025 at 15:18 Loïc Valéry lvalery using outlook.fr wrote:
> 
> > 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+DEV using proton.me
> > Envoyé : vendredi 16 mai 2025 18:14
> > À : Loïc Valéry lvalery using outlook.fr
> > Cc : r-sig-geo r-sig-geo using 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 using r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
> 
> [[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo