This project is read-only.

Find all Polygons a Linestring intersects

Dec 26, 2013 at 2:40 AM
Good Morning (well, it is a the moment)
My question related to performance/best approach: I have a grid of 10,000 polygons each polygon is 50x50 meters. I have a set of linestrings (around 220 ATM) the cross over the polygons.

The result I'm looking for is: for each polygon I want to know the number of Linestrings that intersect it. At the moment I'm using two loops to check which grids each linestring intersects. Then increasing a counter for that polygon, plus collecting some identification data so I can find the linestring again if needed. But as you can imagine this is taking a very very long time 10,000 * 220 = 2200000 iterations and that's only my start, I actually want to have a large grid of over 1,000,000 squares. And I want to load different day sets of data (so the linestrings for each day will differ)

I was hoping there is a better way of achieving my outcome

Maybe there is a way to estimate (using the direction of the points in the linesting) to identify only the "most likely" polygons to be penetrated thus greatly reducing the number of polygons to be tested?

Any ideas pointers appreciated.

Regards Bruce Quinton
Dec 27, 2013 at 9:16 AM
I guess this problem is of no easy solution.

In the past I had a problem where I needed to understand for 20.000 point shapes which polygon they lied in (among 500).
I had to cycle through all the points and for each point cycle through all the polygons and apply a "point in polygon" function...which is quite slow... so very time consuming.

I achieved a little better performance by checking the single polygon extent (east, west, north and south) to see if the point under analysis could be in...
if the point coordinate are lower or higher than one of the four extent coordinate of the polygon I excluded such polygon from the comparison.

Your problem is that you have line strings, which is not a point or a straight line in principle...so this would be very time consuming.

Since the intersection function is quite computationally demanding, and since vb.net or c#.net are not very computationally efficient languages, I am wandering if it would be better for you to use some scientific language such as c++ or fortran to achieve better performance in this problem. Personally I would prefer fortran also because I know it better but I think c++ should be fine and also with c++ you can build a dll that can be called directly from your .net code as a library.
Since you need to pass all the required data to the routine and read the results back in your .net based application, I guess c++ would be better because you can send your data to the routine as arguments, while if you use fortran you need to work with an external application.

This is how I would proceed...but maybe someone that knows dotspatial better than me would be able to advice you better.
However if you find a dotspatial based solution let us know

Oscar
Dec 27, 2013 at 10:03 PM
Edited Dec 28, 2013 at 2:04 PM
You should go with NetTopologySuite's PreparedGeometry functionality. There is a conversion project between DotSpatial geometries/shapes and GeoAPI/NTS geometries if you rely on other DotSpatial functionality.

the code would then look somewhat like this:
//get line strings from somewhere
var dsLineStrings = new List<DotSpatial.Topology.ILineString>();
// get polygons from somewhere else
var dsPolygons = new Dictionary<DotSpatial.Topology.IPolygon, int>();
var polygons = new Dictionary<GeoAPI.Geometries.IPolygon, int>();

// build a quadtree index on the polygons
var quadTree = new NetTopologySuite.Index.Quadtree.Quadtree<IPolygon>();
var i = 0;
foreach (var dsPolygon in dsPolygons.Keys)
{
        var polygon = dsPolygon.ToGeoAPI(); 
        polygons.Add(polygon, 0);
    quadTree.Insert(polygon.EnvelopeInternal, polygon);
}

// iterate over all line strings
foreach (var dsLineString in dsLineStrings)
{
    // convert to GeoAPI linestring
        var lineString = dsLineString.ToGeoAPI();
        // build a prepared geometry, which is optimized for a one against many comparison
    var prepLineString = NetTopologySuite.Geometries.Prepared.PreparedGeometryFactory.Prepare(lineString);
        // search quad tree for polygon candidates
    foreach (var candidate in quadTree.Query(lineString.EnvelopeInternal))
    {
                // Do intersection test and increment counter
        if (prepLineString.Intersects(candidate))
        {
                        // Does this really work, I have not tried.
            polygons[candidate]++;
        }
    }
}
NOTE: Code has not seen a compiler, there may be typos and other flaws but you should get the picture.
Dec 30, 2013 at 1:55 AM
Thank you both for your suggestions, FObermaier I will give the Code a go shortly and test the results against the code I ended with myself.

My approach was,
Using the outer edge of the Gird, I intersection each Linestring, finding the starting "grid square' I intersection'd that each pair of points along the line string, stepping to the next grid the line touched (using the direction of the line string points and finding the edge of the last intersecting point).. I have written this description poorly I think it has made the procedure a lot more complicated than needed.

My code is able to process 6000 line strings (average length of 40 nm) over 2.5 million squaregrid (50mx50x) within 5 minutes.

Thanks again, and I will let you know how the code above goes.
Jan 2, 2014 at 9:00 AM
Hi Bruce

Maybe I am missing something. Are you saying you have a regular grid of square polygons (with sides oriented north-south and east-west) and you want to see which of these grid polygons (which are in principle the same as raster "cells") are crossed by your linestrings?

If this is the case maybe there is a rather faster option

furthermore, what is the average length of your line strings? 40 nm? nanometers?

Oscar