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

Tue Jan 31 15:50:29 CET 2023

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["Oblique Stereographic",
             ID["EPSG",9809]],
         PARAMETER["Latitude of natural origin",52.1561605555556,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8801]],
         PARAMETER["Longitude of natural origin",5.38763888888889,
             ANGLEUNIT["degree",0.0174532925199433],
             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 be 
suitable when the output is never sent to others.

Roger

>
> 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://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 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://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 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://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 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://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 list
>>>>>>> R-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 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://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 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 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 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