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 <lval...@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+...@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 > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
