# How to determine if point is surrounded using raster processing?

I am attempting to improve a currently extremely cumbersome vector/python process for a natural hazard model. At the moment we have a lengthy script which generates distance/bearing lines from a given point to determine:

1. the type of polygon that it intersects (e.g.. Forest, grass, swamp, etc)
2. the distance to that polygon
3. how many of these lines intersect polygons, to determine how 'surrounded' it is.

There's a lot more involved but that's the gist of it. I am trying to find a way to improve this and am currently stumped on part 3. The idea is to determine if a point is completely surrounded by polygons, within say 200m So in my attached image, I would want point A to be marked as being at higher risk than point B as it is completely surrounded by my polygons. This is repeated for about 13 million points so is not a small task and I would rather have a surface to derive values from, rather than running our script. I'm thinking there's got to be a variation of hydrology tools or cost-path to do this but I can't seem to get my head around it.

Some assistance/suggestions on the best way to go about this would be MUCH appreciated.

• Viewshed is up to the task but is going to need considerable help when applied to 13 million points! Think first about how to weed out points whose surroundings are easy to check, such as points located in regions external to the polygons that have diameters less than (say) 200m: that could rule out "A" but perhaps not "B". "B" would never be ruled out because its viewshed (taking the polygon areas to be very "high" locations blocking the view) extends more than 200m from B's location. – whuber Jan 13 '14 at 1:25
• A good point @whuber. certainly I am able to reduce the total number of points actually processed by proximity and in fact unique lat-longs as well (i'm talking geocoded addresses so apartment blocks can be reduced from 50 points to 1) but I'll still be looking at several million locations. I can also just break everything up into overlapping blocks if necessary. Will investigate viewshed. Thanks! – Loz Jan 13 '14 at 2:27
• Another quick screen is to compute a focal mean of the 0-1 indicator grid of your polygons using an annular neighborhood: at any cells where its value is 1, your polygons occupy the entire perimeter out at the radius, whence they must be surrounded. This is a fast calculation and might weed out the vast majority of your points, depending on where they are and how convoluted your polygons are. You can also speed up initial screening by first resampling your grid to a coarser resolution, such as 25-50 m or thereabouts. – whuber Jan 13 '14 at 15:53
• Another potential processing step, or pre-processing step, would be to pass your points through a rasterized version of your dataset comparing the statistics of a neighborhood around the points. You could either abstract away your 'surrounded' requirement as a statistic of the neighborhood of points, or if 'surrounded' is necessary, you could find the 'easy' points (i.e. a point completely within a risk area) using a raster neighborhood, parse out the 'easy' points from all points, then use vector analysis for the rest of the points. – DPierce Jan 17 '14 at 19:30
• wow my query has certainly generated a lot of interest! Thank you to everyone who's contributed suggestions and comments. I'm going to work my way though them all and respond but they're all going to take some time for me to test. I promise I will respond eventually! – Loz Feb 19 '14 at 2:40

There is a cost-path solution but you will have to code it yourself. Here is what it might look like when applied to every point in the image in the question (coarsened a little to speed up the calculations): The black cells are parts of the surrounding polygons. The colors, ranging from light orange (short) through blue (long), show the maximum distance (out to a maximum of 50 cells) that can be reached by line-of-sight traversal without intercepting the polygon cells. (Any cell outside the extent of this image is treated as part of the polygons.)

Let's discuss an efficient way to do that using a raster representation of the data. In this representation all the "surrounding" polygonal cells will have, say, nonzero values and any cell that can be "seen through" will have a zero value.

### Step 1: Precomputing a neighborhood data structure

You first have to decide what it means for one cell to block another. One of the fairest rules I can find is this: using integral coordinates for rows and columns (and assuming square cells), let us consider which cells might block cell (i,j) from the view at the origin (0,0). I nominate the cell (i',j') which is closest to the line segment connecting (i,j) to (0,0) among all cells whose coordinates differ from i and j by at most 1. Because this does not always yield a unique solution (for instance, with (i,j) = (1,2) both (0,1) and (1,1) will work equally well), some means to resolve ties is needed. It would be nice for this resolution of ties to respect the symmetries of circular neighborhoods in grids: negating either coordinate or switching the coordinates preserves these neighborhoods. Therefore we can decide which cells block (i,j) for i >=0, j>=0, and i >= j and use these symmetries to determine the remaining blocking relations.

Illustrating this rule is the following prototype code written in `R`. This code returns a data structure that will be convenient for determining the "surroundedness" of arbitrary cells in a grid.

``````screen <- function(k=1) {
#
# Returns a data structure:
#   \$offset is an array of offsets
#   \$screened is a parallel array of screened offset indexes.
#   \$distance is a parallel array of distances.
# The first index always corresponds to (0,0).
#
screened.by <- function(xy) {
uv <- abs(xy)
if (reversed <- uv > uv) {
uv <- rev(uv)
}
i <- which.min(c(uv, abs(uv-uv), uv))
ij <- uv + c(floor((1-i)/3), floor(i/3)-1)
if (reversed) ij <- rev(ij)
return(ij * sign(xy))
}
#
# For each lattice point within the circular neighborhood,
# find the unique lattice point that screens it from the origin.
#
xy <- subset(expand.grid(x=(-k:k), y=(-k:k)),
subset=(x^2+y^2 <= k^2) & (x != 0 | y != 0))
g <- t(apply(xy, 1, function(z) c(screened.by(z), z)))
#
# Sort by distance from the origin.
#
colnames(g) <- c("x", "y", "x.to", "y.to")
ij <- unique(rbind(g[, 1:2], g[, 3:4]))
i <- order(abs(ij[,1]), abs(ij[,2])); ij <- ij[i, , drop=FALSE]
rownames(ij) <- 1:length(i)
#
# Invert the "screened by" relation to produce the "screened" relation.
#
# (Row, column) offsets.
ij.df <- data.frame(ij, i=1:length(i))
#
# Distances from the origin (in cells).
distance <- apply(ij, 1, function(u) sqrt(sum(u*u)))
#
# "Screens" relation (represented by indexes into ij).
g <- merge(merge(g, ij.df), ij.df,
by.x=c("x.to", "y.to"), by.y=c("x","y"))
g <- subset(g, select=c(i.x, i.y))
h <- by(g\$i.y, g\$i.x, identity)

return( list(offset=ij, screened=h, distance=distance) )
}
``````

The value of `screen(12)` was used to produce this depiction of this screening relation: arrows point from cells to those that immediately screen them. The hues are proportioned by distance to the origin, which is at the middle of this neighborhood: This computation is fast and needs to be done only once for a given neighborhood. For instance, when looking out 200 m on a grid with 5 m cells, the neighborhood size will be 200/5 = 40 units.

### Step 2: Applying the computation to selected points

The rest is straightforward: to determine whether a cell located at (x,y) (in row and column coordinates) is "surrounded" with respect to this neighborhood data structure, perform the test recursively beginning with an offset of (i,j) = (0,0) (the neighborhood origin). If the value in the polygon grid at (x,y) + (i,j) is nonzero, then visibility is blocked there. Otherwise, we will have to consider all the offsets that could have been blocked at offset (i,j) (which are found in O(1) time using the data structure returned by `screen`). If there are none that are blocked, we have reached the perimeter and conclude that (x,y) is not surrounded, so we stop the computation (and do not bother to inspect any remaining points in the neighborhood).

We can collect even more useful information by keeping track of the furthest line-of-sight distance reached during the algorithm. If this is less than the desired radius, the cell is surrounded; otherwise it is not.

Here is an `R` prototype of this algorithm. It's longer than it seems, because `R` does not natively support the (simple) stack structure needed to implement the recursion, so a stack has to coded, too. The actual algorithm starts about two-thirds of the way through and needs only a dozen lines or so. (And half of those merely handle the situation around the edge of the grid, checking for out-of-range indexes within the neighborhood. This could be made more efficient simply by expanding the polygon grid by `k` rows and columns around its perimeter, eliminating any need for index range checking at the cost of a little more RAM to hold the polygon grid.)

``````#
# Test a grid point `ij` for a line-of-sight connection to the perimeter
# of a circular neighborhood.
#   `xy` is the grid.
#   `counting` determines whether to return max distance or count of stack ops.
#   `perimeter` is the assumed values beyond the extent of `xy`.
#
# Grid values of zero admit light; all others block visibility
# Returns maximum line-of-sight distance found within `nbr`.
#
panvisibility <- function(ij, xy, nbr=screen(), counting=FALSE, perimeter=1) {
#
# Implement a stack for the algorithm.
#
count <- 0 # Stack count
stack <- list(ptr=0, s=rep(NA, dim(nbr\$offset)))
push <- function(x) {
n <- length(x)
count <<- count+n         # For timing
stack\$s[1:n + stack\$ptr] <<- x
stack\$ptr <<- stack\$ptr+n
}
pop <- function() {
count <<- count+1         # For timing
if (stack\$ptr <= 0) return(NULL)
y <- stack\$s[stack\$ptr]
#stack\$s[stack\$ptr] <<- NA # For debugging
stack\$ptr <<- stack\$ptr - 1
return(y)
}
#
# Initialization.
#
m <- dim(xy); n <- dim(xy)
push(1) # Stack the *indexes* of nbr\$offset and nbr\$screened.
dist.max <- -1
#
# The algorithm.
#
while (!is.null(i <- pop())) {
cell <- nbr\$offset[i, ] + ij
if (cell <= 0 || cell > m || cell <= 0 || cell > n) {
value <- perimeter
} else {
value <- xy[cell, cell]
}
if (value==0) {
if (nbr\$distance[i] > dist.max) dist.max <- nbr\$distance[i]
s <- nbr\$screened[[paste(i)]]
if (is.null(s)) {
#exited = TRUE
break
}
push(s)
}
}
if (counting) return ( count )
return(dist.max)
}
`````` In this example, polygonal cells are black. Colors give the maximum line-of-sight distance (out to 50 cells) for non-polygonal cells, ranging from light orange for short distances to dark blue for the longest distances. (The cells are one unit wide and high.) The visibly evident streaks are created by the small polygon "islands" in the middle of the "river": each one blocks a long line of other cells.

### Analysis of the algorithm

The stack structure implements a depth-first search of the neighborhood visibility graph for evidence that a cell is not surrounded. Where cells are far from any polygon, this search will require inspection of only O(k) cells for a radius-k circular neighborhood. The worst cases occur when there are a small number of scattered polygon cells within the neighborhood but even so the boundary of the neighborhood cannot quite be reached: these require inspecting almost all the cells in each neighborhood, which is a O(k^2) operation.

The following behavior is typical of what will be encountered. For small values of k, unless the polygons fill most of the grid, most non-polygonal cells will be obviously unsurrounded and the algorithm scales like O(k). For intermediate values, the scaling starts looking like O(k^2). As k gets really large, most cells will be surrounded and that fact can be determined well before the entire neighborhood is inspected: the algorithm's computational effort thereby reaches a practical limit. This limit is attained when the neighborhood radius approaches the diameter of the largest connected non-polygonal regions in the grid.

As an example, I using the `counting` option coded into the prototype of `screen` to return the number of stack operations used in each call. This measures the computational effort. The following graph plots the mean number of stack ops as a function of neighborhood radius. It exhibits the predicted behavior. We can use this to estimate the computation needed to evaluate 13 million points on a grid. Suppose that a neighborhood of k = 200/5 = 40 is used. Then a few hundred stack operations will be needed on average (depending on the complexity of the polygon grid and where the 13 million points are located relative to the polygons), implying that in an efficient compiled language, at most a few thousand simple numeric operations will be required (add, multiply, read, write, offset, etc.). Most PCs will be able to evaluate the surroundedness of about a million points at that rate. (The `R` implementation is much, much slower than that, because it is poor at this kind of algorithm, which is why it can only be considered a prototype.) Accordingly, we might hope that an efficient implementation in a reasonably efficient and appropriate language--C++ and Python come to mind--could complete the evaluation of 13 million points in a minute or less, assuming the entire polygon grid resides in RAM.

When a grid is too large to fit into RAM, this procedure can be applied to tiled portions of the grid. They only have to overlap by `k` rows and columns; take the maxima at the overlaps when mosaicing the results.

### Other applications

The "fetch" of a body of water is closely related to the "surroundedness" of its points. In fact, if we use a neighborhood radius equal to or greater than the waterbody's diameter, we will create a grid of the (non-directional) fetch at every point in the waterbody. By using a smaller neighborhood radius we will at least obtain a lower bound for the fetch at all the highest-fetch points, which in some applications may be good enough (and can substantially reduce the computational effort). A variant of this algorithm that limits the "screened by" relation to specific directions would be one way to compute fetch efficiently in those directions. Note that such variants require modifying the code for `screen`; the code for `panvisibility` does not change at all.

I can definitely see how one might want to do this with a raster solution, but given even a reduced # of points, I would expect a very large/high-resolution and therefore hard to process grid or set of grids. Given that, I wonder if exploiting topology in a gdb might be more efficient. You could find all the internal voids with something like:

``````arcpy.env.workspace = 'myGDB'
arcpy.CreateTopology_management('myGDB', 'myTopology', '')
arcpy.AddRuleToTopology_management ('myToplogy', 'Must Not Have Gaps (Area)', 'myFeatures', '', '', '')
arcpy.ValidateTopology_management('myTopology', 'Full_Extent')
arcpy.ExportTopologyErrors_management('myTopology', 'myGDB', 'topoErrors')
arcpy.FeatureToPolygon_management('topoErrors_line','topoErrorsVoidPolys', '0.1')`
``````

you can then work with `topoErrorsVoidPolys` in your normal pattern of `Intersect_analysis()` or whatever. You may need to mess around with extracting the polys of interest from `topoErrorsVoidPolys`. @whuber has a number of pretty excellent posts on this kind of stuff elsewhere here on gis.stackexchange.com.

• That's an interesting pre-processing idea. I think it could be readily adapted to incorporate the 200m limit (by buffering and intersection, etc.) Your point about the grids getting pretty large is certainly a real concern. There's no rule in GIS that says a solution to a problem has to be purely raster-based or purely vector-based (although there is a principle that says you should have a pretty good reason to convert from one representation to the other in the middle of an analysis; here, as you suggest, there may be substantial benefit from doing exactly that). – whuber Jan 17 '14 at 16:45

If you really want to go raster...I would do something along the lines of this pseudo code (don't cringe just 'cause it's obvious I'm an AML throwback! :p )

1. rasterize points ("pts_g") and polys ("polys_g"(
2. voids = regiongroup(con(isnull(polys_g), 1))
3. might need to do something to refine voids to eliminate unwanted external polygon/open universe area
4. pts_surrounded = con(voids, pts_g)

Just kind of making that up, so might need refinement.

• Your solution makes no reference to the limiting distance (of, say, 200m), so it does not seem to respond correctly to the question. – whuber Jan 17 '14 at 16:42
• you're right. This applies to my other answer, too. I suppose one could use `Expand()`, but at that point I would think that the answer from @radouxju would be functionally equivalent and probably faster. (nothing against viewshed, just don't use it much). – Roland Jan 17 '14 at 16:44
• was trying to edit when time ran out. want to expand on the `Expand()` to say do that on `pts_g` and just use `Con()` to intersect with `polys_g`. – Roland Jan 17 '14 at 16:52

If you use a threshold distance value (here you talk about 200 m), then the best solution is to use vector analysis:

1) create a 200 m buffer around each point (in black on the illustration)

2) use the intersect tool (analysis) between the buffer and the polygons (in blue on the illustration). It will look nicer if you do this between the boundaries of the surrounding polygons and the buffer, but the final result is the same.

3) use feature to polygon (management) to create polygons where your points are completely surrounded (in red on the illustration)

4) select layers by location (management) or spatial join (analysis) to identify the points that are surrounded. The use of spatial join allows you to have an information about the embedding polygon (area of the polygon, zonal statistics...) which could be useful for further processing.

Alternatives 2b) Depending on your needs, you could select by location the surrounding polygons within a 200 m distance, you can then identify some kinds of "enclosure" but not as strictly as in 2). Considering the "maze case", this could help : evaluate how long it is to "escape" from the location.

• You can already exclude from analysis the points that are fully included or fully free

• then you convert your obstacles to a raster and set the values to NoData where you have a polygon, and to the cell size in meter where you don't (this will make your cost raster).

• third, you can compute the cost distance using the newly generated cost raster

• finally, you use a zonal statistics as table based on the buffer boundaries converted to raster (forming an annulus). If you can escape in all direction, the minimum should be approximatevely 200 (depending on the cell size of your analysis). But if you are in a maze, the maximum will be larger than 200. So the maximum of the zonal statistics minus 200 will be a continuous value indicating how difficult it is to "escape".

• Please clarify your definition of "surrounded." The description in the question suggests that a point should be considered "surrounded" when some part of the polygon is visible in all directions around that point (out to a distance of 200 m). How do you test that in step (3), exactly? (It isn't easy using a vector analysis!) – whuber Jan 17 '14 at 16:44
• I've added a little illustration, it is easier to explain like that. If the buffer does not intersect a polygon in all directions, then the loop won't be closed. And if the loop isn't close, this will not make a polygon. – radouxju Jan 17 '14 at 19:03
• I'm not sure what you mean by "loop" or "closed." Note that a point can be "surrounded" even when no circle of radius r (less than 200 m) around it is entirely contained within the polygon. Think of a maze: the polygon is everything except the corridors in the maze. One can escape from the maze starting at any point within it, but most points will be "surrounded" in the sense that the exterior of the maze will not be visible from them. – whuber Jan 17 '14 at 19:48
• from my understanding, surounded mean somewhere you cannot escape. On the ilustration, you can escape from B but not from A. On the other hand, B would seem to be surrounded if you use viewshed (well, maybe not at 200 m as there is no scale bar on the image, but you would see the polygon boundaries when looking in all directions). I think that we need more details from @Loz – radouxju Jan 17 '14 at 20:13
• This would not be a difficult question at all if "cannot escape" were the criterion to check: just regiongroup the complement of the polygon, keep only the unique external component, and check for inclusion within it. I think a close reading of the question--especially its references to looking at all possible bearings--clarifies the sense in which "surrounded" is intended, although I agree it's stated pretty vaguely. – whuber Jan 17 '14 at 20:16