Either rasterizing a shapefile or fast determination of which polygon contains a point?

Developer
Apr 26, 2012 at 4:48 PM

I am converting code from MapWindow 4 that overlays a set of layers (grids and shapefiles).

In MW4, there is a function called MapWinGIS.Shapefile.PointInShapefile that can fairly quickly return the index of the shape in a polygon shapefile that contains a given point. There does not seem to be an equivalent function in DotSpatial, so we have written a function that loops through each Feature in a FeatureSet and uses the .BasicGeometry function on each Feature, then uses .Intersects on the BasicGeometry to figure out which feature intersects the point. This probably works correctly, but it is taking hours to run so I am looking for a faster way.

A fast version of this function would be a straightforward way to go. Another option would be to make a raster version of the shapefile that lines up exactly with the grid layer we are comparing to. From there, overlay would be easy and fast.

I have found DotSpatial.Analysis.VectorToRaster, but it uses GDI+ and that imposes limitations and I have also not gotten it to produce a raster containing anything other than all zero.

Has anyone used DotSpatial to do either of these things?

Developer
Apr 26, 2012 at 11:38 PM

Unfortunately anything built using NTS Topology is ridiculously slow.  However, if you are loading the shapefile directly from disk, instead of looking at the FeatureSet at all (basically try never to touch the "Features" property at all.  Instead, there should be something like Shapes or ShapeList or some such thing.  You can cycle through through that and each item should be a "ShapeRange" object.  The ShapeRange supports and "Intersects" function.  So if the shapes are polygons it directs it to the code in DotSpatial.Data.PolygonShape which actually does the calculation.  It should be loads faster than using features.

By default, just remember anything NTS is super slow.  The only thing I would ever really use it for is overlay calculations where you are creating new polygons from intersecting old polygons or something.  Simple intersection testing doesn't need something that heavy and slow.  You can do almost everything from using ShapeRanges.

Ted

Editor
Apr 27, 2012 at 6:19 AM

I assume that the use of an spatial index ([Quad|R]Tree) on the FeatureSet has more impact on the calculation than using either Intersects algorithm.

Regarding the complaint about NTS being super slow. The version used in DotSpatial is 1.7+, NTS is now at version 1.12, and I know there are not only robustness improvements. You may want to checkout http://nuget.org/packages/NetTopologySuite.DotSpatial.Converter.

Hth FObermaier

Developer
Apr 27, 2012 at 4:57 PM

I was never sure if it was an NTS thing or just something that went wrong with the way I implemented it in DotSpatial.Topology.  The "Feature" object was always my creation, and it always ended up disappointing me in both memory and speed tests.  However when Joe tried to do a side by side trial with the real NTS, he had trouble setting up the objects and so I never got around to getting a clean comparison on performance.  No doubt the modern version has been improved.  I always envisioned that the only way it could be as slow as it is is if it was doing something unintended.  Alternately, maybe it could be something that was perfectly fast in Java's JTS that didn't translate verbatim and so became slower in C#?  If either is the case, then maybe it has been fixed by now.  It would be worth investigating.  Definitely NTS is much more fully featured and robust, so if it is faster now, then that is the preferred tool to use.

Ted   

Developer
Apr 30, 2012 at 2:57 PM

I'm using the IFeatureSource file-based implementation that Ted and I worked on a while back to search for shapes.  I am happy with the performance results. If you are not interested in attributes, you can just use the IShapeSource implementation which is a subset of the IFeatureSource approach.  If you are doing repetitive searches, there is a Quad Tree option to speed up your searches.  If you are interested in this, I will try to give some more info on the usage.  This approach also has the advantage of low memory usage.

Kyle

Jul 20, 2012 at 9:02 AM

Hello,

I'm trying to bulid an application that (amongst other things) opens a polygon shapefile, converts it to a raster (grid, even in Ram or as an array of integer values might do) and passes it to a .dll which processes it and so on.

I first inteded to convert the shp2grid plugin for MapWindow 4 into DotSpatial. but now I found this thread which deals with a similar problem.

vatavian, did aou already find a satisfactory way to go for it?

I'd much appreciate your feedback before waisting too much time in a similar problem.

 

Thanks, Claudia

Developer
Jul 20, 2012 at 5:47 PM

My current code is in my class that is a wrapper around a DotSpatial featureset.

This function returns the index of the shape that contains the given point or -1 if no shape does.

Keeping track of the most recent shape index is crucial to performance because we can test the latest one first and that is most often the right one as we scan through the grid cells.

Private pLastCoordinateIndex As Integer = -1

Function CoordinatesInShapefile(ByVal aX As Double, ByVal aY As Double) As Integer

        With Me.AsFeatureSet()
            If .Extent.Intersects(aX, aY) Then
                Dim lCoordinate As New DotSpatial.Topology.Coordinate(aX, aY)
                Dim lShapeIndices = .ShapeIndices

                If pLastCoordinateIndex > -1 AndAlso lShapeIndices(pLastCoordinateIndex).Intersects(lCoordinate) Then
                    Return pLastCoordinateIndex
                End If

                Dim lIndex As Integer = 0
                For Each lShapeRange As DotSpatial.Data.ShapeRange In .ShapeIndices
                    If lIndex <> pLastCoordinateIndex AndAlso lShapeRange.Intersects(lCoordinate) Then
                        pLastCoordinateIndex = lIndex
                        Return lIndex
                    End If
                    lIndex += 1
                Next
            End If
        End With

       Return -1

End Function