[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
Josiah Parry
jo@|@h@p@rry @end|ng |rom gm@||@com
Sat May 17 00:23:48 CEST 2025
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]]
More information about the R-sig-Geo
mailing list