Re-scale units of coordinates in sf geometry?

Thank you very much for your example! I briefly checked the examples
reported in ?st_transform and the docs at the PROJ website around
"Unit conversion"
(https://proj.org/operations/conversions/unitconvert.html) and I found
that the following code also seems to work (although with a warning
message that I'm not sure I understand):

library(sp)
demo(meuse, echo = FALSE, ask = FALSE)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE
sf_proj_network(TRUE)
#> [1] "https://cdn.proj.org"
st_as_sf(meuse) |> st_bbox()
#>   xmin   ymin   xmax   ymax
#> 178605 329714 181390 333611
meuse2 <- st_as_sf(meuse) |>
   st_transform(pipeline = "+proj=pipeline +step +proj=unitconvert
+xy_in=m +xy_out=km")
#> Warning in st_transform.sfc(st_geometry(x), crs, ...): pipeline not found in
#> PROJ-suggested candidate transformations
st_bbox(meuse2)
#>    xmin    ymin    xmax    ymax
#> 178.605 329.714 181.390 333.611

I think this might be easier than manually adjusting the WKT
representation of the CRS, but I'm not sure that this is really a
recommended way to transform units since, for some reason, it does not
modify the CRS of the objects which makes any analysis very difficult:

all.equal(st_crs(meuse), st_crs(meuse2))
#> [1] TRUE
library(mapview)
mapview(meuse) # seems right
mapview(meuse2) # completely off

Andrea

Il giorno mar 31 gen 2023 alle ore 12:33 Edzer Pebesma
<edzer.pebesma at uni-muenster.de> ha scritto:

On 31/01/2023 12:12, Andrea Gilardi wrote:

   st_transform can do transformations and conversions between different CRSs, unit conversions is a special case of that.

Dear all, I'm sorry if this is a trivial or stupid question but I was
wondering if you could provide an example of using st_transform to
apply a unit conversion transformation. For example, if I define
something like

library(sf)
pt <- st_sfc(st_point(c(1500000, 5000000)), crs = 3003) # some place
in North Italy
st_crs(3003) # units are expressed in metres

can I use st_transform to convert the unit measurement of pt from
metres to km and preserve that information in the CRS?

It seems so, using proj4strings (not recommended, but simple):

library(sp)
demo(meuse, echo = FALSE, ask = FALSE)
library(sf)
# Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
st_as_sf(meuse) |> st_bbox()
#   xmin   ymin   xmax   ymax
# 178605 329714 181390 333611
st_crs(meuse)$proj4string
# [1] "+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
st_as_sf(meuse) |>
    st_transform("+proj=sterea +lat_0=52.1561605555556
+lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000
+ellps=bessel +units=km +no_defs") |>
    st_bbox()
#    xmin    ymin    xmax    ymax
# 178.605 329.714 181.390 333.611

The better way would be to modify WKT representations of CRSs.

Thanks
Andrea


Il giorno mer 25 gen 2023 alle ore 18:58 Edzer Pebesma
<edzer.pebesma at uni-muenster.de> ha scritto:

On 25/01/2023 18:42, Josiah Parry wrote:

I wonder if `units::set_units()` is a better fit for the job
https://r-quantities.github.io/units/articles/measurement_units_in_R.html

It could definitely tell you which number to choose when going from,
say, US feet to km. Coordinates in sf geometries are not stored as units
objects, unit info is encoded in the CRS. st_transform can do
transformations and conversions between different CRSs, unit conversions
is a special case of that.

On Wed, Jan 25, 2023 at 12:37 PM Steve Gutreuter <sgutreuter at gmail.com>
wrote:

Is it safe to re-scale sf geometry coordinates from meters to
kilometers using, for example:

sfobj$geometry <- sfobj$geometry / 1000

It seems to work, but I understand too little about spatial data to
know whether that practice is actually safe.  I am working with spatial
data in a continental-scale equidistant conic projection, and the units
are meters.  It seems that kilometers would be better suited stochastic
partial differential equation modeling on a finite-element mesh, but
maybe that is a moot point.

Any and all advice will be much appreciated.

Thanks!
--
Steve Gutreuter

           [[alternative HTML version deleted]]

        [[alternative HTML version deleted]]

--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081