changes in 'raster::buffer' cause R CMD check failures in 'inRadius'

[achubaty]

Recent builds failing due to what looks like a change in the raster package that affects inRadius.

Sarah tracked down the issue to raster::buffer:

> w1 <- createWorld(minPxcor = 0, maxPxcor = 4, minPycor = 0, maxPycor = 4)
> agents <- patch(w1, 0, 0)
> agentsSP <- SpatialPoints(coords = agents)
> agentsSP
class       : SpatialPoints 
features    : 1 
extent      : 0, 0, 0, 0  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
> buffer(agentsSP, dissolve = FALSE, width = 2)
class       : SpatialPolygons 
features    : 1 
extent      : -1.796631e-05, 1.796631e-05, -1.808739e-05, 1.808739e-05  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
Warning messages:
1: In couldBeLonLat(...) : CRS is NA. Assuming it is longitude/latitude
2: In .local(x, ...) : crs unknown, assuming lonlat

It looks like there was a recent change to extent and buffer (see https://cran.r-project.org/web/packages/raster/ChangeLog) in the raster package.

[SarahBauduin]

When using buffer(), with the width = 2, I want (and had before the new raster package version) a SpatialPolygons of extent -2, 2, -2, 2. Now because there is no CRS in the SpatialPoints "agentsSP", it assumes it's lat/long and it's not working properly/as before. The way around I found is to add a projected CRS to the SpatialPoints:

proj4string(agentsSP) <-CRS("+proj=utm +zone=10 +datum=WGS84")
raster::buffer(agentsSP, dissolve = FALSE, width = 2)
class       : SpatialPolygons 
features    : 1 
extent      : -2, 2, -2, 2  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=10 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

But I don't like it as it makes the SpatialPoints agentsSP (fictional point) to a specific location in the real world. It seems that to use the over() function (from sp) later in the inRadius() function, a projected CRS is also necessary for all the Spatial objects in order to work. Is there a way to make the buffer() function work like before without specifying a CRS?

Then, there's a second problem, with over() this time.

pBuffer <- raster::buffer(agentsSP, dissolve = FALSE, width = 2)
sp1 <- SpatialPoints(coords = patches(w1))
proj4string(sp1) <-CRS("+proj=utm +zone=10 +datum=WGS84")
over(pBuffer, sp1, returnList = TRUE)
$`1`
[1] 16 17 21 22

16, 17, 21, 22 are the (rank) numbers of the points from sp1 which are INSIDE pBuffer whereas before, were selected the points inside and intersecting. Now only the points [0,1], [1,1], [0,0] and [1,0] are selected whereas before the points [0,2] and [2,0] would also be selected.

[achubaty]

Hmmm...I think this highlights a potential issue: here we are using rasters in a non-georeferenced way, which is not how they are intended to be used (i.e., in the GIS world). I don't know of, and can't easily find, a CRS that is not georeferenced. I think what we really want in this instance is a matrix/array -- these are non-spatial.

More generally, this gets at the fact that their are two types of spatial simulation: those that are georeferenced and those that aren't. How does NetLogo treat this distinction (if at all)? Once we are in the R-spatial world we need to deal with it. This may be a point worth highlighting in the MS.

[SarahBauduin]

There is no georeferencing at all in NetLogo, the world is treated like a matrix/array made up of squares and there is no way of linking some sort of coordinates or coordinate system. Which is what we reproduced in NetLogoR, but now it's causing issues apparently as we use Spatial objects that need coordinate systems. Spatial objects are only used inside functions and I think are never returned as result of a function. I know it's not pretty but can we use a random CRS in the function (like I did in the example)?

Note: I don't think it changes the problem but we only use Spatial objects (from sp) in the functions. We don't use raster (not that I recall) as we only use the worldMatrix/Array classes we created.

[eliotmcintire]

I have made a few changes that address both problems. Please give your thoughts:

Problem:

1. buffer needing a CRS

   Use CRS("+proj=utm") This give the units to be meters, which is critical, yet doesn't place it anywhere on the real planet Earth.

2. buffer now performing a slightly different operation, i.e,. it rounds incorrectly compared to previously, which is equivalent to the problem of < vs <= . inRadius previously was using <= width. With the latest changes to raster::buffer, this seems to now perform a < width, which is not what we want.

   Since this appears to be a floating point issue (though I didn't actually check that) or will become one if we are needing to calculate the difference between < and <=, we can resolve this by either using <= instead of the < or by adding a small number. Since we can't change the operator used in the internals of raster::buffer where the calculation occurs, I suggest we add a very small number. This is suboptimal, but it should be adequate in most cases. We can indicate this in the documentation for inRadius so users know.

It is passing all tests again on my computer. I will push to development shortly.

[eliotmcintire]

Passing on Travis-CI.

[achubaty]

Eliot's changes made in d46b99e3e8cdfd3edfe3ec0515babecec4709c7d