Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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:

enter image description here

share|improve this question

2 Answers 2

up vote 3 down vote accepted

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.

Polygon indicator

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.

Focal sum

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

Solution

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.

share|improve this answer
    
+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
 * it under the terms of the GNU General Public License as published by
 * 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)
            }
        }
        final Thread t = new Thread(r)
        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
        }
        // read the input parameters
        String inputFile = args[0]
        String outputFile = args[1]
        double neighbourhoodRadius = Double.parseDouble(args[2])
        double radiusSquared = neighbourhoodRadius * neighbourhoodRadius
        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]

                pointsTree.addPoint(entry, (double)f)
            }

            // 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
                results = pointsTree.neighborsWithinRange(entry, neighbourhoodRadius)
                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
            }
        }


        output.addMetadataEntry("Created by the "
            + descriptiveName + " tool.")
        output.addMetadataEntry("Created on " + new Date())
        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)
            }
        }
        final Thread t = new Thread(r)
        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...

enter image description here

enter image description here

share|improve this answer
1  
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
1  
@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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.