How to choose raster interpolation method?

Question

I work with a dataset of ~1.3 million residential buildings; each of them is described by a normally distributed variable ranging from 0 to 100. The dataset covers whole Switzerland so there are obviously areas with both very dense and very sparse point density.

The limitation imposed by our data provider is that we cannot release datasets more detailed than 5 buildings.

I was thinking about converting this point data into a raster surface. My main goal in this exercise would be to create the surface that minimizes the error that will arise when somebody tries to get the point value from this raster by extracting data in their's point locations of interest (as compared to extracting it from the original dataset of point locations).

I'd like to know what kind of method should be best suited to do that. There are quite few to choose from in the Spatial Analyst toolbox and even more in Geostatistical one so I'd very much appreciate some help in where to start.

Answer 1

It appears this question is related to an earlier one that asks about disguising such data using an irregular grid. If we accept that a regular grid will be used, then it seems that

• Most cells should be large enough to cover five or more buildings and

• When cells do not cover five buildings, their values should be changed in unpredictable (but controlled) ways.

How you measure the error determines the best solution. Let the value to be computed in a cell be y and let the values of the buildings located within (or at least overlapping) that cell be x1, x2, ..., xk. Further, suppose that each building has a non-negative "level of interest" (which might be proportional to the number of residents in it), to be used as a surrogate for the expected number of times your raster will be used to estimate that building's value. Let's call these levels w1, ..., wk respectively and let w designate their (nonzero) sum.

1. The average absolute error is the arithmetic mean of the error sizes |y - xi| as i ranges over the building indexes. This is minimized by choosing y to be the median of the xi.

2. The maximum error is the largest among max(|y - xi|) as i ranges over the building indexes. This is minimized by choosing y to be the midrange of (xi) (average of the max and min). This is heavily influenced by just a single outlying value, though, so the median might be preferable.

3. The expected error is the weighted average of |y - xi| with the weights given by wi / w. This is minimized by taking y to be a weighted median of the xi (but no GIS will do this calculation for you--you will have to use a statistical or mathematical package like R or Mathematica for such work.)

4. The expected squared error is the weighted average of (y - xi)^2. It is minimized by taking y to be the weighted mean of the xi, equal to the sum of wi xi / w.

You might be happy with (1) or (2) due to their simplicity and direct interpretation; I included (3) and (4) to give some sense of the options. To implement (1), you could begin by gridding all the data with a cell size so small that each building occupies its own cell. (At around 200 by 300 Km in extent, a cellsize of, say, 5 m would require an enormous grid of 40,000 by 60,000 cells, but only about a million of them would be occupied, thereby requiring only around 10 MB disk storage in native Arc format if you take care to store the values as integers.) Aggregate this grid to a larger cellsize using the Median option. (The cellsize of the aggregated grid would likely be around 100m, giving a country-wide grid of 2000 by 3000 cells: sufficiently small to make the procedures described below not only practicable, but rather quick to execute.)

You should also aggregate a binary indicator grid of the buildings--requesting the Sum this time--in order to count the numbers of buildings per cell. For aggregate cells with counts less than 5, randomly perturb the median. Do this with a Con operation. An effective, albeit slightly complicated, choice for the perturbation would be to add normally-distributed noise to the logit of the value (scaled from 0 to 1 rather than 0 to 100): this will guarantee a result that still lies between 0 and 100. You might also slightly perturb all the cells so that nobody can distinguish perturbed cells from unperturbed cells by inspecting the least significant digits.

The workflow for this "logistic perturbation" procedure, then, is as follows. It depends on two parameters: "sigma" is the amount of perturbation of the cells that need it and "epsilon" is the minimal amount to perturb all cells. Both are non-negative numbers. Experiment on small subgrids, starting with sigma=0.15 and epsilon = 0.01, and vary these parameters until results are satisfactory. (Setting epsilon to zero will eliminate the perturbation altogether for such cells.)

1. Begin with a grid [Z] of median values (all in the range from 0 to 100) and another grid [N] counting the numbers of buildings in each cell, both produced by Aggregate.

2. Create a grid for the amounts of perturbation using a Con command like

   Con["N" < 5, sigma, epsilon]
   

3. Generate the Normally-distributed perturbation by taking a grid of unit normal variates (use CreateNormalRaster and multiplying it by the previous grid. Call the result "e", say.

4. Compute the perturbed logits of the values as

   [Logit] = log("z" / (100 - "z")) + "e"
   

5. Convert back to values in the 0..100 range:

   100 / (1 + exp(-"logit"))
   

To illustrate, here is some R code to create small sample aggregate grids, perturb them, and compare the perturbed to the original values.

ncol <- 30; nrow <- 20
seed.random <- 17
x <- rpois(ncol * nrow, 5)
y <- floor(100 / (1 + exp(-(rnorm(ncol * nrow, mean = -2, sd = 1/sqrt(x))))))

sigma <- 0.15
epsilon <- 0.01

e <- rnorm(ncol*nrow, sd = ((x < 5)*sigma + (x >= 5)*epsilon))
logit <- log(y / (100 - y)) + e
y0 <- 100 / (1 + exp(-logit))

library(raster)
z <- matrix(y, ncol=ncol)
n <- matrix(x, ncol=ncol)
z0 <- matrix(y0, ncol=ncol)

par(mfrow=c(2,2))
n.r <- raster(n)
plot(n.r, main="Counts of residences [N]")

z.r <- raster(z)
plot(z.r, main="Median values [Z]")

z0.r <- raster(z0)
plot(z0.r, main="Perturbed median values")

plot(y, y0, type="n", xlab="Original medians", ylab="Perturbed medians",
     main="Perturbed vs. original medians")
points(y[x < 5], y0[x < 5], col="Red")
points(y[x >= 5], y0[x >= 5], pch=19)