[R-sig-Geo] Re-scale units of coordinates in sf geometry?
Roger Bivand
Roger@B|v@nd @end|ng |rom nhh@no
Tue Jan 31 20:37:54 CET 2023
On Tue, 31 Jan 2023, Steve Gutreuter wrote:
> 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?
Sorry, there is no easy way both to deviate from an authorised CRS, such
as "EPSG:28992", and still be conformant, unless the authoriser has
covered that case. Only three +datum= values still exist in PROJ strings
for PROJ > 5, WGS84, NAD83 and NAD27. There are no +towgs84= values at
all. All the values are stored in proj.db, as are candidate transformation
pipelines. The same applies to current rgdal, the underlying PROJ and GDAL
code is the same in sf, terra and rgdal.
This operation is a conversion rather than a transformation, and is
lossless, but it isn't possible to approach through the legacy PROJ
string. It might be possible to edit a PROJJSON representation of
EPSG:28992 to convert the units and the false easting and northing values,
but the output will not be EPSG:28992, so the "id" component would have to
be deleted, see projinfo -o PROJJSON "EPSG:28992" if you also have PROJ
applications installed. The idea would be to create a valid PROJJSON
string and use it as the crs= argument to sf::st_transform().
For too much detail on CRS etc., possibly look at
https://www.youtube.com/watch?v=o0yMeb_UdE4&list=PLzREt6r1NenmWEidssmLm-VO_YmAh4pq9&index=2
or just the slides
https://rsbivand.github.io/csds_jan23/csds_crs_workshop_230119.html, but
note that the presented version in that talk is in error at 74:27 as
described in
https://rsbivand.github.io/csds_jan23/csds_crs_workshop_230131.html - the
latter should definitely be preferred.
Roger
> 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
>
--
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