# Generating Kernel density raster/stars with same dimensions

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)
``````

``````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")
``````

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.