[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
Michael Sumner
md@umner @end|ng |rom gm@||@com
Sat May 17 14:13:13 CEST 2025
s2 is actually less accurate, because s2 is spherical (terra always uses
ellipsoidal Karney when crs is set, and sf+longlat-s2 uses Karney, as does
geodist, but sf uses "the map projection" when crs is nonlonlat or NA,
because this is an appropriate coordinate system in which to calculate
distance (within this region)). It's a fun complex to grok, not helpful for
inexpert use sadly, or with a lack of deep knowledge of the inter-package
history. Fun stuff.
There are three numbers in play in this code (I've worked on illustrating
this in different ways over the years, trying to get attention on it to no
avail, this is a great example and here is some code as a start at
unpicking it:
The numbers in play are 129.285384, 129.248935, 129.09815
https://gist.github.com/mdsumner/5daedfcbee851ca31a7aa9006dbf74cf
Choose your tools very carefully, it's pretty treacherous space here.
Best, Mike
On Sat, May 17, 2025 at 8:29 AM 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
>
--
Michael Sumner
Research Software Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsumner using gmail.com
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list