[R-sig-Geo] Re-scale units of coordinates in sf geometry?
Andrea Gilardi
@ndre@@g||@rd| @end|ng |rom un|m|b@|t
Tue Jan 31 18:36:18 CET 2023
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
<Roger.Bivand using 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.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 using 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 using 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 using 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 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 using 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]]
> >>>>>>>>
> >>>>>>>> _______________________________________________
> >>>>>>>> R-sig-Geo mailing list
> >>>>>>>> R-sig-Geo using 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 list
> >>>>>>> R-sig-Geo using 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 Pebesma
> >>>>>> Institute for Geoinformatics
> >>>>>> Heisenbergstrasse 2, 48151 Muenster, Germany
> >>>>>> Phone: +49 251 8333081
> >>>>>>
> >>>>>> _______________________________________________
> >>>>>> R-sig-Geo mailing list
> >>>>>> R-sig-Geo using 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 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
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo using 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 Bivand
> Emeritus Professor
> Department of Economics, Norwegian School of Economics,
> Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
> e-mail: Roger.Bivand using nhh.no
> https://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
More information about the R-sig-Geo
mailing list