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

2023-01-31 Thread Roger Bivand

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=PLzREt6r1NenmWEidssmLm-VO_YmAh4pq9=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.biv...@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.156160556 +lon_0=5.387639+k=0.079
+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.156160556, ANGLEUNIT["degree",0.01745329
25199433], ID["EPSG",8801]], PARAMETER["Longitu
de of natural
origin",5.387639, ANGLEUNIT["degree",0.01745329
25199433], ID["EPSG",8802]], PARAMETER["Scale
factor at natural
origin",0.079, 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.pebe...@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=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=fVmKdcLGpNo7%2F%2Bn%2FEOiMs3qOX8lssLIh%2BW%2FLZkGjywk%3D=0
) and I foundthat the following code also seems to work
(although with a 

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

2023-01-31 Thread Steve Gutreuter
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.biv...@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.156160556 +lon_0=5.387639+k=0.079
> > +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.156160556, ANGLEUNIT["degree",0.01745329
> > 25199433], ID["EPSG",8801]], PARAMETER["Longitu
> > de of natural
> > origin",5.387639, ANGLEUNIT["degree",0.01745329
> > 25199433], ID["EPSG",8802]], PARAMETER["Scale
> > factor at natural
> > origin",0.079, 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.pebe...@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=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=fVmKdcLGpNo7%2F%2Bn%2FEOiMs3qOX8lssLIh%2BW%2FLZkGjywk%3D=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=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=rJdG7JgaY786%2BqOgeRTufTMdoGzmw97dvLb%2BUa6eYOc%3D=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)#>xminyminxmaxy
> > > > > max#> 178.605 329.714 181.390 333.611
> > > > > I think this might be easier than manually adjusting the
> > > 

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

2023-01-31 Thread Andrea Gilardi
Thanks, Roger. That's not a common problem (at least for me), but it
is typically useful to adjust the unit of measurement when I need to
integrate sf objects with other packages (e.g. spatstat) where the
numerical algorithm might benefit from rescaling.

Andrea

Il giorno mar 31 gen 2023 alle ore 15:50 Roger Bivand
 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.156160556 +lon_0=5.387639
> +k=0.079 +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.156160556,
>  ANGLEUNIT["degree",0.0174532925199433],
>  ID["EPSG",8801]],
>  PARAMETER["Longitude of natural origin",5.387639,
>  ANGLEUNIT["degree",0.0174532925199433],
>  ID["EPSG",8802]],
>  PARAMETER["Scale factor at natural origin",0.079,
>  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
> >  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=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=fVmKdcLGpNo7%2F%2Bn%2FEOiMs3qOX8lssLIh%2BW%2FLZkGjywk%3D=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=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=rJdG7JgaY786%2BqOgeRTufTMdoGzmw97dvLb%2BUa6eYOc%3D=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)
> >>> #>xminyminxmaxymax
> >>> #> 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 

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

2023-01-31 Thread Roger Bivand

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.156160556 +lon_0=5.387639 
+k=0.079 +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.156160556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",5.387639,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.079,
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
 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=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=fVmKdcLGpNo7%2F%2Bn%2FEOiMs3qOX8lssLIh%2BW%2FLZkGjywk%3D=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=05%7C01%7CRoger.Bivand%40nhh.no%7C9087c8564fa9442f824e08db0398f02f%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638107727639842599%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=rJdG7JgaY786%2BqOgeRTufTMdoGzmw97dvLb%2BUa6eYOc%3D=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)
#>xminyminxmaxymax
#> 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
 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(150, 500)), 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 

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

2023-01-31 Thread Andrea Gilardi
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
 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)
> > #>xminyminxmaxymax
> > #> 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
> >  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(150, 500)), 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.156160556 +lon_0=5.387639
> >> +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
> >> st_as_sf(meuse) |>
> >> st_transform("+proj=sterea +lat_0=52.156160556
> >> +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
> >> +ellps=bessel +units=km +no_defs") |>
> >> st_bbox()
> >> #xminyminxmaxymax
> >> # 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
> >>>  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 
> > 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 

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

2023-01-31 Thread 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)
#>xminyminxmaxymax
#> 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
 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(150, 500)), 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.156160556 +lon_0=5.387639
+k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
st_as_sf(meuse) |>
st_transform("+proj=sterea +lat_0=52.156160556
+lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
+ellps=bessel +units=km +no_defs") |>
st_bbox()
#xminyminxmaxymax
# 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
 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 
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



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

2023-01-31 Thread Josiah Parry
While the package interface isn't the friendliest, I think the crsuggest
package may be of interest for you. It's helpful for finding an appropriate
CRS in the units you need given an input CRS.

https://github.com/walkerke/crsuggest

I've used this to find an appropriate CRS and then transform based on the
top suggestion.

On Tue, Jan 31, 2023 at 8:48 AM 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)
> #>xminyminxmaxymax
> #> 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
>  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(150, 500)), 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.156160556 +lon_0=5.387639
> > +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
> > st_as_sf(meuse) |>
> >st_transform("+proj=sterea +lat_0=52.156160556
> > +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
> > +ellps=bessel +units=km +no_defs") |>
> >st_bbox()
> > #xminyminxmaxymax
> > # 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
> > >  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.
> > 
> >  

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

2023-01-31 Thread Andrea Gilardi
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)
#>xminyminxmaxymax
#> 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
 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(150, 500)), 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.156160556 +lon_0=5.387639
> +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"
> st_as_sf(meuse) |>
>st_transform("+proj=sterea +lat_0=52.156160556
> +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000
> +ellps=bessel +units=km +no_defs") |>
>st_bbox()
> #xminyminxmaxymax
> # 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
> >  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 
> >>> 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 

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

2023-01-31 Thread Edzer Pebesma




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(150, 500)), 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.156160556 +lon_0=5.387639 
+k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs"

st_as_sf(meuse) |>
  st_transform("+proj=sterea +lat_0=52.156160556 
+lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000 
+ellps=bessel +units=km +no_defs") |>

  st_bbox()
#xminyminxmaxymax
# 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
 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 
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

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


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

2023-01-31 Thread Andrea Gilardi
>  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(150, 500)), 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?

Thanks
Andrea


Il giorno mer 25 gen 2023 alle ore 18:58 Edzer Pebesma
 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 
> > 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

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


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

2023-01-25 Thread Edzer Pebesma




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


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

2023-01-25 Thread Edzer Pebesma
I guess the question is: "safe for what?" - if numerical errors are a 
problem without rescaling, then you should do it. You might be able to 
choose the rescaling factor such that these errors are minimized.


On 25/01/2023 18:37, Steve Gutreuter 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!


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


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

2023-01-25 Thread Josiah Parry
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

On Wed, Jan 25, 2023 at 12:37 PM Steve Gutreuter 
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