0

I would like to generate different rasters from kernel density functions with the same resolution, so that I can compare them.

Take this example from here:

library(MASS)
library(raster)
set.seed(123); xy = cbind(rnorm(1000), rnorm(1000))
k = kde2d(xy[,1],xy[,2],h=.4,n=100)
r1 = raster(k)

xy = cbind(rnorm(1000), rnorm(1000))
k = kde2d(xy[,1],xy[,2],h=.4,n=100)
r2 = raster(k)

par(mfrow = c(1,2))
plot(r1); plot(r2)

r1r2

r1 - r2
#Error in compareRaster(e1, e2, extent = FALSE, rowcol = FALSE, crs = TRUE,  : 
#  different origin

I know I can resample the rasters, but I would like to define a standard raster/stars object

I can compare as stars objects.

Is it correct?

s1 <- st_as_stars(r1)
s2 <- st_as_stars(r2)
s3 <- s1 - s2

plot(s3, axes = T, breaks = "equal")

stars diff

0

1 Answer 1

3

To evaluate a kde2d over the same basis grid as a raster r1, use this:

k2 = kde2d(xy[,1],xy[,2],h=.4,
   n=c(ncol(r1),nrow(r1)),
   lims=c(range(coordinates(r1)[,1]),
          range(coordinates(r1)[,2])))

Note the use of coordinates here which gets the cell centres, the extent function could be used but you have to subtract the cell half-width and half-height to get the cell centres. There's probably better ways of doing this (for one thing I'm calling coordinates twice) but this works as a proof of concept.

Testing this on a 200x100 raster with this extent etc:

> r1
class      : RasterLayer 
dimensions : 200, 100, 20000  (nrow, ncol, ncell)
resolution : 0.02020202, 0.0201005  (x, y)
extent     : -1.010101, 1.010101, -2.01005, 2.01005  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 1.011512e-05, 0.2004785  (min, max)

> k2 = kde2d(xy[,1],xy[,2],h=.4,n=c(ncol(r1),nrow(r1)), lims=c(range(coordinates(r1)[,1]), range(coordinates(r1)[,2])))

Produces a raster with the same grid geometry:

> r2 = raster(k2)
> r2
class      : RasterLayer 
dimensions : 200, 100, 20000  (nrow, ncol, ncell)
resolution : 0.02020202, 0.0201005  (x, y)
extent     : -1.010101, 1.010101, -2.01005, 2.01005  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 6.548952e-05, 0.2128547  (min, max)

And these are compatible:

> plot(r2-r1)
>

no errors.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.