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

Tue Jan 31 20:00:59 CET 2023

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.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["Obli
> > que
> > Stereographic",             ID["EPSG",9809]],         PARAMETER["La
> > titude of natural
> > origin",52.1561605555556,             ANGLEUNIT["degree",0.01745329
> > 25199433],             ID["EPSG",8801]],         PARAMETER["Longitu
> > de of natural
> > origin",5.38763888888889,             ANGLEUNIT["degree",0.01745329
> > 25199433],             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 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.pebesma using 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)#>    xmin    ymin    xmax    y
> > > > > max#> 178.605 329.714 181.390 333.611
> > > > > I think this might be easier than manually adjusting the
> > > > > WKTrepresentation of the CRS, but I'm not sure that this is
> > > > > really arecommended way to transform units since, for some
> > > > > reason, it does notmodify the CRS of the objects which makes
> > > > > any analysis very difficult:
> > > > > all.equal(st_crs(meuse), st_crs(meuse2))#> [1]
> > > > > TRUElibrary(mapview)mapview(meuse) # seems
> > > > > rightmapview(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 waswondering if you could provide an
> > > > > > > example of using st_transform toapply a unit conversion
> > > > > > > transformation. For example, if I definesomething like
> > > > > > > library(sf)pt <- st_sfc(st_point(c(1500000, 5000000)),
> > > > > > > crs = 3003) # some placein North Italyst_crs(3003) #
> > > > > > > units are expressed in metres
> > > > > > > can I use st_transform to convert the unit measurement of
> > > > > > > pt frommetres 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 TRUEst_as_sf(meuse) |>
> > > > > > st_bbox()#   xmin   ymin   xmax   ymax# 178605 329714
> > > > > > 181390 333611st_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.
> > > > > > > ThanksAndrea
> > > > > > > 
> > > > > > > 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 unitsobjects, unit info is
> > > > > > > > encoded in the CRS. st_transform can dotransformations
> > > > > > > > and conversions between different CRSs, unit
> > > > > > > > conversionsis 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 tokilometers using, for example:
> > > > > > > > > > sfobj$geometry <- sfobj$geometry / 1000
> > > > > > > > > > It seems to work, but I understand too little about
> > > > > > > > > > spatial data toknow whether that practice is
> > > > > > > > > > actually safe.  I am working with spatialdata in a
> > > > > > > > > > continental-scale equidistant conic projection, and
> > > > > > > > > > the unitsare meters.  It seems that kilometers
> > > > > > > > > > would be better suited stochasticpartial
> > > > > > > > > > differential equation modeling on a finite-element
> > > > > > > > > > mesh, butmaybe that is a moot point.
> > > > > > > > > > Any and all advice will be much appreciated.
> > > > > > > > > > Thanks!--Steve Gutreuter
> > > > > > > > > >            [[alternative HTML version deleted]]
> > > > > > > > > > _______________________________________________R-
> > > > > > > > > > sig-Geo mailing listR-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 listR-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 PebesmaInstitute for
> > > > > > > > GeoinformaticsHeisenbergstrasse 2, 48151 Muenster,
> > > > > > > > GermanyPhone: +49 251 8333081
> > > > > > > > _______________________________________________R-sig-
> > > > > > > > Geo mailing listR-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 PebesmaInstitute for
> > > > > > GeoinformaticsHeisenbergstrasse 2, 48151 Muenster,
> > > > > > GermanyPhone: +49 251 8333081
> > > > 
> > > > --Edzer PebesmaInstitute for GeoinformaticsHeisenbergstrasse 2,
> > > > 48151 Muenster, GermanyPhone: +49 251 8333081
> > > 
> > > _______________________________________________R-sig-Geo mailing 
> > > listR-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 BivandEmeritus ProfessorDepartment 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
> 
> _______________________________________________R-sig-Geo mailing 
> listR-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

	[[alternative HTML version deleted]]