1

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?

1 Answer 1

2

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.

1
  • That works, thank you! I've opened an INLA issue on GitHub - there is also something slightly amiss with inlabru's fm_spTransform at the moment, which they're looking into. May 21, 2021 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.