Comparision of projections in predict.gstat - to strict?
Closed this issue · 1 comments
waternumbers commented
The comparison of the projection information in predict.gstat seems stricter then testing for identical projections with the sp of sf packages. I can believe there is an explanation but can't find it. Example with data attached but in outline:
With two objects of class sf called obs and grd
> st_crs(obs)
Coordinate Reference System:
User input: WGS 84
wkt:
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]]
> st_crs(grd)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
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]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
passes a comparison of projections
> st_crs(obs) == st_crs(grd) ## TRUE
[1] TRUE
and when converted to sp
> sp_obs <- as(obs,"Spatial")
> sp_grd <- as(grd,"Spatial")
> sp::identicalCRS(sp_obs,sp_grd) ## TRUE
[1] TRUE
Running a model fit fails due to different projects
> nm <- names(grd)[1]
> m <- idw(formula(paste0(nm," ~ 1")),obs,grd)
[1] "+proj=longlat +datum=WGS84 +no_defs"
[1] "+proj=longlat +datum=WGS84 +no_defs"
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim, :
var1 : data item in gstat object and newdata have different coordinate reference systems
In addition: Warning messages:
1: In proj4string(d$data) :
CRS object has comment, which is lost in output
2: In proj4string(newdata) :
CRS object has comment, which is lost in output
The test performed by predict.gstat (line 41 in R code) would seem to be equivalent to
> identical( sp_obs@proj4string, sp_grd@proj4string )
[1] FALSE
This fails since the wkt (which as shown above is different) is copied into the comment attribute of the proj4string attribute of the data
> attributes(sp_obs@proj4string)$comment <- attributes(sp_grd@proj4string)$comment
> identical( sp_obs@proj4string, sp_grd@proj4string )
[1] TRUE
Full example data and code
projection_test.zip
edzer commented
Thanks for the reprex!