Dear list members,
I hope this message finds you well.
Here's a possible problem that I won't be able to solve on my own!
When trying to determine the perimeter of a plot of land based on the ellipsoid, I've just noticed that the sf::st_perimeter() function doesn't seem to return the right value.
Actually, the sf::st_perimeter() function returns a different perimeter value from that obtained using other approaches, such as :
- with R: using sf::st_boundary() |> sf::st_length()
or using sf::st_cast() |> sf::st_length()
- with POSTGIS: ROUND(ST_Perimeter(ST_Transform(geom,4326)::geography)::numeric,4)
- with QGIS: round($perimeter,4)
In my particular case, sf::st perimeter() returns the value 129.0982 [m], whereas all the other approaches mentioned above return the value 129.2854 [m].
Of course the difference is very small, but I guess the sf::st_perimeter() function is supposed to return exactly the same value as that return by the other approaches.
Please find below a small REPREX to quickly test the three approaches in R (i.e. sf::st_perimeter() vs. sf::st_boundary() vs. sf::st_cast()).
Thank you in advance for your help
Cheers,
Lo?c
----------------------------
REPREX
# Generate the land plot as an sf object (french projection: EPSG:2154)
plot <- 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(plot)
# The geom of the plot seems O.K.:
sf::st_is_valid(plot)
sf::sf_use_s2(FALSE)
# Calculating ellipsoidal perimeter with st_perimeter() --> 129.0982 [m]
plot |>
sf::st_transform(4326) |>
sf::st_perimeter()
# Calculating ellipsoidal perimeter with st_boundary() --> 129.2854 [m]
plot |>
sf::st_transform(4326) |>
sf::st_boundary() |>
sf::st_length()
# Calculating ellipsoidal perimeter with st_cast() --> 129.2854 [m]
plot |>
sf::st_transform(4326) |>
sf::st_cast("MULTILINESTRING") |>
sf::st_length()
The sf::st_perimeter() function doesn't seem to work properly when trying to determine the perimeter of an sf object on the ellipsoid
7 messages · Micha Silver, Loïc Valéry, Josiah Parry +1 more
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
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 at proton.me>
Envoy? :?vendredi 16 mai 2025 18:14
? :?Lo?c Val?ry <lvalery at outlook.fr>
Cc?:?r-sig-geo <r-sig-geo at 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
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 at 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 at proton.me>
Envoy? : vendredi 16 mai 2025 18:14
? : Lo?c Val?ry <lvalery at outlook.fr>
Cc : r-sig-geo <r-sig-geo at 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 at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
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 at 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 at 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 at proton.me
Envoy? : vendredi 16 mai 2025 18:14
? : Lo?c Val?ry lvalery at outlook.fr
Cc : r-sig-geo r-sig-geo at 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 at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
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 at 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 at 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 at proton.me> Envoy? : vendredi 16 mai 2025 18:14 ? : Lo?c Val?ry <lvalery at outlook.fr> Cc : r-sig-geo <r-sig-geo at 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 at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at 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 at gmail.com [[alternative HTML version deleted]]
Dear Micha, Josiah and Mike, First of all, a very big thank you to all three of you for taking the time to share your thoughts on the problem I raised. Mike, an even bigger thank you for the detailed explanations. It's really appreciated.? The details of the calculations on your Github post are particularly interesting and confirm the doubt I was expressing about the sf::st_perimeter() function when sf::sf_use_s2() is set as FALSE (i.e. the case where the function returns 129.09815 [m] as shown in my reprex). I understand from your Github post that you validate the fact that there is an error somewhere in the function's code, which will most probably be corrected soon (by you or Edzer, I guess). Many thanks again for your particularly enlightening and helpful contribution... and thanks in advance to the person who will make the correction in the code. Cheers, Lo?c ________________________________________ De :?Michael Sumner <mdsumner at gmail.com> Envoy? :?samedi 17 mai 2025 14:13 ? :?Josiah Parry <josiah.parry at gmail.com> Cc?:?Lo?c Val?ry <lvalery at outlook.fr>; r-sig-geo <r-sig-geo at 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 ? 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 at 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 at 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 at proton.me>
Envoy? : vendredi 16 mai 2025 18:14
? : Lo?c Val?ry <lvalery at outlook.fr>
Cc : r-sig-geo <r-sig-geo at 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 at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
? ? ? ? [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo at 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 at gmail.com