[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 15:35:57 CET 2023
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
<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://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 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://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 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://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>>>>>
> >>>>>
> >>>>> [[alternative HTML version deleted]]
> >>>>>
> >>>>> _______________________________________________
> >>>>> R-sig-Geo mailing list
> >>>>> R-sig-Geo using 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 using 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
>
> --
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081
More information about the R-sig-Geo
mailing list