r-spatial / gstat

Spatial and spatio-temporal geostatistical modelling, prediction and simulation
http://r-spatial.github.io/gstat/
GNU General Public License v2.0
195 stars 49 forks source link

Comparision of projections in predict.gstat - to strict? #95

Closed waternumbers closed 2 years ago

waternumbers commented 2 years ago

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 2 years ago

Thanks for the reprex!