| Title: | Determining and Evaluating High-Risk Zones |
|---|---|
| Description: | Functions for determining and evaluating high-risk zones and simulating and thinning point process data, as described in 'Determining high risk zones using point process methodology - Realization by building an R package' Seibold (2012) <http://highriskzone.r-forge.r-project.org/Bachelorarbeit.pdf> and 'Determining high-risk zones for unexploded World War II bombs by using point process methodology', Mahling et al. (2013) <doi:10.1111/j.1467-9876.2012.01055.x>. |
| Authors: | Heidi Seibold <[email protected]>, Monia Mahling <[email protected]>, Sebastian Linne <[email protected]>, Felix Guenther <[email protected]> |
| Maintainer: | Rickmer Schulte <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.4.9 |
| Built: | 2024-12-21 06:44:23 UTC |
| Source: | CRAN |
The package highriskzone provides tools to determine and evaluate high-risk zones of unobserved events by using point process methodology.
Heidi Seibold [email protected], Monia Mahling [email protected] Sebastian Linne [email protected] Felix Guenther [email protected] Maintainer: Felix Guenther [email protected]
Monia Mahling, Michael Hoehle & Helmut Kuechenhoff (2013),
Determining high-risk zones for unexploded World War II bombs by using point process methodology.
Journal of the Royal Statistical Society, Series C 62(2), 181-199.
Monia Mahling (2013),
Determining high-risk zones by using spatial point process methodology.
Ph.D. thesis, Cuvillier Verlag Goettingen,
available online: http://edoc.ub.uni-muenchen.de/15886/
Heidi Seibold (2012), Determining high risk zones using point process methodology - Realization by building an R package. Bachelor Thesis, Ludwig Maximilian University of Munich.
Simulation-based iterative procedure to correct for possible bias with respect to the failure probability alpha
bootcor( ppdata, cutoff, numit = 1000, tol = 0.02, nxprob = 0.1, intens = NULL, covmatrix = NULL, simulate = "intens", radiusClust = NULL, clustering = 5, verbose = TRUE )bootcor( ppdata, cutoff, numit = 1000, tol = 0.02, nxprob = 0.1, intens = NULL, covmatrix = NULL, simulate = "intens", radiusClust = NULL, clustering = 5, verbose = TRUE )
ppdata |
Observed spatial point process of class ppp. |
cutoff |
Desired failure probability alpha, which is the probability of having unobserved events outside the high-risk zone. |
numit |
Number of iterations to perform (per tested value for cutoff). Default value is 1000. |
tol |
Tolerance: acceptable difference between the desired failure probability and the fraction of high-risk zones not covering all events. Default value is 0.02. |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
intens |
(optional) estimated intensity of the observed process (object of class "im",
see |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only meaningful if no intensity is given. If not given, it will be estimated. |
simulate |
The type of simulation, can be one of |
radiusClust |
(optional) radius of the circles around the parent points in which the cluster
points are located. Only used for |
clustering |
a value >= 1 which describes the amount of clustering; the
adjusted estimated intensity of the observed pattern is divided by
this value; it also is the parameter of the Poisson distribution
for the number of points per cluster. Only used for |
verbose |
logical. Should information on tested values/progress be printed? |
For a desired failure probability alpha, the corresponding parameter which is to use
when determining a high-risk zone is found in an iterative procedure. The simulation procedure
is the same as in eval_method. In every iteration,
the number of high-risk zones with at least one unobserved event located outside is
compared with the desired failure probability. If necessary, the value of cutoff is
increased or decreased. The final value alphastar can than be used in
det_hrz.
If there are restriction areas in the observation window, use bootcor_restr
instead.
An object of class bootcorr, which consists of a list of the final value for alpha (alphastar)
and a data.frame course containing information on the simulation course, e.g. the tested values.
Monia Mahling, Michael H?hle & Helmut K?chenhoff (2013), Determining high-risk zones for unexploded World War II bombs by using point process methodology. Journal of the Royal Statistical Society, Series C 62(2), 181-199.
Monia Mahling (2013), Determining high-risk zones by using spatial point process methodology. Ph.D. thesis, Cuvillier Verlag G?ttingen, available online: http://edoc.ub.uni-muenchen.de/15886/ Chapter 6
det_hrz, eval_method, bootcor_restr
## Not run: data(craterB) set.seed(4321) bc <- bootcor(ppdata=craterB, cutoff=0.2, numit=100, tol=0.02, nxprob=0.1) bc summary(bc) plot(bc) hrzbc <- det_hrz(craterB, type = "intens", criterion = "indirect", cutoff = bc$alphastar, nxprob = 0.1) ## End(Not run)## Not run: data(craterB) set.seed(4321) bc <- bootcor(ppdata=craterB, cutoff=0.2, numit=100, tol=0.02, nxprob=0.1) bc summary(bc) plot(bc) hrzbc <- det_hrz(craterB, type = "intens", criterion = "indirect", cutoff = bc$alphastar, nxprob = 0.1) ## End(Not run)
Simulation-based iterative procedure to correct for possible bias with respect to the failure probability alpha
bootcor_restr( ppdata, cutoff, numit = 100, tol = 0.001, nxprob = 0.1, hole = NULL, obsprobimage = NULL, intens = NULL, covmatrix = NULL, simulate = "intens", radiusClust = NULL, clustering = 5, verbose = TRUE )bootcor_restr( ppdata, cutoff, numit = 100, tol = 0.001, nxprob = 0.1, hole = NULL, obsprobimage = NULL, intens = NULL, covmatrix = NULL, simulate = "intens", radiusClust = NULL, clustering = 5, verbose = TRUE )
ppdata |
Observed spatial point process of class ppp. |
cutoff |
Desired failure probability alpha, which is the probability of having unobserved events outside the high-risk zone. |
numit |
Number of iterations to perform (per tested value for cutoff). Default value is 1000. |
tol |
Tolerance: acceptable difference between the desired failure probability and the fraction of high-risk zones not covering all events. Default value is 0.02. |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
hole |
(optional) an object of class |
obsprobimage |
(optional) an object of class |
intens |
(optional) estimated intensity of the observed process (object of class "im",
see |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only meaningful if no intensity is given. If not given, it will be estimated. |
simulate |
The type of simulation, can be one of |
radiusClust |
(optional) radius of the circles around the parent points in which the cluster
points are located. Only used for |
clustering |
a value >= 1 which describes the amount of clustering; the
adjusted estimated intensity of the observed pattern is divided by
this value; it also is the parameter of the Poisson distribution
for the number of points per cluster. Only used for |
verbose |
logical. Should information on tested values/progress be printed? |
For a desired failure probability alpha, the corresponding parameter which is to use
when determining a high-risk zone is found in an iterative procedure. The simulation procedure
is the same as in eval_method. In every iteration,
the number of high-risk zones with at least one unobserved event located outside is
compared with the desired failure probability. If necessary, the value of cutoff is
increased or decreased. The final value alphastar can than be used in
det_hrz.
The function offers the possibility to take into account so-called restriction areas. This is relevant in
situations where the observed point pattern ppdata is incomplete. If it is known that no observations
can be made in a certain area (for example because of water expanses),
this can be accounted for by integrating a hole in the observation window.
The shape and location of the hole is given by hole. Holes are
part of the resulting high-risk zone.
Another approach consists in weighting the observed events with their reciprocal observation probability when
estimating the intensity. To do so, the observation probability can be specified by using
obsprobsimage (an image of the observation probability). Note that the
observation probability may vary in space.
For further information, see Mahling (2013), Appendix A (References).
If there are no restriction areas in the observation window, bootcor
can be used instead.
An object of class bootcorr, which consists of a list of the final value for alpha (alphastar)
and a data.frame course containing information on the simulation course, e.g. the tested values.
Monia Mahling, Michael H?hle & Helmut K?chenhoff (2013), Determining high-risk zones for unexploded World War II bombs by using point process methodology. Journal of the Royal Statistical Society, Series C 62(2), 181-199.
Monia Mahling (2013), Determining high-risk zones by using spatial point process methodology. Ph.D. thesis, Cuvillier Verlag G?ttingen, available online: http://edoc.ub.uni-muenchen.de/15886/ Chapter 6 and Appendix A
data(craterA) set.seed(4321) # define restriction area restrwin <- spatstat.geom::owin(xrange = craterA$window$xrange, yrange = craterA$window$yrange, poly = list(x = c(1500, 1500, 2000, 2000), y = c(2000, 1500, 1500, 2000))) # create image of observation probability (30% inside restriction area) wim <- spatstat.geom::as.im(craterA$window, value = 1) rim <- spatstat.geom::as.im(restrwin, xy = list(x = wim$xcol, y = wim$yrow)) rim$v[is.na(rim$v)] <- 0 oim1 <- spatstat.geom::eval.im(wim - 0.7 * rim) ## Not run: # perform bootstrap correction bc1 <- bootcor_restr(ppdata=craterA, cutoff=0.4, numit=100, tol=0.02, obsprobimage=oim1, nxprob=0.1) bc1 summary(bc1) plot(bc1) # determine high-risk zone by weighting the observations hrzi1 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = bc1$alphastar, hole=NULL, obsprobs=NULL, obsprobimage=oim1, nxprob = 0.1) # perform bootstrap correction set.seed(4321) bc2 <- bootcor_restr(ppdata=craterA, cutoff=0.4, numit=100, tol=0.02, hole=restrwin, nxprob=0.1) bc2 summary(bc2) plot(bc2) # determine high-risk zone by accounting for a hole hrzi2 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = bc2$alphastar, hole=restrwin, obsprobs=NULL, obsprobimage=NULL, nxprob = 0.1) ## End(Not run)data(craterA) set.seed(4321) # define restriction area restrwin <- spatstat.geom::owin(xrange = craterA$window$xrange, yrange = craterA$window$yrange, poly = list(x = c(1500, 1500, 2000, 2000), y = c(2000, 1500, 1500, 2000))) # create image of observation probability (30% inside restriction area) wim <- spatstat.geom::as.im(craterA$window, value = 1) rim <- spatstat.geom::as.im(restrwin, xy = list(x = wim$xcol, y = wim$yrow)) rim$v[is.na(rim$v)] <- 0 oim1 <- spatstat.geom::eval.im(wim - 0.7 * rim) ## Not run: # perform bootstrap correction bc1 <- bootcor_restr(ppdata=craterA, cutoff=0.4, numit=100, tol=0.02, obsprobimage=oim1, nxprob=0.1) bc1 summary(bc1) plot(bc1) # determine high-risk zone by weighting the observations hrzi1 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = bc1$alphastar, hole=NULL, obsprobs=NULL, obsprobimage=oim1, nxprob = 0.1) # perform bootstrap correction set.seed(4321) bc2 <- bootcor_restr(ppdata=craterA, cutoff=0.4, numit=100, tol=0.02, hole=restrwin, nxprob=0.1) bc2 summary(bc2) plot(bc2) # determine high-risk zone by accounting for a hole hrzi2 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = bc2$alphastar, hole=restrwin, obsprobs=NULL, obsprobimage=NULL, nxprob = 0.1) ## End(Not run)
Simulation-based iterative procedure to correct for possible bias with respect to the failure probability alpha
bootcorr( ppdata, cutoff, numit = 1000, tol = 0.02, nxprob = 0.1, intens = NULL, covmatrix = NULL, simulate = "intens", radiusClust = NULL, clustering = 5, verbose = TRUE )bootcorr( ppdata, cutoff, numit = 1000, tol = 0.02, nxprob = 0.1, intens = NULL, covmatrix = NULL, simulate = "intens", radiusClust = NULL, clustering = 5, verbose = TRUE )
ppdata |
Observed spatial point process of class ppp. |
cutoff |
Desired failure probability alpha, which is the probability of having unobserved events outside the high-risk zone. |
numit |
Number of iterations to perform (per tested value for cutoff). Default value is 1000. |
tol |
Tolerance: acceptable difference between the desired failure probability and the fraction of high-risk zones not covering all events. Default value is 0.02. |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
intens |
(optional) estimated intensity of the observed process (object of class "im",
see |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only meaningful if no intensity is given. If not given, it will be estimated. |
simulate |
The type of simulation, can be one of |
radiusClust |
(optional) radius of the circles around the parent points in which the cluster
points are located. Only used for |
clustering |
a value >= 1 which describes the amount of clustering; the
adjusted estimated intensity of the observed pattern is divided by
this value; it also is the parameter of the Poisson distribution
for the number of points per cluster. Only used for |
verbose |
logical. Should information on tested values/progress be printed? |
For a desired failure probability alpha, the corresponding parameter which is to use
when determining a high-risk zone is found in an iterative procedure. The simulation procedure
is the same as in eval_method. In every iteration,
the number of high-risk zones with at least one unobserved event located outside is
compared with the desired failure probability. If necessary, the value of cutoff is
increased or decreased. The final value alphastar can than be used in
det_hrz.
If there are restriction areas in the observation window, use bootcor_restr
instead.
An object of class bootcorr, which consists of a list of the final value for alpha (alphastar)
and a data.frame course containing information on the simulation course, e.g. the tested values.
Monia Mahling, Michael H?hle & Helmut K?chenhoff (2013), Determining high-risk zones for unexploded World War II bombs by using point process methodology. Journal of the Royal Statistical Society, Series C 62(2), 181-199.
Monia Mahling (2013), Determining high-risk zones by using spatial point process methodology. Ph.D. thesis, Cuvillier Verlag G?ttingen, available online: http://edoc.ub.uni-muenchen.de/15886/ Chapter 6
det_hrz, eval_method, bootcor_restr
## Not run: data(craterB) set.seed(4321) bc <- bootcor(ppdata=craterB, cutoff=0.2, numit=100, tol=0.02, nxprob=0.1) bc summary(bc) plot(bc) hrzbc <- det_hrz(craterB, type = "intens", criterion = "indirect", cutoff = bc$alphastar, nxprob = 0.1) ## End(Not run)## Not run: data(craterB) set.seed(4321) bc <- bootcor(ppdata=craterB, cutoff=0.2, numit=100, tol=0.02, nxprob=0.1) bc summary(bc) plot(bc) hrzbc <- det_hrz(craterB, type = "intens", criterion = "indirect", cutoff = bc$alphastar, nxprob = 0.1) ## End(Not run)
For each argument it is checked if it is of a correct value or class.
check_det_hrz_input( ppdata, type, criterion, cutoff, distancemap, intens, nxprob, covmatrix )check_det_hrz_input( ppdata, type, criterion, cutoff, distancemap, intens, nxprob, covmatrix )
ppdata |
Observed spatial point process of class ppp. |
type |
Method to use, can be one of |
criterion |
criterion to limit the high-risk zone, can be one of
|
cutoff |
Value of criterion (area, radius, quantile, alpha or threshold). Depending on criterion and type: If criterion = "direct" and type = "intens", cutoff is the maximum intensity of unexploded bombs outside the risk zone. If type = "dist" instead, cutoff is the radius of the circle around each exploded bomb. "If criterion = "indirect", cutoff is the quantile for the quantile-based method and the failure probability alpha for the intensity-base method. If criterion = "area", cutoff is the area the high-risk zone should have. |
distancemap |
(optional) distance map: distance of every pixel to the nearest observation
of the point pattern; only needed for |
intens |
(optional) estimated intensity of the observed process (object of class "im"),
only needed for type="intens". If not given,
it will be estimated using |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only needed for
|
For each argument it is checked if it is of a correct value or class.
check_det_hrz_restr_input( ppdata, type, criterion, cutoff, hole, integratehole, obsprobs, obsprobimage, distancemap, intens, nxprob, covmatrix, returnintens )check_det_hrz_restr_input( ppdata, type, criterion, cutoff, hole, integratehole, obsprobs, obsprobimage, distancemap, intens, nxprob, covmatrix, returnintens )
ppdata |
Observed spatial point process of class ppp. |
type |
Method to use, can be one of |
criterion |
criterion to limit the high-risk zone, can be one of
|
cutoff |
Value of criterion (area, radius, quantile, alpha or threshold). Depending on criterion and type. |
hole |
(optional) an object of class |
integratehole |
Should the |
obsprobs |
(optional) Vector of observation probabilities associated with the observations contained in |
obsprobimage |
(optional) an object of class |
distancemap |
(optional) distance map: distance of every pixel to the nearest observation
of the point pattern; only needed for |
intens |
(optional) estimated intensity of the observed process (object of class "im",
see |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only needed for
|
returnintens |
Should the image of the estimated intensity be returned? Defaults to |
Bomb crater Point Pattern
data(craterA)data(craterA)
An object of class "ppp" representing a point pattern of bomb craters. The Cartesian coordinates are in meters.
See ppp.object for details of the format of a point pattern object.
Bomb crater Point Pattern
data(craterB)data(craterB)
An object of class "ppp" representing a point pattern of bomb craters. The Cartesian coordinates are in meters.
See ppp.object for details of the format of a point pattern object.
This function is used for the intensity-based method. It determines the probability to have at least one unobserved event outside the high-risk zone. A Poisson distribution is used for the number of unobserved events in a certain area or field. Used in functions det_threshold, det_thresholdfromarea.
det_alpha(intens, threshold, nxprob = 0.1)det_alpha(intens, threshold, nxprob = 0.1)
intens |
estimated intensity of the observed process (object of class "im", see |
threshold |
threshold c: The high-risk zone is the field in which the estimated intensity exceeds this value. |
nxprob |
probability of having unobserved events |
value of alpha
Determination of failure probability within evaluation area
det_alpha_eval_ar(intens, eval_ar, threshold, nxprob = 0.1)det_alpha_eval_ar(intens, eval_ar, threshold, nxprob = 0.1)
intens |
estimated intensity |
eval_ar |
evaluation area |
threshold |
given threshold |
nxprob |
constant probability of non-explosion |
This function is used for the intensity-based method. Calculation of the area of the high-risk zone given the observation window, the intensity matrix and the threshold c. Used in function det_thresholdfromarea.
det_area(win, intensmatrix, threshold)det_area(win, intensmatrix, threshold)
win |
observation window |
intensmatrix |
matrix of the estimated intensity of the observed process ( |
threshold |
threshold c: The high-risk zone is the field in which the estimated intensity exceeds this value |
A numerical value giving the area of the high-risk zone.
This function is used for the intensity-based method with a hole restriction area. Calculation of the area of the high-risk zone given the observation window, the intensity matrix, the threshold c and a hole. Used in function det_thresholdfromarea_hole.
det_area_hole(win, intensmatrix, threshold, hole, integratehole = TRUE)det_area_hole(win, intensmatrix, threshold, hole, integratehole = TRUE)
win |
observation window |
intensmatrix |
matrix of the estimated intensity of the observed process ( |
threshold |
threshold c: The high-risk zone is the field in which the estimated intensity exceeds this value |
hole |
specified hole |
integratehole |
Should the |
A numerical value giving the area of the high-risk zone.
det_guard_width determines the necessary width of a guard region in which
the existence of additional observed bomb craters could change a intensity based estimated
highriskzone within the evaluation area of interest.
Within the evaluation area, the high risk zone consists of all points at which the estimated
intensity of unexploded bombs exceeds a certain, specified or estimated threshold c. At a given
point s, the intensity of unexploded bombs is given by the sum of all evaluated bivariate normal kernels
centered at the observed bomb craters multiplied by a constant nxprob/1-nxprob.
If the estimated intensity of unexploaded bombs is zero at a point at the boarder of the evaluation area
an additional observation outside the area could lift the intensity only above the determined threshold if the
distance to the boarder is small enough so that the density of the normal kernel (which is centered at the additional
observation) is bigger than the threshold at the boarder (assuming that the estimated kernel doesn't change due to the
additional observation). The function returns the biggest distance in which it is possible that the density of
the bivariate normal kernel of the intensity of the supplied highriskzone exceeds thresh_const times the threshold
of the highriskzone. If thresh_const is set to 1, the guard region is the smallest region with constant width around
the evaluation area in which a single additional observation could (but not necessarily does) increase the
highriskzone within the evaluation area at a point at the boarder if the intensity of unexploaded bombs was zero at this point before.
If the intensity was >0 at a point at the boarder of the evaluation area, or more than 1 additional observations
are found nearby outside of the evaluation area, the highriskzone within the evaluation area could already expand
by addditional observations with a bigger distance from the boarder. This can be considered by setting thresh_const < 1,
which intuitively means that 1/thresh_const crater observation at the same point could expand the highriskzone within
the evaluation area in the direction of the additional observations, or that a point the boarder becomes part of the highriskzone
by the observation of a single additional crater if the intensity at this point was thresh_cont times the highriskzone threshold
based on all crater observations within the evaluation area.
det_guard_width(highriskzone, thresh_const = 0.5)det_guard_width(highriskzone, thresh_const = 0.5)
highriskzone |
the estimated highriskzone for the evaluation area |
thresh_const |
the constant multiplied with the determined threshold, 0 < thresh_const < 1. |
For more infos on the construction of guard zones see Mahling (2013, Appendix B, Approach 2)
The constant width of the guard region.
## change npixel to 1000 to obtain nicer plots spatstat.geom::spatstat.options(npixel=100) data(craterA) # reduce number of observations for faster computation thin.craterA <- craterA[1:50] hrzi1 <- det_hrz(thin.craterA, type = "intens", criterion = "area", cutoff = 100000, nxprob = 0.1) det_guard_width(hrzi1, thresh_const = .25)## change npixel to 1000 to obtain nicer plots spatstat.geom::spatstat.options(npixel=100) data(craterA) # reduce number of observations for faster computation thin.craterA <- craterA[1:50] hrzi1 <- det_hrz(thin.craterA, type = "intens", criterion = "area", cutoff = 100000, nxprob = 0.1) det_guard_width(hrzi1, thresh_const = .25)
det_hrz determines the high-risk zone through the method of fixed radius
(type = "dist" and criterion = "direct"), the quantile-based method (type = "dist" and
criterion = "area"/"indirect") and the intensity-based method (type = "intens").
det_hrz( ppdata, type, criterion, cutoff, distancemap = NULL, intens = NULL, nxprob = 0.1, covmatrix = NULL )det_hrz( ppdata, type, criterion, cutoff, distancemap = NULL, intens = NULL, nxprob = 0.1, covmatrix = NULL )
ppdata |
Observed spatial point process of class ppp. |
type |
Method to use, can be one of |
criterion |
criterion to limit the high-risk zone, can be one of
|
cutoff |
Value of criterion (area, radius, quantile, alpha or threshold). Depending on criterion and type: If criterion = "direct" and type = "intens", cutoff is the maximum intensity of unexploded bombs outside the risk zone. If type = "dist" instead, cutoff is the radius of the circle around each exploded bomb. "If criterion = "indirect", cutoff is the quantile for the quantile-based method and the failure probability alpha for the intensity-base method. If criterion = "area", cutoff is the area the high-risk zone should have. |
distancemap |
(optional) distance map: distance of every pixel to the nearest observation
of the point pattern; only needed for |
intens |
(optional) estimated intensity of the observed process (object of class "im"),
only needed for type="intens". If not given,
it will be estimated using |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only needed for
|
There are different methods implemented to determine a high-risk zone.
In this method, the high-risk zone is determined by drawing a circle around each
observed event with a fixed radius. This method will be used when type = "dist"
and criterion = "direct". cutoff then is the radius.
This method is a development of the above. Here the radius is not fixed. It uses
the distance of every observed event to the nearest other event, which is calculated by the
nearest-neighbour distance. The radius is assessed by the p-quantile of the empirical
distribution function of the nearest-neighbour distance. This method will be used when
type = "dist" and criterion = "indirect" or "area". If
criterion = "indirect", then cutoff is the quantile that should be used.
If criterion = "area" then cutoff is the area that the high-risk zone
has to have at the end and from that the quantile/the radii are determined. When the
calculation is done via the area, it can not really be classified to the quantile-based
method. It is rather a third "distance-based" method.
The first step of this method is to estimate the intensity of the observed events.
Based on the estimated intensity and the specified probability of unobserved bombs nxprob
it is possible to estimate the intensity of unobserved/unexploded bombs.
The high-risk zone is then the area in which the estimated intensity of unexploded bombs exceeds a
certain value. This value is called threshold c.
The method will be used when type = "intens". There are three different ways to
construct a high-risk zone:
Fixing the threshold c: criterion = "direct"
Fixing the area of the high-risk zone: criterion = "area"
Fixing the failure probability alpha, which is the probability of having
unobserved events outside the high-risk zone: criterion = "indirect"
Here, the point process is assumed to be an inhomogeneous Poisson process.
For further information see Mahling et al. (2013) (References).
If there are restriction areas in the observation window, use det_hrz_restr
instead. For estimation of intensity based highrikszones with a bigger observation area than area of interest
(evaluation area) use det_hrz_eval_ar.
An object of class "highriskzone", which is a list of
typehrz, criterion, cutoff, nxprob
|
see arguments |
zone |
Determined high-risk zone: Object of class "owin" based on a binary mask.
See |
threshold |
determined threshold. If type = "dist" and criterion = "direct" it is the specified radius. If criterion = "indirect" or "area" the determined radius used to construct a risk zone fulfilling the specified criterion and cutoff. If type = "dist" it is the specified or calculated threshold c, the maximum intensitiy of unexploded bombs outside the risk zone. |
calccutoff |
determined cutoff-value. For type="dist" and criterion="area", this is the quantile of the nearest-neighbour distance. For type="intens" and criterion="area" or "direct", it is the failure probability alpha. For all other criterions it is NA. |
covmatrix |
If not given (and |
Monia Mahling, Michael Hoehle & Helmut Kuechenhoff (2013), Determining high-risk zones for unexploded World War II bombs by using point process methodology. Journal of the Royal Statistical Society, Series C 62(2), 181-199.
Monia Mahling (2013), Determining high-risk zones by using spatial point process methodology. Ph.D. thesis, Cuvillier Verlag Goettingen, available online: http://edoc.ub.uni-muenchen.de/15886/
distmap, eval.im, owin,
eval_method, det_hrz_restr
data(craterA) ## change npixel to 1000 to obtain nicer plots spatstat.geom::spatstat.options(npixel=100) ## type: dist hrzd1 <- det_hrz(craterA, type = "dist", criterion = "area", cutoff = 1000000, nxprob = 0.1) hrzd2 <- det_hrz(craterA, type = "dist", criterion = "indirect", cutoff = 0.9, nxprob = 0.1) hrzd3 <- det_hrz(craterA, type = "dist", criterion = "direct", cutoff = 100, nxprob = 0.1) op <- par(mfrow = c(2, 2)) plot(craterA) plot(hrzd1, zonecol = 2, win = craterA$window, plotwindow = TRUE) plot(hrzd2, zonecol = 3, win = craterA$window, plotwindow = TRUE) plot(hrzd3, zonecol = 4, win = craterA$window, plotwindow = TRUE) par(op) ## Not run: # or first calculate the distancemap and use it: distm <- distmap(craterA) hrzd <- det_hrz(craterA, type = "dist", criterion = "direct", cutoff = 100, distancemap = distm, nxprob = 0.1) ## End(Not run) ## type: intens # reduce number of observations for faster computation thin.craterA <- craterA[1:10] hrzi1 <- det_hrz(thin.craterA, type = "intens", criterion = "area", cutoff = 100000, nxprob = 0.1) plot(hrzi1) plot(thin.craterA, add = TRUE) plot(thin.craterA$window, add = TRUE) ## Not run: hrzi2 <- det_hrz(craterA, type = "intens", criterion = "indirect", cutoff = 0.1, nxprob = 0.1) hrzi3 <- det_hrz(craterA, type = "intens", criterion = "direct", cutoff = 0.0001, nxprob = 0.1) plot(hrzi2) plot(hrzi3) ## End(Not run) ## More detailed examples on http://highriskzone.r-forge.r-project.org/data(craterA) ## change npixel to 1000 to obtain nicer plots spatstat.geom::spatstat.options(npixel=100) ## type: dist hrzd1 <- det_hrz(craterA, type = "dist", criterion = "area", cutoff = 1000000, nxprob = 0.1) hrzd2 <- det_hrz(craterA, type = "dist", criterion = "indirect", cutoff = 0.9, nxprob = 0.1) hrzd3 <- det_hrz(craterA, type = "dist", criterion = "direct", cutoff = 100, nxprob = 0.1) op <- par(mfrow = c(2, 2)) plot(craterA) plot(hrzd1, zonecol = 2, win = craterA$window, plotwindow = TRUE) plot(hrzd2, zonecol = 3, win = craterA$window, plotwindow = TRUE) plot(hrzd3, zonecol = 4, win = craterA$window, plotwindow = TRUE) par(op) ## Not run: # or first calculate the distancemap and use it: distm <- distmap(craterA) hrzd <- det_hrz(craterA, type = "dist", criterion = "direct", cutoff = 100, distancemap = distm, nxprob = 0.1) ## End(Not run) ## type: intens # reduce number of observations for faster computation thin.craterA <- craterA[1:10] hrzi1 <- det_hrz(thin.craterA, type = "intens", criterion = "area", cutoff = 100000, nxprob = 0.1) plot(hrzi1) plot(thin.craterA, add = TRUE) plot(thin.craterA$window, add = TRUE) ## Not run: hrzi2 <- det_hrz(craterA, type = "intens", criterion = "indirect", cutoff = 0.1, nxprob = 0.1) hrzi3 <- det_hrz(craterA, type = "intens", criterion = "direct", cutoff = 0.0001, nxprob = 0.1) plot(hrzi2) plot(hrzi3) ## End(Not run) ## More detailed examples on http://highriskzone.r-forge.r-project.org/
det_hrz_eval_ar determines intensity based highriskzones if bomb crater observations are available
for a bigger area than the area of main interest (evaluation area).
All observations are used for intensity estimation, the highriskzone is however constructed only in the
evaluation area. Either based on specifying a failure probability alpha that indicates the probability of
unobserved bombs outside the highriskzone but inside the evaluation area of interest (and not in the
overall observation area) (criterion = "indirect"), or by specifying the threshold (maximum intensity of non-
exploded bombs outside the) highriskzone directly and intersecting the resulting hrz with the
evaluation area (criterion = "direct").
det_hrz_eval_ar( ppdata, eval_ar, criterion = c("indirect", "direct"), cutoff, intens = NULL, nxprob = 0.1, covmatrix = NULL )det_hrz_eval_ar( ppdata, eval_ar, criterion = c("indirect", "direct"), cutoff, intens = NULL, nxprob = 0.1, covmatrix = NULL )
ppdata |
Observed spatial point process of class ppp in the observation area. |
eval_ar |
area of interest specified via an object of class owin |
criterion |
criterion to limit the high-risk zone, can be |
cutoff |
Value of criterion (alpha or threshold) |
intens |
(optional) estimated intensity of the observed process (object of class "im") in (bigger)
observation area, if not given, it will be estimated using |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only needed for
|
An object of class "highriskzone"
set.seed(12412) spatstat.geom::spatstat.options(npixel=300) data(craterB) # reduce number of observations for faster computation thin.craterB <- craterB[sample(1:craterB$n, 40)] # define evaluation area of interest eval.ar <- spatstat.geom::owin(xrange = c(0, 1900), yrange = c(0, 3400), poly = matrix(c(250,250, 1200,1000,250,1000), byrow = TRUE, ncol = 2)) hrzi1 <- det_hrz_eval_ar(thin.craterB, eval_ar = eval.ar, criterion = "direct", cutoff = 3e-6, nxprob = .2) plot(hrzi1) plot(thin.craterB, add = TRUE) plot(eval.ar, add = TRUE) plot(craterB$window, add = TRUE)set.seed(12412) spatstat.geom::spatstat.options(npixel=300) data(craterB) # reduce number of observations for faster computation thin.craterB <- craterB[sample(1:craterB$n, 40)] # define evaluation area of interest eval.ar <- spatstat.geom::owin(xrange = c(0, 1900), yrange = c(0, 3400), poly = matrix(c(250,250, 1200,1000,250,1000), byrow = TRUE, ncol = 2)) hrzi1 <- det_hrz_eval_ar(thin.craterB, eval_ar = eval.ar, criterion = "direct", cutoff = 3e-6, nxprob = .2) plot(hrzi1) plot(thin.craterB, add = TRUE) plot(eval.ar, add = TRUE) plot(craterB$window, add = TRUE)
det_hrz_restr determines the high-risk zone through the method of fixed radius
(type = "dist" and criterion = "direct"), the quantile-based method (type = "dist" and
criterion = "area"/"indirect") and the intensity-based method (type = "intens").
Restriction areas can be taken into account.
det_hrz_restr( ppdata, type, criterion, cutoff, hole = NULL, integratehole = TRUE, obsprobs = NULL, obsprobimage = NULL, distancemap = NULL, intens = NULL, nxprob = 0.1, covmatrix = NULL, returnintens = TRUE )det_hrz_restr( ppdata, type, criterion, cutoff, hole = NULL, integratehole = TRUE, obsprobs = NULL, obsprobimage = NULL, distancemap = NULL, intens = NULL, nxprob = 0.1, covmatrix = NULL, returnintens = TRUE )
ppdata |
Observed spatial point process of class ppp. |
type |
Method to use, can be one of |
criterion |
criterion to limit the high-risk zone, can be one of
|
cutoff |
Value of criterion (area, radius, quantile, alpha or threshold). Depending on criterion and type. |
hole |
(optional) an object of class |
integratehole |
Should the |
obsprobs |
(optional) Vector of observation probabilities associated with the observations contained in |
obsprobimage |
(optional) an object of class |
distancemap |
(optional) distance map: distance of every pixel to the nearest observation
of the point pattern; only needed for |
intens |
(optional) estimated intensity of the observed process (object of class "im",
see |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only needed for
|
returnintens |
Should the image of the estimated intensity be returned? Defaults to |
Used in functions eval_method, sim_clintens, sim_intens.
This function contains the same functionalities as det_hrz.
In addition, it offers the possibility to take into account so-called restriction areas. This is relevant in
situations where the observed point pattern ppdata is incomplete. If it is known that no observations
can be made in a certain area (for example because of water expanses),
this can be accounted for by integrating a hole in the observation window.
The shape and location of the hole is given by hole, whereas integratehole is used to state
whether the hole is to become part of the resulting high-risk zone.
This may also be a reasonable approach if only few observations could be made in a certain area.
Another approach consists in weighting the observed events with their reciprocal observation probability when
estimating the intensity. To do so, the observation probability can be specified by using obsprobs (value of the
observation probability for each event) or obsprobsimage (image of the observation probability). Note that the
observation probability may vary in space.
If there are no restriction areas in the observation window, det_hrz
can be used instead.
Note that for criterion = "area", cutoff specifies the area of the high-risk zone outside the hole. If
integratehole = TRUE, the area of the resulting high-risk zone will exceed cutoff.
For further information, Mahling et al. (2013) and Mahling (2013), Chapters 4 and 8 and Appendix A (References).
An object of class "highriskzone", which is a list of
typehrz, criterion, cutoff, nxprob
|
see arguments |
zone |
Determined high-risk zone: Object of class "owin" based on a binary mask.
See |
threshold |
determined threshold. If type = "dist" and criterion = "direct" it is the specified radius. If criterion = "indirect" or "area" the determined radius used to construct a risk zone fulfilling the specified criterion and cutoff. If type = "dist" it is the specified or calculated threshold c, the maximum intensitiy of unexploded bombs outside the risk zone. |
calccutoff |
determined cutoff-value. For type="dist" and criterion="area", this is the quantile of the nearest-neighbour distance. For type="intens" and criterion="area" or "direct", it is the failure probability alpha. For all other criterions it is NA. |
covmatrix |
If not given (and |
estint |
Estimated intensity. See |
set.seed(1211515) data(craterA) #change npixel = 100 to 1000 to get a nicer picture spatstat.geom::spatstat.options(npixel=100) # reduce number of observations for faster computation craterA <- craterA[sample(1:craterA$n, 150)] # define restriction area restrwin <- spatstat.geom::owin(xrange=craterA$window$xrange, yrange=craterA$window$yrange, poly=list(x=c(1500, 1500, 2000, 2000), y=c(2000, 1500, 1500, 2000))) # create image of observation probability (30% inside restriction area) wim <- spatstat.geom::as.im(craterA$window, value=1) rim <- spatstat.geom::as.im(restrwin, xy=list(x=wim$xcol, y=wim$yrow)) rim$v[is.na(rim$v)] <- 0 oim1 <- spatstat.geom::eval.im(wim - 0.7 * rim) # determine high-risk zone by weighting the observations hrzi1 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = 0.4, hole=NULL, obsprobs=NULL, obsprobimage=oim1, nxprob = 0.1) # determine high-risk zone by accounting for a hole hrzi2 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = 0.4, hole=restrwin, obsprobs=NULL, obsprobimage=NULL, nxprob = 0.1)set.seed(1211515) data(craterA) #change npixel = 100 to 1000 to get a nicer picture spatstat.geom::spatstat.options(npixel=100) # reduce number of observations for faster computation craterA <- craterA[sample(1:craterA$n, 150)] # define restriction area restrwin <- spatstat.geom::owin(xrange=craterA$window$xrange, yrange=craterA$window$yrange, poly=list(x=c(1500, 1500, 2000, 2000), y=c(2000, 1500, 1500, 2000))) # create image of observation probability (30% inside restriction area) wim <- spatstat.geom::as.im(craterA$window, value=1) rim <- spatstat.geom::as.im(restrwin, xy=list(x=wim$xcol, y=wim$yrow)) rim$v[is.na(rim$v)] <- 0 oim1 <- spatstat.geom::eval.im(wim - 0.7 * rim) # determine high-risk zone by weighting the observations hrzi1 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = 0.4, hole=NULL, obsprobs=NULL, obsprobimage=oim1, nxprob = 0.1) # determine high-risk zone by accounting for a hole hrzi2 <- det_hrz_restr(ppdata=craterA, type = "intens", criterion = "indirect", cutoff = 0.4, hole=restrwin, obsprobs=NULL, obsprobimage=NULL, nxprob = 0.1)
Used in function det_radius.
det_nnarea(cutoffval, distancemap, win)det_nnarea(cutoffval, distancemap, win)
cutoffval |
distance used as radius of the discs |
distancemap |
distance map (object of class "im", see |
win |
observation window of class owin |
A numerical value giving the area of the window.
Used in function sim_nsppp.
det_nsintens(ppdata, radius)det_nsintens(ppdata, radius)
ppdata |
observed point pattern whose estimated intensity (adjusted for thinning and divided by "clustering") is used for simulating the parent process |
radius |
radius of the circles around the parent points in which the cluster points are located |
A pixel image (object of class "im"). See density.ppp.
density.ppp, boundingbox,
owin, Hscv
Used in function bootcor_restr.
det_nsintens_restr(ppdata, radius, weights)det_nsintens_restr(ppdata, radius, weights)
ppdata |
observed point pattern whose estimated intensity (adjusted for thinning and divided by "clustering") is used for simulating the parent process |
radius |
radius of the circles around the parent points in which the cluster points are located |
weights |
Vector of observation probabilities associated with the observations contained in |
A pixel image (object of class "im"). See density.ppp.
density.ppp, boundingbox,
owin, Hscv
Used in function det_hrz.
det_radius(ppdata, distancemap, areahrz, win)det_radius(ppdata, distancemap, areahrz, win)
ppdata |
observed spatial point pattern of class ppp. |
distancemap |
distance map (object of class "im", see |
areahrz |
given area of the high-risk zone |
win |
observation window of class owin |
A list of
cutoffdist |
quantile of the nearest-neighbour distance |
thresh |
distance |
The high-risk zone is the field in which the estimated intensity exceeds the threshold c, which is determined here, having the failure probability alpha. This function is for the intensity-based method. Used in function det_hrz.
det_threshold(intens, alpha = 1e-05, nxprob = 0.1)det_threshold(intens, alpha = 1e-05, nxprob = 0.1)
intens |
estimated intensity of the observed process (object of class "im", see |
alpha |
failure probability: probability to have at least one unobserved event outside the high-risk zone |
nxprob |
probability of having unobserved events |
value of the threshold c
Determination of necessary threshold to keep alpha in evaluation area
det_threshold_eval_ar(intens, eval_ar, alpha = 1e-05, nxprob = 0.1)det_threshold_eval_ar(intens, eval_ar, alpha = 1e-05, nxprob = 0.1)
intens |
estimated intensity |
eval_ar |
evaluation area |
alpha |
desired failure probability in eval area |
nxprob |
constant probability of non-explosion |
This function is used for the intensity-based method. Used in function det_hrz.
det_thresholdfromarea(intens, areahrz, win, nxprob = 0.1)det_thresholdfromarea(intens, areahrz, win, nxprob = 0.1)
intens |
estimated intensity of the observed process (object of class "im", see |
areahrz |
area of the high-risk zone |
win |
observation window |
nxprob |
probability of having unbserved events |
A list of
threshold |
Value of the threshold c. The high-risk zone is the field in which the estimated intensity exceeds this value |
calccutoff |
failure probability alpha for given area; probability to have at least unobserved event outside the high-risk zone |
This function is used for the intensity-based method. Used in function det_hrz_restr.
det_thresholdfromarea_rest( intens, areahrz, win, nxprob = 0.1, hole = hole, integratehole = TRUE )det_thresholdfromarea_rest( intens, areahrz, win, nxprob = 0.1, hole = hole, integratehole = TRUE )
intens |
estimated intensity of the observed process (object of class "im", see |
areahrz |
area of the high-risk zone |
win |
observation window |
nxprob |
probability of having unbserved events |
hole |
an object of class |
integratehole |
Should the |
A list of
threshold |
Value of the threshold c. The high-risk zone is the field in which the estimated intensity exceeds this value |
calccutoff |
failure probability alpha for given area; probability to have at least unobserved event outside the high-risk zone |
Estimates the intensity of the point pattern by a kernel method
(See density.ppp).
est_intens(ppdata, covmatrix = NULL, weights = NULL)est_intens(ppdata, covmatrix = NULL, weights = NULL)
ppdata |
data of class ppp |
covmatrix |
(Optional) Covariance matrix of the kernel of a normal distribution |
weights |
(Optional) vector of weights attached to each observation |
A list of
intensest |
Estimated intensity (object of class "im", see |
covmatrix |
Covariance matrix. If |
data(craterA) #change npixel = 50 to 1000 to get a nicer picture spatstat.geom::spatstat.options(npixel=50) # use only ten observations for fast computation thin.craterA <- craterA[1:10] int <- est_intens(thin.craterA) # Plot estimated intensity plot(int$intensest, main = "pixel image of intensity") plot(craterA$window, main = "contour plot of intensity") contour(int$intensest, add =TRUE)data(craterA) #change npixel = 50 to 1000 to get a nicer picture spatstat.geom::spatstat.options(npixel=50) # use only ten observations for fast computation thin.craterA <- craterA[1:10] int <- est_intens(thin.craterA) # Plot estimated intensity plot(int$intensest, main = "pixel image of intensity") plot(craterA$window, main = "contour plot of intensity") contour(int$intensest, add =TRUE)
Estimates the intensity of the point pattern by using the SPDE method from r-INLA.
est_intens_spde( coords, win = NULL, npixel = 50, fine_mesh = FALSE, mesh = NULL, weights = NULL, alpha = 2, ... )est_intens_spde( coords, win = NULL, npixel = 50, fine_mesh = FALSE, mesh = NULL, weights = NULL, alpha = 2, ... )
coords |
ppp object or matrix with x and y coordinates of the observed bombs |
win |
observation window, either of class owin or a matrix with the x and y coordinates of the boundary, not neccessary if coords is a ppp object |
npixel |
number of pixel per dimension (see |
fine_mesh |
logical, if FALSE a coarse mesh will be created, if TRUE a fine mesh will be created, only used if argument mesh is NULL |
mesh |
(optional) a predefined mesh for the spde model |
weights |
(optional) integration weights for the spde model, only used if argument mesh is NULL |
alpha |
(optional) alpha value for the spde model, only used if argument spde is NULL |
... |
additional arguments for the construction of the spde model (see INLA/inla.spde2.matern documentation) |
A list of
intensest |
Pixel image with the estimated intensities of the random field. |
mesh |
The mesh. |
## Not run: data(craterA) est_spde <- est_intens_spde(coords=craterA) image.plot(list(x=est_spde$intensest$xcol, y=est_spde$intensest$yrow, z=log(t(est_spde$intensest$v))), main="estimated logarithmic intensity") points(craterA) ## End(Not run)## Not run: data(craterA) est_spde <- est_intens_spde(coords=craterA) image.plot(list(x=est_spde$intensest$xcol, y=est_spde$intensest$yrow, z=log(t(est_spde$intensest$v))), main="estimated logarithmic intensity") points(craterA) ## End(Not run)
Estimates the intensity of the point pattern by a kernel method
(See density.ppp).
est_intens_weight(ppdata, covmatrix = NULL, weights = NULL)est_intens_weight(ppdata, covmatrix = NULL, weights = NULL)
ppdata |
data of class ppp |
covmatrix |
(Optional) Covariance matrix of the kernel of a normal distribution |
weights |
(Optional) vector of weights attached to each observation |
A list of
intensest |
Estimated intensity (object of class "im", see |
covmatrix |
Covariance matrix. If |
data(craterA) #change npixel = 50 to 1000 to get a nicer picture spatstat.geom::spatstat.options(npixel=50) # use only ten observations for fast computation thin.craterA <- craterA[1:10] # weight first 5 observations twice weights <- c(rep(2, 5), rep(1, 5)) int <- est_intens_weight(thin.craterA, weights = weights) plot(int$intensest, main = "pixel image of intensity") plot(craterA$window, main = "contour plot of intensity") contour(int$intensest, add =TRUE)data(craterA) #change npixel = 50 to 1000 to get a nicer picture spatstat.geom::spatstat.options(npixel=50) # use only ten observations for fast computation thin.craterA <- craterA[1:10] # weight first 5 observations twice weights <- c(rep(2, 5), rep(1, 5)) int <- est_intens_weight(thin.craterA, weights = weights) plot(int$intensest, main = "pixel image of intensity") plot(craterA$window, main = "contour plot of intensity") contour(int$intensest, add =TRUE)
Evaluation of the high-risk zone, which is only possible with simulated or thinned data or if the locations of the unobserved events have been revealed..
eval_hrz(hrz, unobspp, obspp = NULL)eval_hrz(hrz, unobspp, obspp = NULL)
hrz |
High-risk zone of class owin based on a binary mask (see |
unobspp |
Unobserved spatial point process |
obspp |
Observed spatial point process |
An object of class "hrzeval", which is a list of
numbermiss |
number of unobserved events outside the high-risk zone |
numberunobserved |
number of events in the unobserved point pattern |
missingfrac |
fraction of unobserved events outside the high-risk zone (numbermiss/numberunobserved) |
arearegion |
area of the high-risk zone |
numberobs |
number of events in the observed point pattern |
out |
subset of the unobserved events, which are outside the high-risk zone |
insd |
subset of the unobserved events, which are inside the high-risk zone |
data(craterB) # thin data set.seed(100) thdata <- thin(craterB, nxprob=0.1) # determine hrz for the "observed events" hrz <- det_hrz(thdata$observed, type = "dist", criterion = "area", cutoff = 1500000, nxprob = 0.1) # evaluate the hrz evaluation <- eval_hrz(hrz = hrz$zone, unobspp = thdata$unobserved, obspp = thdata$observed) evaluation$missingfrac op <- par(mar=c(1, 4, 1, 6) , xpd=TRUE) plot(evaluation, hrz = hrz, obspp = thdata$observed, plothrz = TRUE, plotobs = TRUE, insidecol = "magenta", outsidecol = "magenta", obscol = "blue", insidepch = 1, outsidepch = 19, main = "Evaluation visualized") legend(2400, 2456.4061, c("observed", "unobs inside", "unobs outside"), col = c("blue", "magenta", "magenta"), yjust=1, pch=c(1, 1, 19), cex=0.8) par(op)data(craterB) # thin data set.seed(100) thdata <- thin(craterB, nxprob=0.1) # determine hrz for the "observed events" hrz <- det_hrz(thdata$observed, type = "dist", criterion = "area", cutoff = 1500000, nxprob = 0.1) # evaluate the hrz evaluation <- eval_hrz(hrz = hrz$zone, unobspp = thdata$unobserved, obspp = thdata$observed) evaluation$missingfrac op <- par(mar=c(1, 4, 1, 6) , xpd=TRUE) plot(evaluation, hrz = hrz, obspp = thdata$observed, plothrz = TRUE, plotobs = TRUE, insidecol = "magenta", outsidecol = "magenta", obscol = "blue", insidepch = 1, outsidepch = 19, main = "Evaluation visualized") legend(2400, 2456.4061, c("observed", "unobs inside", "unobs outside"), col = c("blue", "magenta", "magenta"), yjust=1, pch=c(1, 1, 19), cex=0.8) par(op)
Evaluates the performance of the three methods:
Method of fixed radius
Quantile-based method
Intensity-based method
For further details on the methods, see det_hrz or the paper of Mahling et al. (2013)(References).
There are three ways to simulate data for the evaluation.
eval_method( ppdata, type, criterion, cutoff, numit = 100, nxprob = 0.1, distancemap = NULL, intens = NULL, covmatrix = NULL, simulate, radiusClust = NULL, clustering = 5, pbar = TRUE )eval_method( ppdata, type, criterion, cutoff, numit = 100, nxprob = 0.1, distancemap = NULL, intens = NULL, covmatrix = NULL, simulate, radiusClust = NULL, clustering = 5, pbar = TRUE )
ppdata |
Observed spatial point process of class ppp. |
type |
Method to use, can be one of |
criterion |
criterion to limit the high-risk zone, can be one of
|
cutoff |
Value of criterion (area, radius, quantile, alpha or threshold). Depending on criterion and type: If criterion = "direct" and type = "intens", cutoff is the maximum intensity of unexploded bombs outside the risk zone. If type = "dist" instead, cutoff is the radius of the circle around each exploded bomb. "If criterion = "indirect", cutoff is the quantile for the quantile-based method and the failure probability alpha for the intensity-base method. If criterion = "area", cutoff is the area the high-risk zone should have. |
numit |
Number of iterations |
nxprob |
Probability of having unobserved events. Default value is 0.1. |
distancemap |
(optional) distance map: distance of every pixel to the nearest observation
of the point pattern; only needed for |
intens |
(optional) estimated intensity of the observed process (object of class "im"),
only needed for type="intens". If not given,
it will be estimated using |
covmatrix |
(optional) Covariance matrix of the kernel of a normal distribution, only needed for
|
simulate |
The type of simulation, can be one of |
radiusClust |
(Optional) radius of the circles around the parent points in which the cluster
points are located. Only used for |
clustering |
a value >= 1 which describes the amount of clustering; the
adjusted estimated intensity of the observed pattern is divided by
this value; it is also the parameter of the Poisson distribution
for the number of points per cluster. Only used for |
pbar |
logical. Should progress bar be printed? |
The three simulation types are:
Here a given data set is used. The data set is thinned as explained below. Note that this method is very different from the others, since it is using the real data.
Here, an inhomogeneous Poisson process is simulated and then that data is thinned.
Here a Neyman-Scott process is simulated (see sim_nsppp, rNeymanScott)
and this data is then also thinned.
Thinning:
Let be the spatial point process, which is the location of all events and let
be a subset of describing the observed process. The process of unobserved events
then is Z = X \ Y , meaning that and are disjoint and together
forming .
Since is not known, in this function an observed or simulated spatial point pattern
ppdata is taken as the full pattern (which we denote by X') comprising the
observed events Y' as well as the unobserved Z'.
Each event in X' is assigned to one of the two processes Y' or
Z' by drawing independent Bernoulli random numbers.
The resulting process of observed events Y' is used to determine the high-risk zone.
Knowing now the unobserved process, it can be seen how many events are outside and inside the
high-risk zone.
type and criterion may be vectors in this function.
A data.frame with variables
Iteration |
Iterationstep of the result |
Type, Criterion, Cutoff, nxprob
|
see arguments |
threshold |
determined threshold. If criterion="area", it is either the distance (if type="dist") or the threshold c (for type="intens"). If criterion="indirect", it is either the quantile of the nearest-neighbour distance which is used as radius (if type="dist") or the threshold c (for type="intens"). If criterion="direct", it equals the cutoff for both types. |
calccutoff |
determined cutoff-value. For type="dist" and criterion="area", this is the quantile of the nearest-neighbour distance. For type="intens" and criterion="area", it is the failure probability alpha. For all other criterions it is NA. |
covmatrix11, covmatrix12, covmatrix21, covmatrix22
|
values in the covariance matrix. covmatrix11 and covmatrix22 are the diagonal elements (variances). |
numbermiss |
number of unobserved points outside the high-risk zone |
numberunobserved |
number of observations in the unobserved point pattern Z' |
missingfrac |
fraction of unobserved events outside the high-risk zone (numbermiss/numberunobserved) |
arearegion |
area of the high-risk zone |
numberobs |
number of observations in the observed point pattern Y' |
det_hrz, rNeymanScott, thin, sim_nsppp, sim_intens
## Not run: data(craterB) # the input values are mainly the same as in det_hrz, so for more example ideas, # see the documentation of det_hrz. evalm <- eval_method(craterB, type = c("dist", "intens"), criterion = c("area", "area"), cutoff = c(1500000, 1500000), nxprob = 0.1, numit = 10, simulate = "clintens", radiusClust = 300, clustering = 15, pbar = FALSE) evalm_d <- subset(evalm, evalm$Type == "dist") evalm_i <- subset(evalm, evalm$Type == "intens") # pout: fraction of high-risk zones that leave at least one unobserved event uncovered # pmiss: Mean fraction of unobserved events outside the high-risk zone data.frame(pmiss_d = mean(evalm_d$missingfrac), pmiss_i = mean(evalm_i$missingfrac), pout_d = ( sum(evalm_d$numbermiss > 0) / nrow(evalm_d) ), pout_i = ( sum(evalm_i$numbermiss > 0) / nrow(evalm_i) )) ## End(Not run)## Not run: data(craterB) # the input values are mainly the same as in det_hrz, so for more example ideas, # see the documentation of det_hrz. evalm <- eval_method(craterB, type = c("dist", "intens"), criterion = c("area", "area"), cutoff = c(1500000, 1500000), nxprob = 0.1, numit = 10, simulate = "clintens", radiusClust = 300, clustering = 15, pbar = FALSE) evalm_d <- subset(evalm, evalm$Type == "dist") evalm_i <- subset(evalm, evalm$Type == "intens") # pout: fraction of high-risk zones that leave at least one unobserved event uncovered # pmiss: Mean fraction of unobserved events outside the high-risk zone data.frame(pmiss_d = mean(evalm_d$missingfrac), pmiss_i = mean(evalm_i$missingfrac), pout_d = ( sum(evalm_d$numbermiss > 0) / nrow(evalm_d) ), pout_i = ( sum(evalm_i$numbermiss > 0) / nrow(evalm_i) )) ## End(Not run)
Plot a visualization of the bootstrap correction for a high-risk zone. The different values tested for alpha are plotted.
## S3 method for class 'bootcorr' plot(x, ...)## S3 method for class 'bootcorr' plot(x, ...)
x |
bootstrap correction for a high-risk zone (object of class " |
... |
extra arguments passed to the generic |
This is the plot method for the class bootcorr.
plot, print.bootcorr, summary.bootcorr
Plot a high-risk zone.
## S3 method for class 'highriskzone' plot( x, ..., pattern = NULL, win = NULL, plotpattern = FALSE, plotwindow = FALSE, windowcol = "white", usegpclib = FALSE, zonecol = "grey" )## S3 method for class 'highriskzone' plot( x, ..., pattern = NULL, win = NULL, plotpattern = FALSE, plotwindow = FALSE, windowcol = "white", usegpclib = FALSE, zonecol = "grey" )
x |
high-risk zone (object of class " |
... |
extra arguments passed to the generic |
pattern |
spatial point pattern for which the highriskzone was determined. |
win |
observation winodw |
plotpattern |
logical flag; if |
plotwindow |
logical flag; if |
windowcol |
the color used to plot the observation window |
usegpclib |
logical flag; if |
zonecol |
the colour used to plot the high-risk zone. |
This is the plot method for the class highriskzone.
plot, for examples see det_hrz
Plot a visualization of the evaluation of a high-risk zone. At least the observation window and the unobserved events inside and outside the high-risk zone are plotted.
## S3 method for class 'hrzeval' plot( x, ..., hrz = NULL, obspp = NULL, plothrz = FALSE, plotobs = FALSE, windowcol = "white", insidecol = "blue", outsidecol = "red", insidepch = 20, outsidepch = 19, zonecol = "grey", obscol = "black", obspch = 1 )## S3 method for class 'hrzeval' plot( x, ..., hrz = NULL, obspp = NULL, plothrz = FALSE, plotobs = FALSE, windowcol = "white", insidecol = "blue", outsidecol = "red", insidepch = 20, outsidepch = 19, zonecol = "grey", obscol = "black", obspch = 1 )
x |
evaluation of a high-risk zone (object of class " |
... |
extra arguments passed to the generic |
hrz |
(optional) high-risk zone (object of class " |
obspp |
(optional) observed point pattern |
plothrz |
logical flag; should the high-risk zone be plotted? |
plotobs |
logical flag; should the observed point pattern be plotted? |
windowcol |
the color used to plot the observation window |
insidecol |
the color used to plot the unobserved events inside the high-risk zone |
outsidecol |
the color used to plot the unobserved events outside the high-risk zone |
insidepch |
plotting 'character' of the unobserved events inside the high-risk zone,
i.e., symbol to use. This can either be a single character or an integer code for one of
a set of graphics symbols. The full set of S symbols is available with pch=0:18, see
|
outsidepch |
plotting 'character' of the unobserved events outside the high-risk zone |
zonecol |
the color used to plot the high-risk zone |
obscol |
the color used to plot the observed events |
obspch |
plotting 'character' of the observed events |
This is the plot method for the class hrzeval.
plot, eval_hrz, plot.highriskzone
Prints a very brief description of the bootstrap correction for a high-risk zone.
## S3 method for class 'bootcorr' print(x, ...)## S3 method for class 'bootcorr' print(x, ...)
x |
bootstrap correction for of a high-risk zone (object of class " |
... |
ignored |
A very brief description of the bootstrap correction x for a high-risk zone is printed.
This is a method for the generic function print.
Prints a very brief description of a high-risk zone.
## S3 method for class 'highriskzone' print(x, ...)## S3 method for class 'highriskzone' print(x, ...)
x |
high-risk zone (object of class " |
... |
ignored |
A very brief description of the highriskzone x is printed.
This is a method for the generic function print.
Prints a very brief description of the evaluation of a high-risk zone.
## S3 method for class 'hrzeval' print(x, ...)## S3 method for class 'hrzeval' print(x, ...)
x |
evaluation of a high-risk zone (object of class " |
... |
ignored |
A very brief description of the evaluation x of a high-risk zone is printed.
This is a method for the generic function print.
If xwin or ywin is NULL, the observation window will be a rectangular bounding box. Vertices must be listed anticlockwise; no vertex should be repeated. Only needed for data that is not already of class ppp.
read_pppdata(xppp, yppp, xwin = NULL, ywin = NULL, unitname = NULL)read_pppdata(xppp, yppp, xwin = NULL, ywin = NULL, unitname = NULL)
xppp |
Vector of x coordinates of data points |
yppp |
Vector of y coordinates of data points |
xwin |
Vector of x coordinates of the vertices of a polygon circumscribing the observation window |
ywin |
Vector of y coordinates of the vertices of a polygon circumscribing the observation window |
unitname |
Optional. Name of unit of length. Either a single character string, or a vector of two character strings giving the singular and plural forms, respectively. |
An object of class "ppp" describing a point pattern in the two-dimensional plane.
data(craterA) windowA <- data.frame(x = craterA$window$bdry[[1]]$x, y = craterA$window$bdry[[1]]$y) patternA <- data.frame(x = craterA$x, y = craterA$y) str(patternA) str(windowA) crater <- read_pppdata(xppp = patternA$x, yppp = patternA$y, xwin = windowA$x, ywin = windowA$y) craterdata(craterA) windowA <- data.frame(x = craterA$window$bdry[[1]]$x, y = craterA$window$bdry[[1]]$y) patternA <- data.frame(x = craterA$x, y = craterA$y) str(patternA) str(windowA) crater <- read_pppdata(xppp = patternA$x, yppp = patternA$y, xwin = windowA$x, ywin = windowA$y) crater
Generation of a random point pattern using the inhomogeneous Poisson process (if lambda is not constant) and thinning of this data, to obtain "observed" and "unobserved" events.
sim_intens(ppdata, intensSim, nxprob)sim_intens(ppdata, intensSim, nxprob)
ppdata |
Observed spatial point process of class ppp |
intensSim |
Intensity to use for the simulation |
nxprob |
Probability of having unobserved events |
A list of of observed and unobserved point patterns (see thin)
This algorithm generates a realisation of a Neyman-Scott process whose expected number of points equals the number of observations in a given pattern.
sim_nsppp(ppdata, radius, clustering = 5, thinning = 0)sim_nsppp(ppdata, radius, clustering = 5, thinning = 0)
ppdata |
observed point pattern, whose estimated intensity (adjusted for thinning and divided by "clustering") is used for simulating the parent process |
radius |
radius of the circles around the parent points in which the cluster points are located (Maximum radius of a random cluster) |
clustering |
a value larger or equal 1 which describes the amount of clustering; the adjusted estimated intensity of the observed pattern is divided by this value; it is also the parameter of the Poisson distribution for the number of points per cluster |
thinning |
constant thinning probability (in case the observed pattern is a thinned version of a full pattern); usually equal to the probability of having unobserved events |
First, the algorithm generates a Poisson point process (see rpoispp for
details) of parent points with intensity kappa, which is a pixel image
object of class "im" (see im.object).
This pixel image is derived from the observed pattern using density.ppp.
The bandwidth is not chosen in advance.
If only a thinned version of the original pattern has been observed,
this can be taken into account using the parameter thinning.
Usually, not the estimated intensity itself is used for simulating the
parent process, but its values are divided by a constant named "clustering".
Second, each parent point is replaced by a random cluster of points, created
by calling the function runifdisc. Each cluster consists of a Poisson
distributed number of points (with clustering being the expected number of
points in each cluster) which are located in a disc of a given radius.
These clusters are combined to yield a single point pattern which is
then returned as the result.
The estimation of the intensity (on an adequate window) and the
simulation of the Neyman-Scott process are performed seperately,
so the intensity does not need to be reestimated in every iteration.
The resulting process is a Mat?rn process whose parent process is an
inhomogeneous Poisson point process.
The simulated point pattern (an object of class "ppp").
Additionally, some intermediate results of the simulation are returned as
attributes of this point pattern: see rNeymanScott.
rNeymanScott, rThomas,
rMatClust
## Not run: data(craterA) data(craterB) set.seed(100) sim_pp1 <- sim_nsppp(craterA, radius=300, clustering=15, thinning=0.1) sim_pp2 <- sim_nsppp(craterB, radius=300, clustering=15, thinning=0.1) op <- par(mfrow = c(1, 2)) plot(sim_pp1, main = "simulated cluster process 1") plot(sim_pp2, main = "simulated cluster process 2") par(op) ## End(Not run)## Not run: data(craterA) data(craterB) set.seed(100) sim_pp1 <- sim_nsppp(craterA, radius=300, clustering=15, thinning=0.1) sim_pp2 <- sim_nsppp(craterB, radius=300, clustering=15, thinning=0.1) op <- par(mfrow = c(1, 2)) plot(sim_pp1, main = "simulated cluster process 1") plot(sim_pp2, main = "simulated cluster process 2") par(op) ## End(Not run)
Simulation of the Neyman-Scott process. Only applicable if the intensity was estimated
for an appropriately enlarged window.
More details in sim_nsppp.
sim_nsprocess(ppdata, intens, radius, clustering = 5, thinning = 0)sim_nsprocess(ppdata, intens, radius, clustering = 5, thinning = 0)
ppdata |
observed point pattern whose estimated intensity (adjusted for thinning and divided by "clustering") is used for simulating the parent process |
intens |
estimated intensity |
radius |
radius of the circles around the parent points in which the cluster points are located (Maximum radius of a random cluster) |
clustering |
a value larger or equal 1 which describes the amount of clustering; the adjusted estimated intensity of the observed pattern is divided by this value; it is also the parameter of the Poisson distribution for the number of points per cluster |
thinning |
constant thinning probability (in case the observed pattern is a thinned version of a full pattern); usually equal to the probability of having unobserved events |
The simulated point pattern (an object of class "ppp").
Additionally, some intermediate results of the simulation are returned as
attributes of this point pattern: see rNeymanScott.
Prints a useful summary of the bootstrap correction for a high-risk zone.
## S3 method for class 'bootcorr' summary(object, ...)## S3 method for class 'bootcorr' summary(object, ...)
object |
bootstrap correction for a high-risk zone (object of class " |
... |
ignored |
A useful summary of the bootstrap correction x for a high-risk zone is printed.
This is a method for the generic function summary.
summary, print.bootcorr, plot.bootcorr
Prints a useful summary of a high-risk zone.
## S3 method for class 'highriskzone' summary(object, ...)## S3 method for class 'highriskzone' summary(object, ...)
object |
high-risk zone (object of class " |
... |
ignored |
A useful description of the highriskzone object is printed.
This is a method for the generic function summary.
Prints a useful summary of the evaluation of a high-risk zone.
## S3 method for class 'hrzeval' summary(object, ...)## S3 method for class 'hrzeval' summary(object, ...)
object |
evaluation of a high-risk zone (object of class " |
... |
ignored |
A useful description of the hrzeval object is printed.
This is a method for the generic function summary.
The thinning is done by drawing independently from a Bernoulli distribution. This function is needed for functions eval_method, sim_clintens, sim_intens
thin(full, nxprob)thin(full, nxprob)
full |
all observations of the point pattern |
nxprob |
probability of having unobserved events |
A list of observed and unobserved point patterns. Both of class ppp.
data(craterB) thdata <- thin(craterB, nxprob=0.1) thdata plot(thdata$observed); points(thdata$unobserved, col=4)data(craterB) thdata <- thin(craterB, nxprob=0.1) thdata plot(thdata$observed); points(thdata$unobserved, col=4)