Problems with the "bru_fill_missing" command #200

gpanunzi · 2023-08-28T09:38:46Z

gpanunzi
Aug 28, 2023

Good morning everyone, after the latest updates to INLA and INLABRU, the codes that I was previously using are no longer functioning. The issue specifically pertains to the calculation of covariates. Prior to including them in the model, I used the following code:

Create mesh projection --------------------------------------------------

xx = range(mesh_medit$loc[,1])
yy = range(mesh_medit$loc[,2])
new_xy = as.matrix(expand.grid(seq(xx[1], xx[2], 1),
seq(yy[1], yy[2], 1)))

A = inla.spde.make.A(mesh = mesh_medit, loc = new_xy)

ips <- inlabru::ipoints(domain = mesh_medit)

################COVARIATES##############
########################################

load("data/BASEcovariates.RData")

FIXED COVARIATES

Depth -------------------------------------------------------------------

depth_raster <- rasterFromXYZ(depthFinal)
values = values(depth_raster)
depth_SPDF_s = SpatialPixelsDataFrame(coordinates(depth_raster),
data = data.frame(values),
proj4string = CRS(proj_utm_km))

values <- eval_spatial(data = depth_SPDF_s,
where = ips)

values2 <- inlabru:::bru_fill_missing(data = depth_SPDF_s,
where = ips,
values = values)

depth_SPDF_s = SpatialPixelsDataFrame(points = new_xy,
data = data.frame(depth_s = as.vector(A %*% values2)),
proj4string = inla.sp_get_crs(boundary))

Now I'm encountering the following error:
Error in sf::st_distance(data_coord[data_ok, , drop = FALSE], where[subset, :
st_crs(x) == st_crs(y) is not TRUE

I'm using version 2.8.0 of INLABRU; everything was working with version 2.7.0.

Answered by finnlindgren

Aug 28, 2023

Thanks. The last two show that while the "user input" parts of the crs of the converted objects are the same, they used different methods to construct the WKT (which is the actual information used) that generated equivalent but text-wise different representations. That's annoying...
I think one can work around this by transforming one of the objects into the version of the crs of the other object.
Please try this:

ips_sf <- sf::st_as_sf(ips)
covar_tr <- terra::rast(depth_SPDE_s)
ips_sf <- fm_transform(ips_sf, crs = fm_crs(covar_tr))
out <- bru_fill_missing(data = covar_tr, where = ips_sf, values = values)

Let me know if that works, and I can implement that workaround.

View full answer

finnlindgren · 2023-08-28T10:48:51Z

finnlindgren
Aug 28, 2023
Maintainer

Hi, I'm preparing to release version 2.9.0, so can you please check if the problem has been solved in the latest development version (2.8.0.9015)? Also, what INLA version did you use, and do you have fmesher installed? You will have fmesher if your INLA version is 23.08.18 or later; fmesher is a new package that takes over the role of the fm_* functions from inlabru and all the mesh handling from INLA, see https://inlabru-org.github.io/fmesher/articles/inla_conversion.html
The old inla.mesh.2d() functions will be around, but will call the new fmesher functions internally, so new code should transition to call fmesher directly. Similarly, ipoints is deprecated in favour of fmesher::fm_int() (that can generate both sf and sp output)

The bru_fill_missing problem likely comes from issue comes from inconsistencies in how CRS information is handled by sp, sf, and terra, which can make CRS objects appear different when having been converted between formats, even though it actually describes the same CRS. bru_fill_missing() now handles sp grids by converting them to terra format, and sp points by converting them to sf (this is part of the changes needed to safely handle the retirement of rgdal and rgeos, that will go away in October), so I would recommend converting your own data to sf and terra, since that would avoid the additional conversions inside inlabru, making it easier to troubleshoot these issues. But if there's a simple workaround the package could implement for sp inputs that would still be useful to better support existing code.

Can you show me the output from these commands in your current installation? proj4string(ips), proj4string(depth_SPDF_s), fm_CRS(ips), fm_CRS(depth_SPDE_s), fm_crs(ips), and fm_crs(depth_SPDE_s)

10 replies

finnlindgren Aug 28, 2023
Maintainer

I forgot some important ones as well: fm_crs(mesh_medit) and fm_CRS(mesh_medit) and mesh_medit$crs

finnlindgren Aug 28, 2023
Maintainer

Also, fm_crs(terra::rast(depth_SPDE_s)) and sf::st_crs(sf::st_as_sf(ips))

gpanunzi Aug 28, 2023
Author

Thank you very much for the prompt response. I am using INLA version 23.08.26 and fmesher version 0.1.2. I have updated inlabru to version 2.8.0.9015, but it still doesn't work. Here are the results of the requested information:

proj4string(ips)
[1] "+proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs"

proj4string(depth_SPDF_s)
[1] "+proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs"

fm_CRS(ips)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km
+no_defs
WKT2 2019 representation:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]]]]

fm_CRS(depth_SPDF_s)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km
+no_defs
WKT2 2019 representation:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]]]]

