ClosestPointOnPolyline if SHP is in WGS84 SHP (4326)

[dusko23]

I have just noticed that the function ClosestPointOnPolyline dramatically slows down the calculation if the SHP file is 4326 type. This might be also the case in some other calculations. Probably some functionality that successfully draws different EPSG shape files on SFmap object is missing here. If I convert WGS84 shp to PseudoMercator (3857) ClosestPointOnPolyline works fine. The test was conducted on the files downloaded from here http://download.geofabrik.de/

[wfletcher]

Missed this one sorry, There are multiple ClosestPointOnPolyline methods. Can you elaborate on which method you are calling? Are you referring to GeometryAlgorithms.ClosestPointOnPolyline

[dusko23]

Yes, that is the one. I don't know where the other methods are. Please tell. However, GeometryAlgorithms.ClosestPointOnPolyline causes a heavily delayed response when the WGS84 shape file is tested in the MouseMove event.

[wfletcher]

I can't really see why WGS84 coordinates would cause the method to be slower than using PseudoMercator coordinates , There are 6 overloaded ClosestPointOnPolyline methods so you may need to post sample code showing how you call the method

[dusko23]

OK, I was hoping that you would download one shp file from http://download.geofabrik.de/ and give it a try :-) No problem, I just need to grab some time and create a solution with a big shp file included to demonstrate the mouse move issue. This is my function

[wfletcher]

I have run some tests on some data I downloaded from geofabrik and I am seeing similar performance using WGS84 geodetic coordinates or PseudoMercator (3857) projected coordinates.

[Test]
public void TestClosestPointOnPolyLine()
{
    const int Iterations = 25000;
    string testShapeFileName = System.IO.Path.Combine(Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location), "data", "pmar.shp");

    using (ShapeFile shapeFile = new ShapeFile(testShapeFileName))
    {
        //read the first record geometry

        var geometry = shapeFile.GetShapeDataD(0);

        var points = geometry[0];

        var testPoint = points[6];
        testPoint.X += 0.00008;
                testPoint.Y += 0.00005;

                Console.Out.WriteLine("points.Length:" + points.Length);

        PolylineDistanceInfo pdi = PolylineDistanceInfo.Empty;

        DateTime tick = DateTime.Now;
        for (int n = 0; n < Iterations; ++n)
        {
            GeometryAlgorithms.ClosestPointOnPolyline(points, 0, points.Length, testPoint, out pdi);
        }
               DateTime tock = DateTime.Now;

        Console.Out.WriteLine("WGS84 ClosestPointOnPolyline time:{0}ms",  tock.Subtract(tick).TotalMilliseconds );
                Console.Out.WriteLine("pdi point index:{0}, tValue:{1:0.000}",  pdi.PointIndex, pdi.TVal);

        double distance = EGIS.ShapeFileLib.ConversionFunctions.GeodesicDistanceAndBearingBetweenLatLonPoints(testPoint.Y, testPoint.X, pdi.PolylinePoint.Y, pdi.PolylinePoint.X).Item1;

        Console.Out.WriteLine("distance:" + distance);

        ICRS pseudoMerc = CoordinateReferenceSystemFactory.Default.GetCRSById(CoordinateReferenceSystemFactory.Wgs84PseudoMercatorEpsgCode);

        using (ICoordinateTransformation transformation = CoordinateReferenceSystemFactory.Default.CreateCoordinateTrasformation(shapeFile.CoordinateReferenceSystem, pseudoMerc))
        {
        var transformedPoints = (PointD[])points.Clone();
            transformation.Transform(transformedPoints, TransformDirection.Forward);

            var transformedTestPoint = transformation.Transform(testPoint, TransformDirection.Forward);

            tick = DateTime.Now;
            for (int n = 0; n < Iterations; ++n)
            {
                GeometryAlgorithms.ClosestPointOnPolyline(transformedPoints, 0, transformedPoints.Length, transformedTestPoint, out pdi);
            }
            tock = DateTime.Now;

            Console.Out.WriteLine("Projected Coords ClosestPointOnPolyline time:{0}ms", tock.Subtract(tick).TotalMilliseconds);
            Console.Out.WriteLine("pdi point index:{0}, tValue:{1:0.000}", pdi.PointIndex, pdi.TVal);
        }

    //repeat wgs84 test again
        tick = DateTime.Now;
        for (int n = 0; n < Iterations; ++n)
        {
            GeometryAlgorithms.ClosestPointOnPolyline(points, 0, points.Length, testPoint, out pdi);
        }
        tock = DateTime.Now;

        Console.Out.WriteLine("WGS84 ClosestPointOnPolyline time:{0}ms", tock.Subtract(tick).TotalMilliseconds);
        Console.Out.WriteLine("pdi point index:{0}, tValue:{1:0.000}", pdi.PointIndex, pdi.TVal);

        distance = EGIS.ShapeFileLib.ConversionFunctions.GeodesicDistanceAndBearingBetweenLatLonPoints(testPoint.Y, testPoint.X, pdi.PolylinePoint.Y, pdi.PolylinePoint.X).Item1;
        Console.Out.WriteLine("distance:" + distance);
    }
}

