[R-sig-Geo] Re-scale units of coordinates in sf geometry?

[Josiah Parry]

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 <andrea.gilardi using unimib.it>
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
>
> _______________________________________________
> 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]]