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

[Edzer Pebesma]

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.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

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081