Test results:

points.Length: 117 WGS84 ClosestPointOnPolyline time:18.0154ms pdi point index:5, tValue:0.790 distance:9.81777630148275 Projected Coords ClosestPointOnPolyline time:16.9783ms pdi point index:5, tValue:0.789 WGS84 ClosestPointOnPolyline time:17.989ms pdi point index:5, tValue:0.790 distance:9.81777630148275

17ms vs 18ms for 25000 iterations of checking for closest point on a polyline with 117 points. I don't think that the difference is significant.

[dusko23]

I have checked what is causing the problem when using the WGS84 shape file. I have finally isolated this code:

If a big SHP file contains WGS84 coordinates it takes 2+ seconds to find the relevant record. A correct record is found. Variable myCheckDistance is 5 pixels converted to meters.

Any other EPSG delivers the record instantly. I am seeking a record because I perform some other record points calculations to serve GUI.

[wfletcher]

Try passing you map's CoordinateReferenceSystem to the GetClosestPoint method.

public int GetClosestShape(PointD centre, double radius, ICRS sourceCRS=null)

If you call GetClosestShape with the radius parameter set to 5 you are searching for the closest record within 5 degrees for a WGS84 shapefile, which is a massive search area.

By using the optional 3rd parameter you are specifying that the centre and radius parameter's units are using sourceCRS units.

//example event handler
private void sfMap1_MouseMove(object sender, MouseEventArgs e)
{
    PointD pt = sfMap1.PixelCoordToGisPoint(new Point(e.X, e.Y));   
    double radius = 8 / sfMap1.ZoomLevel; //8 pixels

    for (int n = sfMap1.ShapeFileCount - 1; n >= 0; --n)
    {
        int selectedIndex = sfMap1[n].GetClosestShape(pt, radius, sfMap1.MapCoordinateReferenceSystem);

        //if selectedIndex >=0 do what you need to do
        if (selectedIndex >= 0)
        {
            Console.Out.WriteLine("selectedIndex is :" + selectedIndex);
        }
    }            
}

[dusko23]

Thanks, that worked like a charm. I figured out that the same goes for GetShapeIndiciesIntersectingCircle. Could you please give just a short explanation about which conversions you perform in this case? I think that this is to be used only if the SHP in focus has WGS84 EPSG. If SFmap EPSG is projected (meters) and SHP EPSG is projected (meters) but EPSGs are different, then SFmap cursor coordinate conversion to SHP EPSG is still needed,

[dusko23]

The solution below is working fine for me. Given myCHECKdistance in meters I calculate delta longitude as "degree distance" based on myCHECKdistance and latitude. myCHECKdistanceWGS84 is calculated. This is the formula provided by chatGPT

Usage