gstat vulnerable to forthcoming changes in sp and rgdal

[rsbivand]

The revdep for gstat with PROJ6/GDAL3, sp 1.3-3 (my fork), rgdal 1.5-1 (R-Forge) fails with:

* checking examples ... ERROR
Running examples in ‘gstat-Ex.R’ failed
The error most likely occurred in:

> ### Name: krigeST
> ### Title: Ordinary global Spatio-Temporal Kriging
> ### Aliases: krigeST krige,formula,ST-method krigeSTTg
> ### Keywords: models
> 
> ### ** Examples
> 
> library(sp)
> library(spacetime)
> sumMetricVgm <- vgmST("sumMetric",
+                       space=vgm( 4.4, "Lin", 196.6,  3),
+                       time =vgm( 2.2, "Lin",   1.1,  2),
+                       joint=vgm(34.6, "Exp", 136.6, 12),
+                       stAni=51.7)
> 
> data(air)
> 
> if (!exists("rural"))
+   rural = STFDF(stations, dates, data.frame(PM10 = as.vector(air)))
> 
> rr <- rural[,"2005-06-01/2005-06-03"]
> rr <- as(rr,"STSDF")
> 
> x1 <- seq(from=6,to=15,by=1)
> x2 <- seq(from=48,to=55,by=1)
> 
> DE_gridded <- SpatialPoints(cbind(rep(x1,length(x2)), rep(x2,each=length(x1))), 
+                             proj4string=CRS(proj4string(rr@sp)))
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum WGS_1984 in CRS definition,
 but +towgs84= values preserved
> gridded(DE_gridded) <- TRUE
> DE_pred <- STF(sp=as(DE_gridded,"SpatialPoints"), time=rr@time)
> DE_kriged <- krigeST(PM10~1, data=rr, newdata=DE_pred,
+                      modelList=sumMetricVgm)
Error in krigeST(PM10 ~ 1, data = rr, newdata = DE_pred, modelList = sumMetricVgm) : 
  identical(proj4string(data@sp), proj4string(newdata@sp)) is not TRUE
Calls: krigeST -> stopifnot
Execution halted

[edzer]

This seems to have done it. I now have (manually) set the projargs in the binary air.rda objects in spacetime to +proj=longlat +datum=WGS84, and reset them in the gstat manuals so that they are the actual output of PROJ for this, which depends on version, and then gets propagated/compared.

It would be nice to have a portable and not time varying encoding of WGS84 longlat CRS coding. I'm also thinking of GADM:

> proj4string(raster::getData('GADM', country='FRA', level=1))
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

[rsbivand]

Right, I found GADM when trying to find a widely used source of frozen sp and sf vector objects as rda/rds files. Maybe I need to revisit rgdal::spTransform() to enforce calling sp::CRS() on input crs without comments? Then we are protected but waste if() checks when things stabilise? However, I can't see how to insert the re-instantiated crs back into the objects as read.

[rsbivand]

A thought - I'll be teaching CRS stuff Tuesday 3/12, 09:15-11:00, and the last 45m on adaptations. I'll be (trying to) stream it. I could try to bring in Mike Sumner too (mid-evening for him), would talking about this make sense? Participants unsuspecting PhD course people whom I don't yet know. Hope to put streamed chunks on OpenGeoHub youtube when ready.