On Tue, 31 Jan 2023, Andrea Gilardi wrote:
Thanks again to both of you for the suggestion and the explanation.
Using:
meuse1 <- st_as_sf(meuse)
meuse2 <- st_transform(meuse1, sub("units=m", "units=km",
st_crs(meuse1)$proj4string))
gives
st_crs(meuse2)
Coordinate Reference System:
User input: +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
wkt:
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["Unknown based on Bessel 1841 ellipsoid",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["unknown",
METHOD["Oblique Stereographic",
ID["EPSG",9809]],
PARAMETER["Latitude of natural origin",52.1561605555556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",5.38763888888889,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9999079,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",155,
LENGTHUNIT["kilometre",1000],
ID["EPSG",8806]],
PARAMETER["False northing",463,
LENGTHUNIT["kilometre",1000],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["kilometre",1000,
ID["EPSG",9036]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["kilometre",1000,
ID["EPSG",9036]]]]
which does preserve the units, but does not preserve the datum. So may be
suitable when the output is never sent to others.
Roger
Kind regards
Andrea
Il giorno mar 31 gen 2023 alle ore 15:28 Edzer Pebesma
<edzer.pebesma at uni-muenster.de> ha scritto:
Yes, the pipeline approach bypasses GDAL, and doesn't result in an
object with an appropriate CRS as a consequence.
On 31/01/2023 14:47, Andrea Gilardi wrote:
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://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fproj.org%2Foperations%2Fconversions%2Funitconvert.html&data=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=fVmKdcLGpNo7%2F%2Bn%2FEOiMs3qOX8lssLIh%2BW%2FLZkGjywk%3D&reserved=0) 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://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fcdn.proj.org%2F&data=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=rJdG7JgaY786%2BqOgeRTufTMdoGzmw97dvLb%2BUa6eYOc%3D&reserved=0"
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:
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]]