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.pebe...@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.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 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.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://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. > >>>>>> > >>>>>> 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 > > -- > 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
