+
+
+
+
+
+
+
+
+
2023-08-03
+
+
1. Preparing example data
+
library(SpatialKWD)
+library(sp)
+library(raster)
+library(sdcSpatial)
+
The issues are illustrated by means of an example data set from the sdcSpatial package called ‘dwellings’. The first step is to create the unprotected raster map, which forms the ground truth.
+
The example uses a map of counts (no. of units per grid cell), with 200m \(\times\) 200m grid cells in a 62 \(\times\) 61 cells raster map.
+
data("dwellings")
+
+dw <- sdcSpatial::sdc_raster(dwellings[, c("x", "y")],
+ variable = 1, r = 200,
+ max_risk = .8, min_count = 10)
+
+ras_orig <- dw$value$count # original raster
+
+raster::plot(ras_orig, main = "no. of dwellings") # map
+
+
ras_orig # raster info
+
## class : RasterLayer
+## dimensions : 62, 61, 3782 (nrow, ncol, ncell)
+## resolution : 200, 200 (x, y)
+## extent : 149400, 161600, 457800, 470200 (xmin, xmax, ymin, ymax)
+## crs : NA
+## source : memory
+## names : count
+## values : 1, 329 (min, max)
+
Next, we create the altered raster maps using three protection methods: 1) suppression of sensitive cells (remove_sensitive), 2) quadtree method (protect_quadtree), 3) kernel smoothing (protect_smooth).
+
# apply protection methods
+dw_rm <- remove_sensitive(dw)
+dw_qt <- protect_quadtree(dw, max_zoom = 2)
+dw_sm <- protect_smooth(dw, bw = 200)
+
+# extract protected rasters from sdcSpatial object
+ras_rm <- dw_rm$value$count
+ras_qt <- dw_qt$value$count
+ras_sm <- dw_sm$value$count
+
Following up, we select a square focus area. It is identified by its centroid and radius.
+
# Focus area (square)
+
+fa_ctrd <- c(155100, 464100) # center of focus area
+fa_rads <- 3000 # radius of focus area
+
+
+# maps
+
+par(mfrow = c(2, 2))
+
+plot(ras_orig, main = "original")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_rm, main = "cell suppression")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_qt, main = "quad tree method")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_sm, main = "smoothing method")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+
+
par(mfrow = c(1, 1))
+
+
+
2. FocusArea sometimes returns zero for clear differences
+
We continue by evaluating the difference between altered maps and ground truth with respect to the chosen focus area.
+
This involves the following steps: 1) extract the center coordinates of tiles (raster::xyFromCell), 2) extract the tile-values (weights) (raster::getValues) [where NA-values are treated as zeros].
+
# extract coordinates from raster
+xy <- raster::xyFromCell(ras_orig, 1:ncell(ras_orig))
+
+# extract cell values from rasters
+w_orig <- raster::getValues(ras_orig)
+w_rm <- raster::getValues(ras_rm)
+w_qt <- raster::getValues(ras_qt)
+w_sm <- raster::getValues(ras_sm)
+
+# treat NA-nodes as having zero weight
+w_orig[is.na(w_orig)] <- 0
+w_rm[is.na(w_rm)] <- 0
+w_qt[is.na(w_qt)] <- 0
+w_sm[is.na(w_sm)] <- 0
+
Using these inputs we calculate the first set of KWDs.
+
kwd_rm <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_rm),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
## FocusArea, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_qt <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_qt),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
## FocusArea, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_sm <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_sm),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
## FocusArea, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
Even though the suppression map is clearly different (all cells with < 10 units are removed, see maps above), we get a KWD of zero.
+
c(kwd_rm$distance, kwd_qt$distance, kwd_sm$distance)
+
## [1] 0.000 6514.312 2332.205
+
Second try: The description of the focusArea function tells us that “All the weights outside the focus area should be equal to zero.” We create new rasters, where all tiles outside the focus area are assigned zero values.
+
ex_fa <- extent(fa_ctrd[1] - (fa_rads - 100), fa_ctrd[1] + (fa_rads - 100),
+ fa_ctrd[2] - (fa_rads - 100), fa_ctrd[2] + (fa_rads - 100))
+clls_fa <- cellsFromExtent(ras_orig, ex_fa)
+
+ras_rm_fa <- ras_rm
+ras_qt_fa <- ras_qt
+ras_sm_fa <- ras_sm
+ras_rm_fa[-clls_fa] <- 0
+ras_qt_fa[-clls_fa] <- 0
+ras_sm_fa[-clls_fa] <- 0
+
+
+par(mfrow = c(2, 2))
+
+plot(ras_orig, main = "original")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_rm_fa, main = "cell suppression (focus area)")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_qt_fa, main = "quad tree method (focus area)")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_sm_fa, main = "smoothing method (focus area)")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+
+
par(mfrow = c(1, 1))
+
With these we re-calculate the KWDs.
+
# updated weights (zero outside focus area)
+w_rm_fa <- getValues(ras_rm_fa)
+w_qt_fa <- getValues(ras_qt_fa)
+w_sm_fa <- getValues(ras_sm_fa)
+w_rm_fa[is.na(w_rm_fa)] <- 0
+w_qt_fa[is.na(w_qt_fa)] <- 0
+w_sm_fa[is.na(w_sm_fa)] <- 0
+
+kwd_rm_fa <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_rm_fa),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
## FocusArea, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_qt_fa <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_qt_fa),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
## FocusArea, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_sm_fa <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_sm_fa),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
## FocusArea, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
We still get a KWD of zero for the ‘removal’ method. Furthermore, we now also get zero for smoothing.
+
c(kwd_rm_fa$distance, kwd_qt_fa$distance, kwd_sm_fa$distance)
+
## [1] 0.000 2748.876 0.000
+
+
+
3. How to apply a fixed cost to missing units of mass?
+
Suppose now we look at the whole map rather than focus areas. Then, we might get a mass mismatch, especially with removal:
+
c(sum(w_orig), sum(w_rm))
+
## [1] 90603 88530
+
There are ca. 2000 units missing after removal. We want to try solving the problem by assigning each a cost of half the map diagonal.
+
# compute length of the diagonal in cells
+diag_cells <- sqrt(ncol(ras_orig)^2 + nrow(ras_orig)^2)
+# ... in meters
+diag_metrs <- diag_cells * 200
+
+# solve with extra artificial bin (?)
+kwd_rm_bal <- compareOneToOne(Coordinates = xy, Weights = cbind(w_orig, w_rm),
+ L = 3, recode = TRUE,
+ unbalanced = TRUE, unbal_cost = diag_metrs/2)
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
Again, this yields a distance of zero.
+
kwd_rm_bal$distance
+
## [1] 0
+
Do we get a distance of zero, because all non-zero tiles after removal match exactly? (I.e. for any given tile the mass is missing completely or not at all.)
+
Given that this is the case, we could, in principle, assign the costs manually (2000 \(\times c\) ). Still, intuitively the function should not return zero (?)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Preparing example data
+
+
+
+
library(SpatialKWD)
+library(sp)
+library(raster)
+library(sdcSpatial)
+
+
+
+
The issues are illustrated by means of an example data set from the sdcSpatial package called ‘dwellings’. The first step is to create the unprotected raster map, which forms the ground truth.
+
The example uses a map of counts (no. of units per grid cell), with 200m \(\times\) 200m grid cells in a 62 \(\times\) 61 cells raster map.
+
+
+
+
data("dwellings")
+
+dw <- sdcSpatial::sdc_raster(dwellings[, c("x", "y")],
+ variable = 1, r = 200,
+ max_risk = .8, min_count = 5)
+
+ras_orig <- dw$value$count # original raster
+
+raster::plot(ras_orig) # map
+ras_orig # raster info
+
+
+
+
Next, we create the altered raster maps using three protection methods: 1) suppression of sensitive cells (remove_sensitive), 2) quadtree method (protect_quadtree), 3) kernel smoothing (protect_smooth).
+
+
+
+
# apply protection methods
+dw_rm <- remove_sensitive(dw)
+dw_qt <- protect_quadtree(dw, max_zoom = 2)
+dw_sm <- protect_smooth(dw, bw = 200)
+
+# extract protected rasters from sdcSpatial object
+ras_rm <- dw_rm$value$count
+ras_qt <- dw_qt$value$count
+ras_sm <- dw_sm$value$count
+
+
+
+
Following up, we select a square focus area. It is identified by its centroid and radius.
+
+
+
+
# Focus area (square)
+
+fa_ctrd <- c(155100, 464100) # center of focus area
+fa_rads <- 3000 # radius of focus area
+
+
+par(mfrow = c(2, 2))
+
+plot(ras_orig, main = "original")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_rm, main = "cell suppression")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_qt, main = "quad tree method")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_sm, main = "smoothing method")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+
+par(mfrow = c(1, 1))
+
+
+
+
+
+
+
# ----- KWDs -----
+
+# extract coordinates from raster
+xy <- raster::xyFromCell(ras_orig, 1:ncell(ras_orig))
+
+# extract cell values from rasters
+w_orig <- raster::getValues(ras_orig)
+w_rm <- raster::getValues(ras_rm)
+w_qt <- raster::getValues(ras_qt)
+w_sm <- raster::getValues(ras_sm)
+
+# treat NA-nodes as having zero weight
+w_orig[is.na(w_orig)] <- 0
+w_rm[is.na(w_rm)] <- 0
+w_qt[is.na(w_qt)] <- 0
+w_sm[is.na(w_sm)] <- 0
+
+kwd_rm <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_rm),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+kwd_qt <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_qt),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+kwd_sm <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_sm),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
+# removal gives dist. of 0, although clearly different
+c(kwd_rm$distance, kwd_qt$distance, kwd_sm$distance)
+
+
+# try with weights outside FA set to zero
+ex_fa <- extent(fa_ctrd[1] - (fa_rads - 100), fa_ctrd[1] + (fa_rads - 100),
+ fa_ctrd[2] - (fa_rads - 100), fa_ctrd[2] + (fa_rads - 100))
+clls_fa <- cellsFromExtent(ras_orig, ex_fa)
+
+ras_rm_fa <- ras_rm
+ras_qt_fa <- ras_qt
+ras_sm_fa <- ras_sm
+ras_rm_fa[-clls_fa] <- 0
+ras_qt_fa[-clls_fa] <- 0
+ras_sm_fa[-clls_fa] <- 0
+
+
+par(mfrow = c(2, 2))
+
+plot(ras_orig, main = "original")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_rm_fa, main = "cell suppression (focus area)")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_qt_fa, main = "quad tree method (focus area)")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+plot(ras_sm_fa, main = "smoothing method (focus area)")
+rect(xleft = fa_ctrd[1] - fa_rads, xright = fa_ctrd[1] + fa_rads,
+ ybottom = fa_ctrd[2] - fa_rads, ytop = fa_ctrd[2] + fa_rads, border = "red")
+
+par(mfrow = c(1, 1))
+
+# updated weights (zero outside focus area)
+w_rm_fa <- getValues(ras_rm_fa)
+w_qt_fa <- getValues(ras_qt_fa)
+w_sm_fa <- getValues(ras_sm_fa)
+w_rm_fa[is.na(w_rm_fa)] <- 0
+w_qt_fa[is.na(w_qt_fa)] <- 0
+w_sm_fa[is.na(w_sm_fa)] <- 0
+
+kwd_rm_fa <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_rm_fa),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+kwd_qt_fa <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_qt_fa),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+kwd_sm_fa <- SpatialKWD::focusArea(Coordinates = xy, Weights = cbind(w_orig, w_sm_fa),
+ L = 3, recode = TRUE, area = "linf",
+ x = fa_ctrd[1], y = fa_ctrd[2], radius = fa_rads)
+
+c(kwd_rm_fa$distance, kwd_qt_fa$distance, kwd_sm_fa$distance)
+
+
+
+
+
+
+
LS0tCnRpdGxlOiAiUmVwcm9kdWNpYmxlIGV4YW1wbGUgb2YgY3VycmVudCBvcGVuIGlzc3VlcyB3aXRoIHRoZSBTcGF0aWFsS1dEIHBhY2thZ2UiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMjIFByZXBhcmluZyBleGFtcGxlIGRhdGEKCmBgYHtyfQpsaWJyYXJ5KFNwYXRpYWxLV0QpCmxpYnJhcnkoc3ApCmxpYnJhcnkocmFzdGVyKQpsaWJyYXJ5KHNkY1NwYXRpYWwpCmBgYAoKVGhlIGlzc3VlcyBhcmUgaWxsdXN0cmF0ZWQgYnkgbWVhbnMgb2YgYW4gZXhhbXBsZSBkYXRhIHNldCBmcm9tIHRoZSBzZGNTcGF0aWFsIHBhY2thZ2UgY2FsbGVkCidkd2VsbGluZ3MnLiBUaGUgZmlyc3Qgc3RlcCBpcyB0byBjcmVhdGUgdGhlIHVucHJvdGVjdGVkIHJhc3RlciBtYXAsIHdoaWNoIGZvcm1zIHRoZQpncm91bmQgdHJ1dGguIAoKVGhlIGV4YW1wbGUgdXNlcyBhIG1hcCBvZiBjb3VudHMgKG5vLiBvZiB1bml0cyBwZXIgZ3JpZCBjZWxsKSwgd2l0aCAyMDBtICRcdGltZXMkIAoyMDBtIGdyaWQgY2VsbHMgaW4gYSA2MiAkXHRpbWVzJCA2MSBjZWxscyByYXN0ZXIgbWFwLgoKYGBge3J9CmRhdGEoImR3ZWxsaW5ncyIpCgpkdyA8LSBzZGNTcGF0aWFsOjpzZGNfcmFzdGVyKGR3ZWxsaW5nc1ssIGMoIngiLCAieSIpXSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdmFyaWFibGUgPSAxLCByID0gMjAwLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIG1heF9yaXNrID0gLjgsIG1pbl9jb3VudCA9IDUpCgpyYXNfb3JpZyA8LSBkdyR2YWx1ZSRjb3VudCAjIG9yaWdpbmFsIHJhc3RlcgoKcmFzdGVyOjpwbG90KHJhc19vcmlnKSAjIG1hcApyYXNfb3JpZyAgICAgICAgICAgICAgICMgcmFzdGVyIGluZm8KYGBgCgpOZXh0LCB3ZSBjcmVhdGUgdGhlIGFsdGVyZWQgcmFzdGVyIG1hcHMgdXNpbmcgdGhyZWUgcHJvdGVjdGlvbiBtZXRob2RzOgoxKSBzdXBwcmVzc2lvbiBvZiBzZW5zaXRpdmUgY2VsbHMgKHJlbW92ZV9zZW5zaXRpdmUpLAoyKSBxdWFkdHJlZSBtZXRob2QgKHByb3RlY3RfcXVhZHRyZWUpLAozKSBrZXJuZWwgc21vb3RoaW5nIChwcm90ZWN0X3Ntb290aCkuCgpgYGB7cn0KIyBhcHBseSBwcm90ZWN0aW9uIG1ldGhvZHMKZHdfcm0gPC0gcmVtb3ZlX3NlbnNpdGl2ZShkdykKZHdfcXQgPC0gcHJvdGVjdF9xdWFkdHJlZShkdywgbWF4X3pvb20gPSAyKQpkd19zbSA8LSBwcm90ZWN0X3Ntb290aChkdywgYncgPSAyMDApCgojIGV4dHJhY3QgcHJvdGVjdGVkIHJhc3RlcnMgZnJvbSBzZGNTcGF0aWFsIG9iamVjdApyYXNfcm0gPC0gZHdfcm0kdmFsdWUkY291bnQKcmFzX3F0IDwtIGR3X3F0JHZhbHVlJGNvdW50CnJhc19zbSA8LSBkd19zbSR2YWx1ZSRjb3VudApgYGAKCkZvbGxvd2luZyB1cCwgd2Ugc2VsZWN0IGEgc3F1YXJlIGZvY3VzIGFyZWEuIEl0IGlzIGlkZW50aWZpZWQgYnkgaXRzIGNlbnRyb2lkIGFuZCAKcmFkaXVzLgoKYGBge3J9CiMgRm9jdXMgYXJlYSAoc3F1YXJlKQoKZmFfY3RyZCA8LSBjKDE1NTEwMCwgNDY0MTAwKSAjIGNlbnRlciBvZiBmb2N1cyBhcmVhCmZhX3JhZHMgPC0gMzAwMCAgICAgICAgICAgICAgIyByYWRpdXMgb2YgZm9jdXMgYXJlYQoKCnBhcihtZnJvdyA9IGMoMiwgMikpCgpwbG90KHJhc19vcmlnLCBtYWluID0gIm9yaWdpbmFsIikKcmVjdCh4bGVmdCA9IGZhX2N0cmRbMV0gLSBmYV9yYWRzLCB4cmlnaHQgPSBmYV9jdHJkWzFdICsgZmFfcmFkcywKICAgICB5Ym90dG9tID0gZmFfY3RyZFsyXSAtIGZhX3JhZHMsIHl0b3AgPSBmYV9jdHJkWzJdICsgZmFfcmFkcywgYm9yZGVyID0gInJlZCIpCnBsb3QocmFzX3JtLCBtYWluID0gImNlbGwgc3VwcHJlc3Npb24iKQpyZWN0KHhsZWZ0ID0gZmFfY3RyZFsxXSAtIGZhX3JhZHMsIHhyaWdodCA9IGZhX2N0cmRbMV0gKyBmYV9yYWRzLAogICAgIHlib3R0b20gPSBmYV9jdHJkWzJdIC0gZmFfcmFkcywgeXRvcCA9IGZhX2N0cmRbMl0gKyBmYV9yYWRzLCBib3JkZXIgPSAicmVkIikKcGxvdChyYXNfcXQsIG1haW4gPSAicXVhZCB0cmVlIG1ldGhvZCIpCnJlY3QoeGxlZnQgPSBmYV9jdHJkWzFdIC0gZmFfcmFkcywgeHJpZ2h0ID0gZmFfY3RyZFsxXSArIGZhX3JhZHMsCiAgICAgeWJvdHRvbSA9IGZhX2N0cmRbMl0gLSBmYV9yYWRzLCB5dG9wID0gZmFfY3RyZFsyXSArIGZhX3JhZHMsIGJvcmRlciA9ICJyZWQiKQpwbG90KHJhc19zbSwgbWFpbiA9ICJzbW9vdGhpbmcgbWV0aG9kIikKcmVjdCh4bGVmdCA9IGZhX2N0cmRbMV0gLSBmYV9yYWRzLCB4cmlnaHQgPSBmYV9jdHJkWzFdICsgZmFfcmFkcywKICAgICB5Ym90dG9tID0gZmFfY3RyZFsyXSAtIGZhX3JhZHMsIHl0b3AgPSBmYV9jdHJkWzJdICsgZmFfcmFkcywgYm9yZGVyID0gInJlZCIpCgpwYXIobWZyb3cgPSBjKDEsIDEpKQpgYGAKCgoKYGBge3J9CiMgLS0tLS0gS1dEcyAtLS0tLQoKIyBleHRyYWN0IGNvb3JkaW5hdGVzIGZyb20gcmFzdGVyCnh5IDwtIHJhc3Rlcjo6eHlGcm9tQ2VsbChyYXNfb3JpZywgMTpuY2VsbChyYXNfb3JpZykpCgojIGV4dHJhY3QgY2VsbCB2YWx1ZXMgZnJvbSByYXN0ZXJzCndfb3JpZyA8LSByYXN0ZXI6OmdldFZhbHVlcyhyYXNfb3JpZykKd19ybSAgIDwtIHJhc3Rlcjo6Z2V0VmFsdWVzKHJhc19ybSkKd19xdCAgIDwtIHJhc3Rlcjo6Z2V0VmFsdWVzKHJhc19xdCkKd19zbSAgIDwtIHJhc3Rlcjo6Z2V0VmFsdWVzKHJhc19zbSkKCiMgdHJlYXQgTkEtbm9kZXMgYXMgaGF2aW5nIHplcm8gd2VpZ2h0Cndfb3JpZ1tpcy5uYSh3X29yaWcpXSA8LSAwCndfcm1baXMubmEod19ybSldIDwtIDAKd19xdFtpcy5uYSh3X3F0KV0gPC0gMAp3X3NtW2lzLm5hKHdfc20pXSA8LSAwCgprd2Rfcm0gPC0gU3BhdGlhbEtXRDo6Zm9jdXNBcmVhKENvb3JkaW5hdGVzID0geHksIFdlaWdodHMgPSBjYmluZCh3X29yaWcsIHdfcm0pLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIEwgPSAzLCByZWNvZGUgPSBUUlVFLCBhcmVhID0gImxpbmYiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHggPSBmYV9jdHJkWzFdLCB5ID0gZmFfY3RyZFsyXSwgcmFkaXVzID0gZmFfcmFkcykKa3dkX3F0IDwtIFNwYXRpYWxLV0Q6OmZvY3VzQXJlYShDb29yZGluYXRlcyA9IHh5LCBXZWlnaHRzID0gY2JpbmQod19vcmlnLCB3X3F0KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBMID0gMywgcmVjb2RlID0gVFJVRSwgYXJlYSA9ICJsaW5mIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB4ID0gZmFfY3RyZFsxXSwgeSA9IGZhX2N0cmRbMl0sIHJhZGl1cyA9IGZhX3JhZHMpCmt3ZF9zbSA8LSBTcGF0aWFsS1dEOjpmb2N1c0FyZWEoQ29vcmRpbmF0ZXMgPSB4eSwgV2VpZ2h0cyA9IGNiaW5kKHdfb3JpZywgd19zbSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgTCA9IDMsIHJlY29kZSA9IFRSVUUsIGFyZWEgPSAibGluZiIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeCA9IGZhX2N0cmRbMV0sIHkgPSBmYV9jdHJkWzJdLCByYWRpdXMgPSBmYV9yYWRzKQoKIyByZW1vdmFsIGdpdmVzIGRpc3QuIG9mIDAsIGFsdGhvdWdoIGNsZWFybHkgZGlmZmVyZW50CmMoa3dkX3JtJGRpc3RhbmNlLCBrd2RfcXQkZGlzdGFuY2UsIGt3ZF9zbSRkaXN0YW5jZSkKCgojIHRyeSB3aXRoIHdlaWdodHMgb3V0c2lkZSBGQSBzZXQgdG8gemVybwpleF9mYSA8LSBleHRlbnQoZmFfY3RyZFsxXSAtIChmYV9yYWRzIC0gMTAwKSwgZmFfY3RyZFsxXSArIChmYV9yYWRzIC0gMTAwKSwKICAgICAgICAgICAgICAgIGZhX2N0cmRbMl0gLSAoZmFfcmFkcyAtIDEwMCksIGZhX2N0cmRbMl0gKyAoZmFfcmFkcyAtIDEwMCkpCmNsbHNfZmEgPC0gY2VsbHNGcm9tRXh0ZW50KHJhc19vcmlnLCBleF9mYSkKCnJhc19ybV9mYSA8LSByYXNfcm0KcmFzX3F0X2ZhIDwtIHJhc19xdApyYXNfc21fZmEgPC0gcmFzX3NtCnJhc19ybV9mYVstY2xsc19mYV0gPC0gMApyYXNfcXRfZmFbLWNsbHNfZmFdIDwtIDAKcmFzX3NtX2ZhWy1jbGxzX2ZhXSA8LSAwCgoKcGFyKG1mcm93ID0gYygyLCAyKSkKCnBsb3QocmFzX29yaWcsIG1haW4gPSAib3JpZ2luYWwiKQpyZWN0KHhsZWZ0ID0gZmFfY3RyZFsxXSAtIGZhX3JhZHMsIHhyaWdodCA9IGZhX2N0cmRbMV0gKyBmYV9yYWRzLAogICAgIHlib3R0b20gPSBmYV9jdHJkWzJdIC0gZmFfcmFkcywgeXRvcCA9IGZhX2N0cmRbMl0gKyBmYV9yYWRzLCBib3JkZXIgPSAicmVkIikKcGxvdChyYXNfcm1fZmEsIG1haW4gPSAiY2VsbCBzdXBwcmVzc2lvbiAoZm9jdXMgYXJlYSkiKQpyZWN0KHhsZWZ0ID0gZmFfY3RyZFsxXSAtIGZhX3JhZHMsIHhyaWdodCA9IGZhX2N0cmRbMV0gKyBmYV9yYWRzLAogICAgIHlib3R0b20gPSBmYV9jdHJkWzJdIC0gZmFfcmFkcywgeXRvcCA9IGZhX2N0cmRbMl0gKyBmYV9yYWRzLCBib3JkZXIgPSAicmVkIikKcGxvdChyYXNfcXRfZmEsIG1haW4gPSAicXVhZCB0cmVlIG1ldGhvZCAoZm9jdXMgYXJlYSkiKQpyZWN0KHhsZWZ0ID0gZmFfY3RyZFsxXSAtIGZhX3JhZHMsIHhyaWdodCA9IGZhX2N0cmRbMV0gKyBmYV9yYWRzLAogICAgIHlib3R0b20gPSBmYV9jdHJkWzJdIC0gZmFfcmFkcywgeXRvcCA9IGZhX2N0cmRbMl0gKyBmYV9yYWRzLCBib3JkZXIgPSAicmVkIikKcGxvdChyYXNfc21fZmEsIG1haW4gPSAic21vb3RoaW5nIG1ldGhvZCAoZm9jdXMgYXJlYSkiKQpyZWN0KHhsZWZ0ID0gZmFfY3RyZFsxXSAtIGZhX3JhZHMsIHhyaWdodCA9IGZhX2N0cmRbMV0gKyBmYV9yYWRzLAogICAgIHlib3R0b20gPSBmYV9jdHJkWzJdIC0gZmFfcmFkcywgeXRvcCA9IGZhX2N0cmRbMl0gKyBmYV9yYWRzLCBib3JkZXIgPSAicmVkIikKCnBhcihtZnJvdyA9IGMoMSwgMSkpCgojIHVwZGF0ZWQgd2VpZ2h0cyAoemVybyBvdXRzaWRlIGZvY3VzIGFyZWEpCndfcm1fZmEgPC0gZ2V0VmFsdWVzKHJhc19ybV9mYSkKd19xdF9mYSA8LSBnZXRWYWx1ZXMocmFzX3F0X2ZhKQp3X3NtX2ZhIDwtIGdldFZhbHVlcyhyYXNfc21fZmEpCndfcm1fZmFbaXMubmEod19ybV9mYSldIDwtIDAKd19xdF9mYVtpcy5uYSh3X3F0X2ZhKV0gPC0gMAp3X3NtX2ZhW2lzLm5hKHdfc21fZmEpXSA8LSAwCgprd2Rfcm1fZmEgPC0gU3BhdGlhbEtXRDo6Zm9jdXNBcmVhKENvb3JkaW5hdGVzID0geHksIFdlaWdodHMgPSBjYmluZCh3X29yaWcsIHdfcm1fZmEpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIEwgPSAzLCByZWNvZGUgPSBUUlVFLCBhcmVhID0gImxpbmYiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHggPSBmYV9jdHJkWzFdLCB5ID0gZmFfY3RyZFsyXSwgcmFkaXVzID0gZmFfcmFkcykKa3dkX3F0X2ZhIDwtIFNwYXRpYWxLV0Q6OmZvY3VzQXJlYShDb29yZGluYXRlcyA9IHh5LCBXZWlnaHRzID0gY2JpbmQod19vcmlnLCB3X3F0X2ZhKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBMID0gMywgcmVjb2RlID0gVFJVRSwgYXJlYSA9ICJsaW5mIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB4ID0gZmFfY3RyZFsxXSwgeSA9IGZhX2N0cmRbMl0sIHJhZGl1cyA9IGZhX3JhZHMpCmt3ZF9zbV9mYSA8LSBTcGF0aWFsS1dEOjpmb2N1c0FyZWEoQ29vcmRpbmF0ZXMgPSB4eSwgV2VpZ2h0cyA9IGNiaW5kKHdfb3JpZywgd19zbV9mYSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgTCA9IDMsIHJlY29kZSA9IFRSVUUsIGFyZWEgPSAibGluZiIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeCA9IGZhX2N0cmRbMV0sIHkgPSBmYV9jdHJkWzJdLCByYWRpdXMgPSBmYV9yYWRzKQoKYyhrd2Rfcm1fZmEkZGlzdGFuY2UsIGt3ZF9xdF9mYSRkaXN0YW5jZSwga3dkX3NtX2ZhJGRpc3RhbmNlKQpgYGAKCgoKCg==
+
+
+
+
+
+
+
+
+
+
+
+
+
We use pre-computed raster maps from 02_sdcSpatial_DE.R, see GitHub repo.
+
library(sdcSpatial)
+library(SpatialKWD)
+library(sp)
+library(raster)
+library(viridisLite)
+library(viridis)
+
+load("sdcRasters_SHARE.RData")
+
+pal <- rev(viridis(10))
+
+
Making focus areas
+
Selecting four focus areas by: centroid (center point in ETRS89-LAEA) and size (measured in cell width).
+
fa1_ctrd <- c(4122250, 3162750)
+fa2_ctrd <- c(4195750, 2993750)
+fa3_ctrd <- c(4535250, 3474750)
+fa4_ctrd <- c(4328750, 2751250)
+
+cellsize <- 500
+
+# focus area sizes
+fa_wcells <- c(55, 41, 75, 55)
+
+# find cell index of centroid in larger raster
+fa1_cntr <- c(raster::colFromX(pop_ras$value, x = fa1_ctrd[1]),
+ raster::rowFromY(pop_ras$value, y = fa1_ctrd[2]))
+fa2_cntr <- c(raster::colFromX(pop_ras$value, x = fa2_ctrd[1]),
+ raster::rowFromY(pop_ras$value, y = fa2_ctrd[2]))
+fa3_cntr <- c(raster::colFromX(pop_ras$value, x = fa3_ctrd[1]),
+ raster::rowFromY(pop_ras$value, y = fa3_ctrd[2]))
+fa4_cntr <- c(raster::colFromX(pop_ras$value, x = fa4_ctrd[1]),
+ raster::rowFromY(pop_ras$value, y = fa4_ctrd[2]))
+
+# create extent objects for chosen positions & sizes
+fa1_extt <- extent(pop_ras$value, fa1_cntr[2] - (fa_wcells[1] - 1)/2,
+ fa1_cntr[2] + (fa_wcells[1] - 1)/2,
+ fa1_cntr[1] - (fa_wcells[1] - 1)/2,
+ fa1_cntr[1] + (fa_wcells[1] - 1)/2)
+
+fa2_extt <- extent(pop_ras$value, fa2_cntr[2] - (fa_wcells[2] - 1)/2,
+ fa2_cntr[2] + (fa_wcells[2] - 1)/2,
+ fa2_cntr[1] - (fa_wcells[2] - 1)/2,
+ fa2_cntr[1] + (fa_wcells[2] - 1)/2)
+
+fa3_extt <- extent(pop_ras$value, fa3_cntr[2] - (fa_wcells[3] - 1)/2,
+ fa3_cntr[2] + (fa_wcells[3] - 1)/2,
+ fa3_cntr[1] - (fa_wcells[3] - 1)/2,
+ fa3_cntr[1] + (fa_wcells[3] - 1)/2)
+
+fa4_extt <- extent(pop_ras$value, fa4_cntr[2] - (fa_wcells[4] - 1)/2,
+ fa4_cntr[2] + (fa_wcells[4] - 1)/2,
+ fa4_cntr[1] - (fa_wcells[4] - 1)/2,
+ fa4_cntr[1] + (fa_wcells[4] - 1)/2)
+
+# crop overall raster to focus area for plotting & pop. counts
+fa1 <- fa2 <- fa3 <- fa4 <- pop_ras
+fa1$value <- crop(fa1$value, fa1_extt)
+fa2$value <- crop(fa2$value, fa2_extt)
+fa3$value <- crop(fa3$value, fa3_extt)
+fa4$value <- crop(fa4$value, fa4_extt)
+
+par(mfrow = c(2, 2))
+ plot(fa1, "count", sensitive = FALSE, col = pal, main = "Ruhr valley")
+ plot(fa2, "count", sensitive = FALSE, col = pal, main = "MZ + WI")
+ plot(fa3, "count", sensitive = FALSE, col = pal, main = "Strelasund")
+ plot(fa4, "count", sensitive = FALSE, col = pal, main = "Allgäu")
+
+
par(mfrow = c(1, 1))
+
+# calculate populations for subregions
+p_fa <- c(sum(getValues(fa1$value$count), na.rm = TRUE),
+ sum(getValues(fa2$value$count), na.rm = TRUE),
+ sum(getValues(fa3$value$count), na.rm = TRUE),
+ sum(getValues(fa4$value$count), na.rm = TRUE))
+
+fa_sizes <- cbind(fa_wcells,
+ fa_wcells^2,
+ fa_wcells^2 * cellsize^2 / 1e+6,
+ p_fa)
+colnames(fa_sizes) <- c("width (ncells)", "size (ncells)", "size (km^2)", "pop.")
+rownames(fa_sizes) <- c("fa1", "fa2", "fa3", "fa4")
+fa_sizes
+
## width (ncells) size (ncells) size (km^2) pop.
+## fa1 55 3025 756.25 1745062
+## fa2 41 1681 420.25 558055
+## fa3 75 5625 1406.25 107863
+## fa4 55 3025 756.25 106283
+
+
+
Risk and bin-by-bin measure for focus areas
+
We start by protecting focus areas, using the sdcSpatial algorithms. For this, wrapper functions are used:
+
wrap_qt <- function(x, max_zoom = Inf, enlarge = TRUE) {
+
+ if(enlarge) {
+ pad <- min(2^max_zoom - 1, 10) # maximum padding is 10 (arbitrary)
+ x$value <- extend(x$value, c(pad, pad))
+ x <- sdcSpatial::protect_quadtree(x, max_zoom)
+ } else {
+ x <- sdcSpatial::protect_quadtree(x, max_zoom)
+ }
+ x
+}
+
+wrap_sm <- function(x, bw = res(x$value)[1], enlarge = TRUE) {
+
+ if(enlarge) {
+ pad <- 3 * (bw / res(x$value)[1])
+ x$value <- extend(x$value, c(pad, pad))
+ x <- sdcSpatial::protect_smooth(x, bw)
+ } else {
+ x <- sdcSpatial::protect_smooth(bw)
+ }
+ x
+}
+
The idea behind these wrappers is to add some empty cells around the borders of the focus areas in order to capture distribution mass that ‘spills out’. This mass can then be taken into account later-on (esp. for KWD computation).
+
The necessary enlargement depends on the protection mechanism and can be inferred by looking into the functions sdcSpatial::protect_quadtree() and sdcSpatial::protect_smooth().
+
For the quadtree method we need to enlarge the grid by \[2^{\text{max_zoom} - 1} - 1\] rows / columns on each side. For Gaussian smoothing it is \[3 \cdot \frac{\text{bandwidth}}{\text{cellwidth}}\] rows / columns on each side.
+
We call the (wrapped) protection algorithms separately for each focus area (4 methods times 4 areas = 16 outcomes).
+
par(mfrow = c(2, 2))
+
+# removal
+fa1_rm <- remove_sensitive(fa1)
+fa2_rm <- remove_sensitive(fa2)
+fa3_rm <- remove_sensitive(fa3)
+fa4_rm <- remove_sensitive(fa4)
+
+# quad tree I
+fa1_qt1 <- wrap_qt(fa1, max_zoom = 2)
+fa2_qt1 <- wrap_qt(fa2, max_zoom = 2)
+fa3_qt1 <- wrap_qt(fa3, max_zoom = 2)
+fa4_qt1 <- wrap_qt(fa4, max_zoom = 2)
+
+# quad tree II
+fa1_qt2 <- wrap_qt(fa1, max_zoom = 3)
+fa2_qt2 <- wrap_qt(fa2, max_zoom = 3)
+fa3_qt2 <- wrap_qt(fa3, max_zoom = 3)
+fa4_qt2 <- wrap_qt(fa4, max_zoom = 3)
+
+# smoothing
+fa1_sm <- wrap_sm(fa1, bw = cellsize)
+fa2_sm <- wrap_sm(fa2, bw = cellsize)
+fa3_sm <- wrap_sm(fa3, bw = cellsize)
+fa4_sm <- wrap_sm(fa4, bw = cellsize)
+
+
+# show results
+
+plot(fa1_rm, "count", sensitive = FALSE, col = pal, main = "Ruhr valley (rm)")
+plot(fa2_rm, "count", sensitive = FALSE, col = pal, main = "MZ + WI (rm)")
+plot(fa3_rm, "count", sensitive = FALSE, col = pal, main = "Strelasund (rm)")
+plot(fa4_rm, "count", sensitive = FALSE, col = pal, main = "Allgäu (rm)")
+
+
plot(fa1_qt1, "count", sensitive = FALSE, col = pal, main = "Ruhr valley (qt1)")
+plot(fa2_qt1, "count", sensitive = FALSE, col = pal, main = "MZ + WI (qt1)")
+plot(fa3_qt1, "count", sensitive = FALSE, col = pal, main = "Strelasund (qt1)")
+plot(fa4_qt1, "count", sensitive = FALSE, col = pal, main = "Allgäu (qt1)")
+
+
plot(fa1_qt2, "count", sensitive = FALSE, col = pal, main = "Ruhr valley (qt2)")
+plot(fa2_qt2, "count", sensitive = FALSE, col = pal, main = "MZ + WI (qt2)")
+plot(fa3_qt2, "count", sensitive = FALSE, col = pal, main = "Strelasund (qt2)")
+plot(fa4_qt2, "count", sensitive = FALSE, col = pal, main = "Allgäu (qt2)")
+
+
plot(fa1_sm, "count", sensitive = FALSE, col = pal, main = "Ruhr valley (sm)")
+plot(fa2_sm, "count", sensitive = FALSE, col = pal, main = "MZ + WI (sm)")
+plot(fa3_sm, "count", sensitive = FALSE, col = pal, main = "Strelasund (sm)")
+plot(fa4_sm, "count", sensitive = FALSE, col = pal, main = "Allgäu (sm)")
+
+
par(mfrow = c(1, 1))
+
Note that for qt1, qt2 and sm the mass spills slightly beyond the original borders. We can show that with the exception of the removal method, all mass has been preserved:
+
mass_rm <- c(sum(getValues(fa1_rm$value$count), na.rm = TRUE),
+ sum(getValues(fa2_rm$value$count), na.rm = TRUE),
+ sum(getValues(fa3_rm$value$count), na.rm = TRUE),
+ sum(getValues(fa4_rm$value$count), na.rm = TRUE))
+
+mass_qt1 <- c(sum(getValues(fa1_qt1$value$count), na.rm = TRUE),
+ sum(getValues(fa2_qt1$value$count), na.rm = TRUE),
+ sum(getValues(fa3_qt1$value$count), na.rm = TRUE),
+ sum(getValues(fa4_qt1$value$count), na.rm = TRUE))
+
+mass_qt2 <- c(sum(getValues(fa1_qt2$value$count), na.rm = TRUE),
+ sum(getValues(fa2_qt2$value$count), na.rm = TRUE),
+ sum(getValues(fa3_qt2$value$count), na.rm = TRUE),
+ sum(getValues(fa4_qt2$value$count), na.rm = TRUE))
+
+mass_sm <- c(sum(getValues(fa1_sm$value$count), na.rm = TRUE),
+ sum(getValues(fa2_sm$value$count), na.rm = TRUE),
+ sum(getValues(fa3_sm$value$count), na.rm = TRUE),
+ sum(getValues(fa4_sm$value$count), na.rm = TRUE))
+
+masses <- cbind(p_fa, mass_rm, mass_qt1, mass_qt2, mass_sm)
+masses <- cbind(masses, round(masses[, 2:5] - p_fa))
+colnames(masses) <- c("original", "removal", "quad tree I", "quad tree II", "smoothing",
+ "diff_removal", "diff_qt1", "diff_qt2", "diff_sm")
+rownames(masses) <- c("fa1", "fa2", "fa3", "fa4")
+masses
+
## original removal quad tree I quad tree II smoothing diff_removal diff_qt1
+## fa1 1745062 1744773 1745062 1745062 1745062 -289 0
+## fa2 558055 557830 558055 558055 558055 -225 0
+## fa3 107863 107175 107863 107863 107863 -688 0
+## fa4 106283 105116 106283 106283 106283 -1167 0
+## diff_qt2 diff_sm
+## fa1 0 0
+## fa2 0 0
+## fa3 0 0
+## fa4 0 0
+
We can now calculate risk measures for the focus areas.
+
# --- share of cells at risk ---
+# initial
+sens <- c(sensitivity_score(fa1), sensitivity_score(fa2), sensitivity_score(fa3), sensitivity_score(fa4))
+# after SDC
+sens_rm <- c(sensitivity_score(fa1_rm), sensitivity_score(fa2_rm), sensitivity_score(fa3_rm), sensitivity_score(fa4_rm))
+sens_qt1 <- c(sensitivity_score(fa1_qt1), sensitivity_score(fa2_qt1), sensitivity_score(fa3_qt1), sensitivity_score(fa4_qt1))
+sens_qt2 <- c(sensitivity_score(fa1_qt2), sensitivity_score(fa2_qt2), sensitivity_score(fa3_qt2), sensitivity_score(fa4_qt2))
+sens_sm <- c(sensitivity_score(fa1_sm), sensitivity_score(fa2_sm), sensitivity_score(fa3_sm), sensitivity_score(fa4_sm))
+
+risk_cell <- cbind(sens, sens_rm, sens_qt1, sens_qt2, sens_sm)
+colnames(risk_cell) <- c("cell_initial", "cell_removal", "cell_quad tree I", "cell_quad tree II", "cell_smoothing")
+rownames(risk_cell) <- c("fa1", "fa2", "fa3", "fa4")
+round(risk_cell, 3)
+
## cell_initial cell_removal cell_quad tree I cell_quad tree II cell_smoothing
+## fa1 0.033 0 0.007 0 0
+## fa2 0.077 0 0.089 0 0
+## fa3 0.186 0 0.142 0 0
+## fa4 0.213 0 0.077 0 0
+
# --- share of pop. at risk ---
+# initial
+sensp <- c(sum(fa1$value$count[is_sensitive(fa1)]), sum(fa2$value$count[is_sensitive(fa2)]),
+ sum(fa3$value$count[is_sensitive(fa3)]), sum(fa4$value$count[is_sensitive(fa4)]))
+# after SDC
+sensp_rm <- c(sum(fa1_rm$value$count[is_sensitive(fa1_rm)]), sum(fa2_rm$value$count[is_sensitive(fa2_rm)]),
+ sum(fa3_rm$value$count[is_sensitive(fa3_rm)]), sum(fa4_rm$value$count[is_sensitive(fa4_rm)]))
+sensp_qt1 <- c(sum(fa1_qt1$value$count[is_sensitive(fa1_qt1)]), sum(fa2_qt1$value$count[is_sensitive(fa2_qt1)]),
+ sum(fa3_qt1$value$count[is_sensitive(fa3_qt1)]), sum(fa4_qt1$value$count[is_sensitive(fa4_qt1)]))
+sensp_qt2 <- c(sum(fa1_qt2$value$count[is_sensitive(fa1_qt2)]), sum(fa2_qt2$value$count[is_sensitive(fa2_qt2)]),
+ sum(fa3_qt2$value$count[is_sensitive(fa3_qt2)]), sum(fa4_qt2$value$count[is_sensitive(fa4_qt2)]))
+sensp_sm <- c(sum(fa1_sm$value$count[is_sensitive(fa1_sm)]), sum(fa2_sm$value$count[is_sensitive(fa2_sm)]),
+ sum(fa3_sm$value$count[is_sensitive(fa3_sm)]), sum(fa4_sm$value$count[is_sensitive(fa4_sm)]))
+risk_pop <- cbind(sensp/p_fa, sensp_rm/p_fa, sensp_qt1/p_fa, sensp_qt2/p_fa, sensp_sm/p_fa)
+colnames(risk_pop) <- c("pop_initial", "pop_removal", "pop_quad tree I", "pop_quad tree II", "pop_smoothing")
+rownames(risk_pop) <- c("fa1", "fa2", "fa3", "fa4")
+risk_pop
+
## pop_initial pop_removal pop_quad tree I pop_quad tree II pop_smoothing
+## fa1 0.0001656102 0 8.309160e-06 0 0
+## fa2 0.0004031861 0 1.334994e-04 0 0
+## fa3 0.0063784616 0 1.678055e-03 0 0
+## fa4 0.0109801191 0 1.232558e-03 0 0
+
Next, we calculate the normalized Hellinger distances.
+
In accordance with the definition, mass mismatch does not bother the metric, so we consider only the original focus areas without enlargement around the borders. We therefore include a step where if the focus area was enlarged by a wrapper function, it gets cropped to its original size first.
+
getHD <- function(x, orig, value = "count") {
+
+ # resize protected raster if it was enlarged
+ if(extent(x$value) > extent(orig$value)) {
+ x$value <- crop(x$value, extent(orig$value))
+ }
+
+ # extract value raster of interest
+ r_x <- x$value[[value]]
+ r_o <- orig$value[[value]]
+
+ # extract cell values
+ v_x <- raster::getValues(r_x)
+ v_o <- raster::getValues(r_o)
+
+ # set NAs to 0
+ v_x[is.na(v_x)] <- 0
+ v_o[is.na(v_o)] <- 0
+
+ # rescale
+ v_x <- v_x / sum(v_x)
+ v_o <- v_o / sum(v_o)
+
+ HD <- sqrt(sum((sqrt(v_x) - sqrt(v_o))^2)) * (1/sqrt(2))
+ HD
+}
+
+hd_rm <- c(getHD(fa1_rm, fa1), getHD(fa2_rm, fa2), getHD(fa3_rm, fa3), getHD(fa4_rm, fa4))
+hd_qt1 <- c(getHD(fa1_qt1, fa1), getHD(fa2_qt1, fa2), getHD(fa3_qt1, fa3), getHD(fa4_qt1, fa4))
+hd_qt2 <- c(getHD(fa1_qt2, fa1), getHD(fa2_qt2, fa2), getHD(fa3_qt2, fa3), getHD(fa4_qt2, fa4))
+hd_sm <- c(getHD(fa1_sm, fa1), getHD(fa2_sm, fa2), getHD(fa3_sm, fa3), getHD(fa4_sm, fa4))
+
+hd <- cbind(hd_rm, hd_qt1, hd_qt2, hd_sm)
+colnames(hd) <- c("removal", "quad tree I", "quad tree II", "smoothing")
+rownames(hd) <- c("fa1", "fa2", "fa3", "fa4")
+round(hd, 3)
+
## removal quad tree I quad tree II smoothing
+## fa1 0.009 0.070 0.083 0.280
+## fa2 0.014 0.114 0.213 0.366
+## fa3 0.057 0.151 0.205 0.451
+## fa4 0.074 0.167 0.226 0.416
+
+
+
KWD for focus areas
+
KWD requires that both maps have the same mass and extent. By enlarging the focus areas slightly the mass was kept for all methods except removal (see above).
+
In order to account for the slightly larger extent, the ground truth is also slightly enlarged before the comparison. Since all the padding is made up of zero values, it does not distort the measure. This takes care of quadtree and smoothing.
+
getKWD <- function(x, orig, value = "count", balance = "none", plot_map = FALSE, ...) {
+
+ # resize original raster to balance mass for qt/sm
+ if(extent(x$value) > extent(orig$value)) {
+ orig$value <- extend(orig$value, extent(x$value))
+ }
+
+ if(balance == "frame") {
+ x$value <- extend(x$value, c(1, 1), value = 0)
+ orig$value <- extend(orig$value, c(1, 1), value = 0)
+ }
+
+ # extract value raster of interest
+ r_x <- x$value[[value]]
+ r_o <- orig$value[[value]]
+
+ # extract cell values
+ v_x <- raster::getValues(r_x)
+ v_o <- raster::getValues(r_o)
+
+ # set NAs to 0
+ v_x[is.na(v_x)] <- 0
+ v_o[is.na(v_o)] <- 0
+
+ xy <- raster::xyFromCell(r_x, 1:ncell(r_x))
+
+ # balance mass for removal
+ if(balance %in% c("pop_ctrd", "frame")) {
+
+ mdiff <- sum(v_o) - sum(v_x) # mass difference must be positive
+
+ if(balance == "frame") {
+ frameclls <- unique(c(raster::cellFromCol(r_x, c(1, ncol(r_x))),
+ raster::cellFromRow(r_x, c(1, nrow(r_x)))))
+ r_x[frameclls] <- 1/(length(frameclls)) * mdiff
+ }
+
+ if(balance == "pop_ctrd") {
+
+ # get population weighted centroid
+ pop_ctrd <- c(stats::weighted.mean(xy[, 1], v_x),
+ stats::weighted.mean(xy[, 2], v_x))
+ mcell <- cellFromXY(r_x, pop_ctrd)
+ r_x[mcell] <- ifelse(is.na(r_x[mcell]), mdiff, r_x[mcell] + mdiff)
+ }
+
+ # extract new re-distributed mass
+ v_x <- raster::getValues(r_x)
+ v_x[is.na(v_x)] <- 0
+ }
+
+ if(plot_map) {
+ plot(r_x, col = pal)
+ if(balance == "frame") {
+ fr <- extent(r_x, 2, ncol(r_x) - 1, 2, nrow(r_x) - 1)
+ } else if(balance == "pop_ctrd") {
+ fr <- c(xFromCell(r_x, mcell) - res(r_x)[1]/2,
+ xFromCell(r_x, mcell) + res(r_x)[1]/2,
+ yFromCell(r_x, mcell) - res(r_x)[2]/2,
+ yFromCell(r_x, mcell) + res(r_x)[2]/2)
+ }
+ if(balance %in% c("frame", "pop_ctrd")) {
+ rect(xleft = fr[1], xright = fr[2], ybottom = fr[3], ytop = fr[4],
+ border = "red")
+ }
+ }
+
+ KWD <- compareOneToOne(Coordinates = xy, Weights = cbind(v_x, v_o), ...)
+
+ KWD$distance
+}
+
For removal, we need to make assumptions on how to balance mass. Two versions of manually balancing are considered:
+
+The missing mass gets added in a frame around the focus area. Therefore, missing mass in the visual center of the map gets punished more.
The missing mass gets added in a single cell at the population-weighted centroid of the map. The intuition behind this is that this point should be a kind of ‘best guess’ as to where it might be located.
+
# Visualizing the two approaches:
+
+par(mfrow = c(2, 2))
+
+getKWD(fa1_rm, fa1, "count", balance = "frame", plot_map = TRUE)
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.002546098
+
getKWD(fa2_rm, fa2, "count", balance = "frame", plot_map = TRUE)
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.00379919
+
getKWD(fa3_rm, fa3, "count", balance = "frame", plot_map = TRUE)
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.1285745
+
getKWD(fa4_rm, fa4, "count", balance = "frame", plot_map = TRUE)
+
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.1314477
+
getKWD(fa1_rm, fa1, "count", balance = "pop_ctrd", plot_map = TRUE)
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.004185462
+
getKWD(fa2_rm, fa2, "count", balance = "pop_ctrd", plot_map = TRUE)
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.006253455
+
getKWD(fa3_rm, fa3, "count", balance = "pop_ctrd", plot_map = TRUE)
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.1740762
+
getKWD(fa4_rm, fa4, "count", balance = "pop_ctrd", plot_map = TRUE)
+
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
## [1] 0.2377131
+
par(mfrow = c(1, 1))
+
We compute now the KWD for all methods and focus areas. For removal, both ideas about mass balancing are compared.
+
kwd_rm_f <- c(getKWD(fa1_rm, fa1, "count", balance = "frame"),
+ getKWD(fa2_rm, fa2, "count", balance = "frame"),
+ getKWD(fa3_rm, fa3, "count", balance = "frame"),
+ getKWD(fa4_rm, fa4, "count", balance = "frame"))
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_rm_c <- c(getKWD(fa1_rm, fa1, "count", balance = "pop_ctrd"),
+ getKWD(fa2_rm, fa2, "count", balance = "pop_ctrd"),
+ getKWD(fa3_rm, fa3, "count", balance = "pop_ctrd"),
+ getKWD(fa4_rm, fa4, "count", balance = "pop_ctrd"))
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_qt1 <- c(getKWD(fa1_qt1, fa1, "count"), getKWD(fa2_qt1, fa2, "count"),
+ getKWD(fa3_qt1, fa3, "count"), getKWD(fa4_qt1, fa4, "count"))
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_qt2 <- c(getKWD(fa1_qt2, fa1, "count"), getKWD(fa2_qt2, fa2, "count"),
+ getKWD(fa3_qt2, fa3, "count"), getKWD(fa4_qt2, fa4, "count"))
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_sm <- c(getKWD(fa1_sm, fa1, "count"), getKWD(fa2_sm, fa2, "count"),
+ getKWD(fa3_sm, fa3, "count"), getKWD(fa4_sm, fa4, "count"))
+
## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+## CompareOneToOne, Solution method: APPROX
+## WARNING: the Xs input coordinates are not consecutives integers.
+## WARNING: the Ys input coordinates are not consecutives integers.
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd <- cbind(kwd_rm_f, kwd_rm_c, kwd_qt1, kwd_qt2, kwd_sm)
+colnames(kwd) <- c("removal (frame)", "removal (centroid)", "quad tree I", "quad tree II", "smoothing")
+rownames(kwd) <- c("fa1", "fa2", "fa3", "fa4")
+round(kwd, 3)
+
## removal (frame) removal (centroid) quad tree I quad tree II smoothing
+## fa1 0.003 0.004 0.012 0.018 0.381
+## fa2 0.004 0.006 0.030 0.144 0.493
+## fa3 0.129 0.174 0.053 0.119 0.626
+## fa4 0.131 0.238 0.067 0.150 0.624
+
The results we get are now a lot more subtle.
+
Removal is not universally worse than other methods (according to KWD), but actually very good in the more densely populated areas, but beaten by quad tree in rural focus areas.
+
+
+
Appendix:How to normalize KWD?
+
Some considerations on how to normalize KWD from the focusArea() function. Suppose we have a total count of 400. We move 10 units of mass by exactly one cell.
+
wdth <- 11
+m <- matrix(data = 1, nrow = wdth, ncol = wdth)
+
+set.seed(100)
+m[sample(1:(11^2), 50)] <- 6
+m[6, 6] <- m[6, 6] + 29
+
+r_exa <- raster(m, xmn = 0, xmx = wdth, ymn = 0, ymx = wdth)
+
+r_exa2 <- r_exa
+r_exa2[6, 6] <- r_exa2[6, 6] - 10
+r_exa2[6, 7] <- r_exa2[6, 7] + 10
+
+par(mfrow = c(1, 2))
+plot(r_exa, breaks = seq(5, 35, 5), col = terrain.colors(7),
+ main = paste("count total:", sum(getValues(r_exa))))
+plot(r_exa2, breaks = seq(5, 35, 5), col = terrain.colors(7),
+ main = paste("count total:", sum(getValues(r_exa2))))
+
+
par(mfrow = c(1, 1))
+
We know the KWD before normalization must be 10, afterwards 0.025.
+
xy <- ceiling(raster::xyFromCell(r_exa, 1:ncell(r_exa)))
+
+kwd_test <- compareOneToOne(Coordinates = xy,
+ Weights = cbind(getValues(r_exa), getValues(r_exa2)),
+ method = "exact")$distance
+
## CompareOneToOne, Solution method: EXACT
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_fa_test <- focusArea(Coordinates = xy,
+ Weights = cbind(getValues(r_exa), getValues(r_exa2)),
+ x = 6, y = 6, radius = 3, area = "linf",
+ method = "exact")$distance
+
## FocusArea, Solution method: APPROX
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
tst_results <- c(kwd_test, kwd_fa_test, kwd_fa_test / sum(getValues(r_exa2)))
+names(tst_results) <- c("full map", "focus area", "focus area / sum(full map)")
+tst_results
+
## full map focus area
+## 0.025 10.000
+## focus area / sum(full map)
+## 0.025
+
Suppose now we select our focus area from a smaller super area (with fewer total counts). In order to get to the same results we need to calibrate.
+
r_exa_new <- crop(r_exa, extent(r_exa, 2, wdth - 1, 2, wdth - 1))
+r_exa2_new <- crop(r_exa2, extent(r_exa2, 2, wdth - 1, 2, wdth - 1))
+
+fa_factor <- sum(getValues(r_exa)) / sum(getValues(r_exa_new))
+
+xy <- ceiling(raster::xyFromCell(r_exa_new, 1:ncell(r_exa_new)))
+
+kwd_test <- compareOneToOne(Coordinates = xy,
+ Weights = cbind(getValues(r_exa_new), getValues(r_exa2_new)),
+ method = "exact")$distance
+
## CompareOneToOne, Solution method: EXACT
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
kwd_fa_test <- focusArea(Coordinates = xy,
+ Weights = cbind(getValues(r_exa_new), getValues(r_exa2_new)),
+ x = 5, y = 5, radius = 3, area = "linf",
+ method = "exact")$distance
+
## FocusArea, Solution method: APPROX
+## INFO: Recoding the input coordinates to consecutive integers.
+## INFO: change <timelimit> to 14400.000000
+## INFO: change <verbosity> to silent
+## INFO: change <opt_tolerance> to 0.000001
+
tst_results <- c(kwd_test, kwd_test / fa_factor,
+ kwd_fa_test / sum(r_exa2_new[(5-3):(5+3), (5-3):(5+3)]),
+ (kwd_fa_test / sum(getValues(r_exa2_new)) / fa_factor))
+names(tst_results) <- c("cropped area", "cropped area (calib.)", "focus area", "focus area (calib.)")
+tst_results
+
## cropped area cropped area (calib.) focus area
+## 0.03571429 0.02500000 0.05617978
+## focus area (calib.)
+## 0.02500000
+
+
+
+
+
+