fm_crs(ips)
Coordinate Reference System:
User input: +proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs
wkt:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["Unknown based on WGS84 ellipsoid using towgs84=0,0,0,0,0,0,0",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Position Vector transformation (geog2D domain)",
ID["EPSG",9606]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]],
PARAMETER["X-axis rotation",0,
ID["EPSG",8608]],
PARAMETER["Y-axis rotation",0,
ID["EPSG",8609]],
PARAMETER["Z-axis rotation",0,
ID["EPSG",8610]],
PARAMETER["Scale difference",1,
ID["EPSG",8611]]]]

fm_crs(depth_SPDF_s)
Coordinate Reference System:
User input: +proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs
wkt:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["Unknown based on WGS84 ellipsoid using towgs84=0,0,0,0,0,0,0",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Position Vector transformation (geog2D domain)",
ID["EPSG",9606]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]],
PARAMETER["X-axis rotation",0,
ID["EPSG",8608]],
PARAMETER["Y-axis rotation",0,
ID["EPSG",8609]],
PARAMETER["Z-axis rotation",0,
ID["EPSG",8610]],
PARAMETER["Scale difference",1,
ID["EPSG",8611]]]]

fm_crs(mesh_medit)
Coordinate Reference System:
User input: +proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs
wkt:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]]]]

fm_CRS(mesh_medit)
Coordinate Reference System:
Deprecated Proj.4 representation:
+proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km
+no_defs
WKT2 2019 representation:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]]]]

mesh_medit$crs
Coordinate Reference System:
User input: +proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs
wkt:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]]]]

fm_crs(terra::rast(depth_SPDF_s))
Coordinate Reference System:
User input: +proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs
wkt:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["Unknown based on WGS84 ellipsoid using towgs84=0,0,0,0,0,0,0",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Position Vector transformation (geog2D domain)",
ID["EPSG",9606]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]],
PARAMETER["X-axis rotation",0,
ID["EPSG",8608]],
PARAMETER["Y-axis rotation",0,
ID["EPSG",8609]],
PARAMETER["Z-axis rotation",0,
ID["EPSG",8610]],
PARAMETER["Scale difference",1,
ID["EPSG",8611]]]]

sf::st_crs(sf::st_as_sf(ips))
Coordinate Reference System:
User input: +proj=utm +zone=33 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=km +no_defs
wkt:
BOUNDCRS[
SOURCECRS[
PROJCRS["unknown",
BASEGEOGCRS["unknown",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]]],
CONVERSION["UTM zone 33N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",15,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]],
ID["EPSG",16033]],
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]]]]],
TARGETCRS[
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
METHOD["Geocentric translations (geog2D domain)",
ID["EPSG",9603]],
PARAMETER["X-axis translation",0,
ID["EPSG",8605]],
PARAMETER["Y-axis translation",0,
ID["EPSG",8606]],
PARAMETER["Z-axis translation",0,
ID["EPSG",8607]]]]

finnlindgren Aug 28, 2023
Maintainer

Thanks. The last two show that while the "user input" parts of the crs of the converted objects are the same, they used different methods to construct the WKT (which is the actual information used) that generated equivalent but text-wise different representations. That's annoying...
I think one can work around this by transforming one of the objects into the version of the crs of the other object.
Please try this:

ips_sf <- sf::st_as_sf(ips)
covar_tr <- terra::rast(depth_SPDE_s)
ips_sf <- fm_transform(ips_sf, crs = fm_crs(covar_tr))
out <- bru_fill_missing(data = covar_tr, where = ips_sf, values = values)

Let me know if that works, and I can implement that workaround.

Answer selected by gpanunzi

gpanunzi Aug 28, 2023
Author

Yes, it's working now. Thank you so much!

finnlindgren Aug 28, 2023
Maintainer

I've implemented a potential workaround in the function itself now. Please try version 2.8.0.9016 with your original code and let me know if that also works.

gpanunzi Aug 28, 2023
Author

Yes, it's also working with the old code now.

finnlindgren Aug 28, 2023
Maintainer

That's great, thanks!

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Problems with the "bru_fill_missing" command #200

{{title}}

Replies: 1 comment 10 replies

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

Problems with the "bru_fill_missing" command #200

gpanunzi Aug 28, 2023

Create mesh projection --------------------------------------------------

FIXED COVARIATES

Depth -------------------------------------------------------------------

Replies: 1 comment · 10 replies

finnlindgren Aug 28, 2023 Maintainer

finnlindgren Aug 28, 2023 Maintainer

finnlindgren Aug 28, 2023 Maintainer

gpanunzi Aug 28, 2023 Author

finnlindgren Aug 28, 2023 Maintainer

gpanunzi Aug 28, 2023 Author

finnlindgren Aug 28, 2023 Maintainer

gpanunzi Aug 28, 2023 Author

finnlindgren Aug 28, 2023 Maintainer

gpanunzi
Aug 28, 2023

Replies: 1 comment 10 replies

finnlindgren
Aug 28, 2023
Maintainer

finnlindgren Aug 28, 2023
Maintainer

finnlindgren Aug 28, 2023
Maintainer

gpanunzi Aug 28, 2023
Author

finnlindgren Aug 28, 2023
Maintainer

gpanunzi Aug 28, 2023
Author

finnlindgren Aug 28, 2023
Maintainer

gpanunzi Aug 28, 2023
Author

finnlindgren Aug 28, 2023
Maintainer