Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
On Tue, 31 Jan 2023, Steve Gutreuter wrote:
The need to re-scale the measurement units of the coordinates may not
be commonplace, but it is also not unreasonable. So, given that R rgdal
is to be retired in 2023, what is the "correct" way to preserve the new
units and the datum?
Sorry, there is no easy way both to deviate from an authorised CRS, such
as "EPSG:28992", and still be conformant, unless the authoriser has
covered that case. Only three +datum= values still exist in PROJ strings
for PROJ > 5, WGS84, NAD83 and NAD27. There are no +towgs84= values at
all. All the values are stored in proj.db, as are candidate transformation
pipelines. The same applies to current rgdal, the underlying PROJ and GDAL
code is the same in sf, terra and rgdal.
This operation is a conversion rather than a transformation, and is
lossless, but it isn't possible to approach through the legacy PROJ
string. It might be possible to edit a PROJJSON representation of
EPSG:28992 to convert the units and the false easting and northing values,
but the output will not be EPSG:28992, so the "id" component would have to
be deleted, see projinfo -o PROJJSON "EPSG:28992" if you also have PROJ
applications installed. The idea would be to create a valid PROJJSON
string and use it as the crs= argument to sf::st_transform().
For too much detail on CRS etc., possibly look at
https://www.youtube.com/watch?v=o0yMeb_UdE4&list=PLzREt6r1NenmWEidssmLm-VO_YmAh4pq9&index=2
or just the slides
https://rsbivand.github.io/csds_jan23/csds_crs_workshop_230119.html, but
note that the presented version in that talk is in error at 74:27 as
described in
https://rsbivand.github.io/csds_jan23/csds_crs_workshop_230131.html - the
latter should definitely be preferred.
Roger
Thanks to all,Steve
On Tue, 2023-01-31 at 18:36 +0100, Andrea Gilardi wrote:
Thanks, Roger. That's not a common problem (at least for me), but
itis typically useful to adjust the unit of measurement when I need
tointegrate sf objects with other packages (e.g. spatstat) where
thenumerical algorithm might benefit from rescaling.
Andrea
Il giorno mar 31 gen 2023 alle ore 15:50 Roger Bivand<
roger.biv...@nhh.no> ha scritto:
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.156160556 +lon_0=5.387639+k=0.079
+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.156160556, ANGLEUNIT["degree",0.01745329
25199433], ID["EPSG",8801]], PARAMETER["Longitu
de of natural
origin",5.387639, ANGLEUNIT["degree",0.01745329
25199433], ID["EPSG",8802]], PARAMETER["Scale
factor at natural
origin",0.079, 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
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
The need to re-scale the measurement units of the coordinates may not
be commonplace, but it is also not unreasonable. So, given that R rgdal
is to be retired in 2023, what is the "correct" way to preserve the new
units and the datum?
Thanks to all,Steve
On Tue, 2023-01-31 at 18:36 +0100, Andrea Gilardi wrote:
> Thanks, Roger. That's not a common problem (at least for me), but
> itis typically useful to adjust the unit of measurement when I need
> tointegrate sf objects with other packages (e.g. spatstat) where
> thenumerical algorithm might benefit from rescaling.
> Andrea
> Il giorno mar 31 gen 2023 alle ore 15:50 Roger Bivand<
> roger.biv...@nhh.no> ha scritto:
> > 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.156160556 +lon_0=5.387639+k=0.079
> > +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.156160556, ANGLEUNIT["degree",0.01745329
> > 25199433], ID["EPSG",8801]], PARAMETER["Longitu
> > de of natural
> > origin",5.387639, ANGLEUNIT["degree",0.01745329
> > 25199433], ID["EPSG",8802]], PARAMETER["Scale
> > factor at natural
> > origin",0.079, 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)#>xminyminxmaxy
> > > > > max#> 178.605 329.714 181.390 333.611
> > > > > I think this might be ea
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
Thanks, Roger. That's not a common problem (at least for me), but it
is typically useful to adjust the unit of measurement when I need to
integrate sf objects with other packages (e.g. spatstat) where the
numerical algorithm might benefit from rescaling.
Andrea
Il giorno mar 31 gen 2023 alle ore 15:50 Roger Bivand
ha scritto:
>
> 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.156160556 +lon_0=5.387639
> +k=0.079 +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.156160556,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8801]],
> PARAMETER["Longitude of natural origin",5.387639,
> ANGLEUNIT["degree",0.0174532925199433],
> ID["EPSG",8802]],
> PARAMETER["Scale factor at natural origin",0.079,
> 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
> > 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)
> >>> #>xminyminxmaxymax
> >>> #> 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 rig
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
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.156160556 +lon_0=5.387639
+k=0.079 +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.156160556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",5.387639,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.079,
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
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)
#>xminyminxmaxymax
#> 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
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(150, 500)), 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
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
Thanks again to both of you for the suggestion and the explanation.
Kind regards
Andrea
Il giorno mar 31 gen 2023 alle ore 15:28 Edzer Pebesma
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://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)
> > #>xminyminxmaxymax
> > #> 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
> > 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(150, 500)), 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.156160556 +lon_0=5.387639
> >> +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
> >> st_as_sf(meuse) |>
> >> st_transform("+proj=sterea +lat_0=52.156160556
> >> +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
> >> +ellps=bessel +units=km +no_defs") |>
> >> st_bbox()
> >> #xminyminxmaxymax
> >> # 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
> >>> 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
> > 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 app
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
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://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)
#>xminyminxmaxymax
#> 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
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(150, 500)), 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.156160556 +lon_0=5.387639
+k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
st_as_sf(meuse) |>
st_transform("+proj=sterea +lat_0=52.156160556
+lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
+ellps=bessel +units=km +no_defs") |>
st_bbox()
#xminyminxmaxymax
# 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
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
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]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
-
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
While the package interface isn't the friendliest, I think the crsuggest
package may be of interest for you. It's helpful for finding an appropriate
CRS in the units you need given an input CRS.
https://github.com/walkerke/crsuggest
I've used this to find an appropriate CRS and then transform based on the
top suggestion.
On Tue, Jan 31, 2023 at 8:48 AM 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://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)
> #>xminyminxmaxymax
> #> 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
> 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(150, 500)), 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.156160556 +lon_0=5.387639
> > +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
> > st_as_sf(meuse) |>
> >st_transform("+proj=sterea +lat_0=52.156160556
> > +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
> > +ellps=bessel +units=km +no_defs") |>
> >st_bbox()
> > #xminyminxmaxymax
> > # 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
> > > 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 <
> sgutreu...@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.
> >
> >
Re: [R-sig-Geo] 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)
#>xminyminxmaxymax
#> 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
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(150, 500)), 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.156160556 +lon_0=5.387639
> +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
> st_as_sf(meuse) |>
>st_transform("+proj=sterea +lat_0=52.156160556
> +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
> +ellps=bessel +units=km +no_defs") |>
>st_bbox()
> #xminyminxmaxymax
> # 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
> > 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
> >>> 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]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> >>>
> >>>[[alternative HTML version deleted]]
> >>>
> >>> ___
> >>> R-sig-Geo mailing list
> >>> R-sig-Geo@r-project.org
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >> --
> >> Edzer Pebesma
> >> Institute for Geoinformatics
> >> Heisenbergstrasse 2, 48151 Muenster, Germany
> >> Phone: +49 251
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
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(150, 500)), 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.156160556 +lon_0=5.387639
+k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
st_as_sf(meuse) |>
st_transform("+proj=sterea +lat_0=52.156160556
+lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
+ellps=bessel +units=km +no_defs") |>
st_bbox()
#xminyminxmaxymax
# 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
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
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]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
> 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(150, 500)), 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? Thanks Andrea Il giorno mer 25 gen 2023 alle ore 18:58 Edzer Pebesma 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 > > 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]] > >> > >> ___ > >> R-sig-Geo mailing list > >> R-sig-Geo@r-project.org > >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > >> > > > > [[alternative HTML version deleted]] > > > > ___ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > -- > Edzer Pebesma > Institute for Geoinformatics > Heisenbergstrasse 2, 48151 Muenster, Germany > Phone: +49 251 8333081 > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
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 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]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
I guess the question is: "safe for what?" - if numerical errors are a problem without rescaling, then you should do it. You might be able to choose the rescaling factor such that these errors are minimized. On 25/01/2023 18:37, Steve Gutreuter 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! -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
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 On Wed, Jan 25, 2023 at 12:37 PM Steve Gutreuter 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]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
