# Error while comparing two raster layer in R

I have got a geographic point distribution with N prehistoric settlements with a size of x and want to compare the distribution with an economic simulation, which produces M settlements (Usally N > M). I've done about 110 runs with different parameters and now I'm looking for a best-fit model. Therefore I calculate the Kernel Density Estimation for the real settlemts (x_dens), the simulation results (y_dens) and for the M-biggest real settlements (xM_dens) and want to compare the cell values of two kernel density estimations in R.

As a first step, i am using the function density.ppp from the Package "spatstat".

``````x_dens <- density(x_ppp, sd, eps=rw, edge=TRUE, at="pixels")
x_ras <- raster(x_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

y_dens <- density(y_ppp, sd, eps=rw, edge=TRUE, at="pixels")
y_ras <- raster(y_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")

xM_dens <- density(xM_ppp, sd, eps=rw, edge=TRUE, at="pixels")
xM_ras <- raster(xM_dens, crs="+proj=utm +zone=35 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0")
``````

The next step includes the calculation of the cell differences for x with xM, to see which would be minimum error for my M simulated settlements:

``````xxM_comp <- sqrt(((x_ras - xM_ras)^2 / x_ras))
xxM_comp [!is.finite(xxM_comp )] <- 0
xxM_comp_result <- sum(xxM_comp @data@values)
``````

As a next step, i'm doing the same for my Y-simulated settlements:

``````xy_comp <- sqrt(((x_ras - y_ras)^2 / x_ras))
xy_comp [!is.finite(xy_comp )] <- 0
xy_comp_result <- sum(xy_comp @data@values)
``````

To evaluate the best-fit model, i caluclate:

``````results <- xy_comp_results/xxM_comp_result
``````

xy_comp shall result in a new raster layer, which shows the differences in both point pattern, but the resulting cell values are enormously high at some points and as a consequence xy_comp_result has a illogically high value. It is due to the fact, that some pixel values of x_ras much to small compared with y_ras. Do you have got any suggestions, how to solve this problem. I've tried to remove outliers, but I did not get proper results.

The methodology is described in D. Stelder, Where do Cities Form? A Geographical Agglomeration Model for Europa, Journal of Regional Science 45 (4), 2005, 657-679 (http://www.regroningen.nl/stelder/doc/JRS_nov2005_c.pdf) Starting with page 669.

• What about extent and resolution of your rasters (x_ras,y_ras)? You need x_ras,y_ras to have the same extent and resolution, and then if you only want the differences you can try diff_rast = x_ras - y_ras. Commented Jan 26, 2016 at 14:41
• Thank you very much for your answer! x_ras and y_ras have got the same resolution as well as the same extent. I need the result of x_comp_result for further calculation. I want to compare a prehistoric settlement pattern with the results of an (economic) simulation. Commented Jan 26, 2016 at 15:00
• Have you thought of standard deviation? You can create a stack with your rasters: stack_rasters=stack(x_ras,y_ras). Then you can use the standard deviation like this: difference=calc(stack_rasters, fun=sd). (not tested). Commented Jan 26, 2016 at 15:19
• Might be a very good idea, but now i get following error: "Error in (function (classes, fdef, mtable): unable to find an inherited method for function ‘calc’ for signature ‘"RasterStack", "numeric"’" Commented Jan 26, 2016 at 15:45
• I believe that you are, in fact, using a ppp object and the density function from spatstat, not raster. This is not a standard KDE and should not be treated as such. The results are an isotropic density and represent the intensity function and an expectation of the spatial process. Across scales, these intensity functions are not comparable. Commented Jan 26, 2016 at 21:14