# How to determine the location that has the most intersections with other polygons in a certain distance/buffer?

I have a very scattered map with a lot of very differently shaped polygons. Now I want to determine the location where a circle polygon of a certain size intersects the most surface of these scattered polygons.

It would be even better to be able to determine not only one best location but maybe even some more locations that are the most intersecting to get something like a heatmap.

Here is an exemplary picture:

-

Have you considered using a raster approach?

Firstly Convert shaped polygons into a raster data set with cell values of 1 where a cell is covered by a polygon.

The focal statistics tool calculates statistics of surrounding cells for each cell in the raster.

• Set Neighborhood type to circle
• Set the radius of the circle (in map units) based on the radius of your circle
• Set the statistic type to SUM

The resulting raster will show a heat map for the amount of intersected polygons for a circle with centroid in each raster cell.

Any locations where the maximum is attained (or almost attained) are the centers of interest:

The focal sum values range from 0 to 16697 cells. The little bit of the grid appearing in the lower right shows all the places where the sums are 16000 and greater. The circle is drawn with a center at the one cell attaining the maximum of 16697. Circles centered at the other high-focal-sum points will be almost as good.

This method involves a certain amount of approximation depending on how you rasterize the initial polygon layer.

-
+1 This is the ideal approach. – whuber Sep 12 '13 at 17:08
Great explanation! Thank you very much! – Gideon Sep 12 '13 at 17:57
Ok, I tried it, but I think in the Point statistics tool you can only work with point or multipoint vectordata. I am not able to choose my raster data there. – Gideon Sep 13 '13 at 12:23
I guess you were refering to the Focal Statistics tool instead of the Point Statistics. I tried to use it but unfortunately, the tool seems to be buggy as already mentioned here: support.esri.com/en/bugs/nimbus/TklNMDU5NzI2 So this was my result: :( abload.de/img/issuesqgugt.jpg – Gideon Sep 13 '13 at 13:40
The examples were produced with ArcGIS 10.1 on a grid in which the circle radius was 125 cells. If there is indeed a bug, then a workaround is to downsample the raster (so that the focal neighborhood uses fewer cells). As the example shows, high resolution usually is not needed to identify an optimal position or set of positions. But since Spatial Analyst has proved itself unreliable with this calculation, I would recommend using different software. Any software that can perform convolutions, Fourier transforms, or flexible image blurring will work. – whuber Sep 13 '13 at 14:37

You can do this in Whitebox Geospatial Analysis Tools without the need to convert the vector polygons file to a raster using the following Groovy script. Download the latest version of Whitebox (https://code.google.com/p/whitebox-geospatial-analysis-tools/downloads/list), open the Scripter, paste the script and save it as IntersectionDensity.groovy. You'll find that the script appears as a plugin tool in the Vector Tools toolbox. It takes a shapefile of POLYLINE or POLYGON base ShapeType as the input, along with the buffer size and the output raster grid resolution. It wouldn't be difficult to modify it to also output a vector point coincident with the highest intersection point if that is what you would like (contact me if you need help with it--jlindsay(at)uoguelph.ca). I'll add the tool to the next official release of Whitebox.

``````/*
* Copyright (C) 2013 Dr. John Lindsay <jlindsay@uoguelph.ca>
*
* This program is free software: you can redistribute it and/or modify
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

import java.awt.event.ActionListener
import java.awt.event.ActionEvent
import java.util.concurrent.Future
import java.util.concurrent.*
import java.util.Date
import java.util.ArrayList
import whitebox.interfaces.WhiteboxPluginHost
import whitebox.geospatialfiles.ShapeFile
import whitebox.geospatialfiles.shapefile.*
import whitebox.geospatialfiles.WhiteboxRaster
import whitebox.geospatialfiles.WhiteboxRasterBase
import whitebox.geospatialfiles.WhiteboxRasterBase.DataType
import whitebox.ui.plugin_dialog.ScriptDialog
import whitebox.utilities.StringUtilities
import whitebox.structures.KdTree
import groovy.transform.CompileStatic

// The following four variables are required for this
// script to be integrated into the tool tree panel.
// Comment them out if you want to remove the script.
def name = "IntersectionDensity"
def descriptiveName = "Intersection Density"
def description = "Calculates the spatial pattern of polygon intersection with a buffer distance"
def toolboxes = ["VectorTools"]

public class IntersectionDensity implements ActionListener {
private WhiteboxPluginHost pluginHost
private ScriptDialog sd;
private String descriptiveName

public IntersectionDensity(WhiteboxPluginHost pluginHost,
String[] args, def descriptiveName) {
this.pluginHost = pluginHost
this.descriptiveName = descriptiveName

if (args.length > 0) {
final Runnable r = new Runnable() {
@Override
public void run() {
execute(args)
}
}
t.start()
} else {
// Create a dialog for this tool to collect user-specified
// tool parameters.
sd = new ScriptDialog(pluginHost, descriptiveName, this)

// Specifying the help file will display the html help
// file in the help pane. This file should be be located
// in the help directory and have the same name as the
// class, with an html extension.
def helpFile = "IntersectionDensity"
sd.setHelpFile(helpFile)

// Specifying the source file allows the 'view code'
// button on the tool dialog to be displayed.
def pathSep = File.separator
def scriptFile = pluginHost.getResourcesDirectory() + "plugins" + pathSep + "Scripts" + pathSep + "FilterUserDefinedWeights.groovy"
sd.setSourceFile(scriptFile)

// add some components to the dialog
sd.addDialogFile("Input file", "Input Vector Polygon File:", "open", "Vector Files (*.shp), SHP", true, false)
sd.addDialogFile("Output file", "Output Raster File:", "close", "Raster Files (*.dep), DEP", true, false)
sd.addDialogDataInput("Buffer Distance:", "Enter a buffer distance", "1000", true, false)
sd.addDialogDataInput("Output Grid Resolution:", "Enter a grid resolution", "50", true, false)

// resize the dialog to the standard size and display it
sd.setSize(800, 400)
sd.visible = true
}
}

// The CompileStatic annotation can be used to significantly
// improve the performance of a Groovy script to nearly
// that of native Java code.
@CompileStatic
private void execute(String[] args) {
try {
int i, f, progress, oldProgress
double x1, x2, y1, y2, dist1, dist2
if (args.length != 4) {
pluginHost.showFeedback("Incorrect number of arguments given to tool.")
return
}
String inputFile = args[0]
String outputFile = args[1]
double cellSize = Double.parseDouble(args[3])

def input = new ShapeFile(inputFile)

// make sure that input is of a POLYLINE or POLYGON base shapetype
if (input.getShapeType().getBaseType() != ShapeType.POLYLINE
&& input.getShapeType().getBaseType() != ShapeType.POLYGON) {
pluginHost.showFeedback("Input shapefile must be of a POLYLINE or POLYGON base shapetype.")
return
}

double[][] points

double x
double y
double z = 0.0
List<KdTree.Entry<Double>> results
int recNum = 0
int numFeatures = input.getNumberOfRecords()

// find out how many points there are to add to the kd-tree
def numPoints = 0
progress = 0
oldProgress = -1
for (ShapeFileRecord record : input.records) {
recNum++
f = record.getRecordNumber() - 1

points = record.getGeometry().getPoints()

numPoints += record.getGeometry().getPoints().length

// update progress
progress = (int)(100.0 * recNum / numFeatures)
if (progress != oldProgress) {
oldProgress = progress
pluginHost.updateProgress("Loop 1 of 3:", progress)
}
// check to see if the user has requested a cancellation
if (pluginHost.isRequestForOperationCancelSet()) {
pluginHost.showFeedback("Operation cancelled")
return
}
}

KdTree<Double> pointsTree = new KdTree.SqrEuclid<Double>(2, new Integer(numPoints))

progress = 0
oldProgress = -1
for (ShapeFileRecord record : input.records) {
recNum++
f = record.getRecordNumber() - 1

points = record.getGeometry().getPoints()

for (i = 0; i < points.length; i++) {
def entry = new double[2]
entry[1] = points[i][0]
entry[0] = points[i][1]

}

// update progress
progress = (int)(100.0 * recNum / numFeatures)
if (progress != oldProgress) {
oldProgress = progress
pluginHost.updateProgress("Loop 2 of 3:", progress)
}
// check to see if the user has requested a cancellation
if (pluginHost.isRequestForOperationCancelSet()) {
pluginHost.showFeedback("Operation cancelled")
return
}
}

// create the output raster
def rows = (int) (Math.ceil((input.getyMax() - input.getyMin()) / cellSize))
def cols = (int) (Math.ceil((input.getxMax() - input.getxMin()) / cellSize))
def north = input.getyMax()
def south = input.getyMax() - rows * cellSize
def east = input.getxMin() + cols * cellSize
def west = input.getxMin()
double noData = -32768.0
def dataType = WhiteboxRasterBase.DataType.FLOAT

def output = new WhiteboxRaster(outputFile, north, south, east, west,
rows, cols, WhiteboxRasterBase.DataScale.CONTINUOUS,
dataType, 0.0, noData);
output.setPreferredPalette("spectrum.plt")

// interpolate into the output raster
def numNeighbouringFeatures = 0
oldProgress = -1
for (int row = 0; row < rows; row++) {
for (int col = 0; col < cols; col++) {
x = output.getXCoordinateFromColumn(col)
y = output.getYCoordinateFromRow(row)
double[] entry = new double[2]
entry[0] = y
entry[1] = x
if (results.size() > 0) {
def neighbouringFeatures = new boolean[numFeatures]
for (KdTree.Entry entry2 : results) {
i = (int)entry2.value
neighbouringFeatures[i] = true
}

int intersections = 0
for (f = 0; f < numFeatures; f++) {
if (neighbouringFeatures[f]) {
ShapeFileRecord record = input.getRecord(f)
points = record.getGeometry().getPoints()

// scan each line segment for an intersection
for (i = 1; i < points.length; i++) {
x1 = points[i - 1][0]
y1 = points[i - 1][1]
x2 = points[i][0]
y2 = points[i][1]
if (lineSegmentInstersectsCircle(x1,
y1, x2, y2, x, y, neighbourhoodRadius)) {
intersections++
}
}
}
}

output.setValue(row, col, intersections)

}

}
// update progress
progress = (int) (100f * row / (rows - 1))
if (progress != oldProgress) {
oldProgress = progress
pluginHost.updateProgress("Loop 3 of 3:", progress)
}
// check to see if the user has requested a cancellation
if (pluginHost.isRequestForOperationCancelSet()) {
pluginHost.showFeedback("Operation cancelled")
return
}
}

+ descriptiveName + " tool.")
output.close()

// display the output image
pluginHost.returnData(outputFile)

// reset the progress bar
pluginHost.updateProgress(0)
} catch (Exception e) {
pluginHost.showFeedback(e.getMessage())
}
}

public static boolean lineSegmentInstersectsCircle(double x1, double y1,
double x2, double y2, double cX, double cY, double r) {
double a = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y1 - y1);
double b  = 2.0 * ((x2 - x1) * (x1 - cX) +(y2 - y1) * (y1 - cY));
double cc = cX * cX + cY * cY + x1 * x1 + y1 * y1 - 2.0 * (cX * x1 + cY * y1) - r * r;
double deter = b * b - 4.0 * a * cc;
if (deter <= 0 ) {
return false;
} else {
double e = Math.sqrt(deter);
double u1 = ( - b + e ) / (2 * a );
double u2 = ( - b - e ) / (2 * a );
if ((u1 < 0 || u1 > 1) && (u2 < 0 || u2 > 1)) {
return false;
} else {
return true;
}
}
}

@Override
public void actionPerformed(ActionEvent event) {
if (event.getActionCommand().equals("ok")) {
final def args = sd.collectParameters()
sd.dispose()
final Runnable r = new Runnable() {
@Override
public void run() {
execute(args)
}
}
t.start()
}
}
}

if (args == null) {
pluginHost.showFeedback("Plugin arguments not set.")
} else {
def f = new IntersectionDensity(pluginHost, args, descriptiveName)
}
``````

This is an example of the result of the script, applied to a vector stream network...

-
First of all: Thank you for this interesting approach! Nevertheless, the script is not working for me. As soon as I run the script after selecting the input shapefile and the output raster data, on the lower right corner it will say "Loop 2 of 3:" and shortly after that a message window will pop up including "2" and the OK-Button. The lower right expression remains the same and nothing happens. A .dep file will be created but it seems to be unusable. Here is a screenshot: abload.de/img/issues76b1i.jpg – Gideon Sep 13 '13 at 11:06
I think the problem was the way that I copied and pasted the script into the forum text box and formatted it. I've updated it. Give it a try and let me know if you have further problems. I will put the tool in Whitebox for the next release. – user21951 Sep 13 '13 at 15:37
@Gideon, just wanted to let you know that I released version 3.0.5 of Whitebox GAT today, which includes this plugin tool if you'd like to give it a try. – user21951 Sep 20 '13 at 20:15