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["Obli
que
Stereographic", ID["EPSG",9809]], PARAMETER["La
titude of natural
origin",52.1561605555556, ANGLEUNIT["degree",0.01745329
25199433], ID["EPSG",8801]], PARAMETER["Longitu
de of natural
origin",5.38763888888889, ANGLEUNIT["degree",0.01745329
25199433], 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 besuitable when the output is never sent to others.
Roger
Kind regardsAndrea
Il giorno mar 31 gen 2023 alle ore 15:28 Edzer Pebesma<
edzer.pebe...@uni-muenster.de> ha scritto:
Yes, the pipeline approach bypasses GDAL, and doesn't result in
anobject 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
examplesreported 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 foundthat the following code also seems to work
(although with a warningmessage 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 TRUEsf_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 333611meuse2 <- 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
transformationsst_bbox(meuse2)#> xmin ymin xmax y
max#> 178.605 329.714 181.390 333.611
I think this might be easier than manually adjusting the
WKTrepresentation of the CRS, but I'm not sure that this is
really arecommended way to transform units since, for some
reason, it does notmodify the CRS of the objects which makes
any analysis very difficult:
all.equal(st_crs(meuse), st_crs(meuse2))#> [1]
TRUElibrary(mapview)mapview(meuse) # seems
rightmapview(meuse2) # completely off
Andrea
Il giorno mar 31 gen 2023 alle ore 12:33 Edzer Pebesma<
edzer.pebe...@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 waswondering if you could provide an
example of using st_transform toapply a unit conversion
transformation. For example, if I definesomething like
library(sf)pt <- st_sfc(st_point(c(1500000, 5000000)),
crs = 3003) # some placein North Italyst_crs(3003) #
units are expressed in metres
can I use st_transform to convert the unit measurement of
pt frommetres 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 TRUEst_as_sf(meuse) |>
st_bbox()# xmin ymin xmax ymax# 178605 329714
181390 333611st_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.
ThanksAndrea
Il giorno mer 25 gen 2023 alle ore 18:58 Edzer Pebesma<
edzer.pebe...@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://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fr-quantities.github.io%2Funits%2Farticles%2Fmeasurement_units_in_R.html&data=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=yBdfWvCeDu9W2Y9yXC5uRKA0QnzVdoZVWOr9KXfjK50%3D&reserved=0
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 unitsobjects, unit info is
encoded in the CRS. st_transform can dotransformations
and conversions between different CRSs, unit
conversionsis a special case of that.
On Wed, Jan 25, 2023 at 12:37 PM Steve Gutreuter <
sgutreu...@gmail.com>wrote:
Is it safe to re-scale sf geometry coordinates from
meters tokilometers using, for example:
sfobj$geometry <- sfobj$geometry / 1000
It seems to work, but I understand too little about
spatial data toknow whether that practice is
actually safe. I am working with spatialdata in a
continental-scale equidistant conic projection, and
the unitsare meters. It seems that kilometers
would be better suited stochasticpartial
differential equation modeling on a finite-element
mesh, butmaybe that is a moot point.
Any and all advice will be much appreciated.
Thanks!--Steve Gutreuter
[[alternative HTML version deleted]]
_______________________________________________R-
sig-Geo mailing listr-sig-...@r-project.org
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=krfLyCtoYG5Iw%2FnktA5LqG%2BfNIcaNXIWmWEDzkUfpK0%3D&reserved=0
[[alternative HTML version deleted]]
_______________________________________________R-sig-
Geo mailing listr-sig-...@r-project.org
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=krfLyCtoYG5Iw%2FnktA5LqG%2BfNIcaNXIWmWEDzkUfpK0%3D&reserved=0
--Edzer PebesmaInstitute for
GeoinformaticsHeisenbergstrasse 2, 48151 Muenster,
GermanyPhone: +49 251 8333081
_______________________________________________R-sig-
Geo mailing listr-sig-...@r-project.org
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=krfLyCtoYG5Iw%2FnktA5LqG%2BfNIcaNXIWmWEDzkUfpK0%3D&reserved=0
--Edzer PebesmaInstitute for
GeoinformaticsHeisenbergstrasse 2, 48151 Muenster,
GermanyPhone: +49 251 8333081
--Edzer PebesmaInstitute for GeoinformaticsHeisenbergstrasse 2,
48151 Muenster, GermanyPhone: +49 251 8333081
_______________________________________________R-sig-Geo mailing
listr-sig-...@r-project.org
https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geo&data=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=krfLyCtoYG5Iw%2FnktA5LqG%2BfNIcaNXIWmWEDzkUfpK0%3D&reserved=0
--Roger BivandEmeritus ProfessorDepartment of Economics, Norwegian
School of Economics,Postboks 3490 Ytre Sandviken, 5045 Bergen,
Norway.e-mail: roger.biv...@nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
