Converting a SpatialPointsDataFrame into km units from m, in R

Question

I read in a UTM shapefile to R using readOGR. I want to use the resulting SpatialPointsDataFrame to build a mesh in the INLA package, but I've been told that to help with numerical stability, I should convert the geometries to km from m.

So far, I've been using INLA:

grid2 <- readOGR(here::here("PNG_grid.shp"))
mycrs <- fm_sp_get_crs(grid2)
mycrs2 <- fm_crs_set_lengthunit(mycrs, "km")
grid3 <- inla.spTransform(grid2, crs0 = mycrs, crs2 = mycrs2)

...but get an error "Invalid target CRS for SpatialPointsDataFrame". Looking at mycrs2 it's pretty weird. The value in mycrs is "+proj=utm +zone=54 +south +ellps=WGS72 +units=m +no_defs", so I thought this might work:

grid3 <- inla.spTransform(grid2, crs0 = "+proj=utm +zone=54 +south +ellps=WGS72 +units=m +no_defs", 
         crs1 = "+proj=utm +zone=54 +south +ellps=WGS72 +units=km +no_defs")

...but I get the same error. Is there anything in sp or rgdal that would help me?

Answer 1

Doing the spTransform with sp package tools works for me:

pts = data.frame(x=seq(1,100000,len=50), y=seq(1,100000,len=50), z=1:50)
coordinates(pts)=~x+y
proj4string(pts) ="+proj=utm +zone=54 +south +ellps=WGS72 +units=m +no_defs"

I get some warning messages but unless you are after very precise results we can let them pass.

Then:

pts2 = spTransform(pts,"+proj=utm +zone=54 +south +ellps=WGS72 +units=km +no_defs")

does not error. Its coordinate are now scaled:

> coordinates(pts)[1:4,]
         x        y
1    1.000    1.000
2 2041.796 2041.796
3 4082.592 4082.592
4 6123.388 6123.388
> coordinates(pts2)[1:4,]
         x        y
1 0.001000 0.001000
2 2.041796 2.041796
3 4.082592 4.082592
4 6.123388 6.123388

If that works for you but the transform in INLA doesn't, then its INLA to blame. If that fails for you then your sp or GDAL/OGR/PROJ installation need looking into. Version numbers of everything would be useful at this point.