Running a spatial query on a layer

Topics: Algorithms
May 29, 2013 at 10:39 AM
Edited May 29, 2013 at 10:43 AM
Hi,

I have a layer with polygons, in .shp format. These polygons represent land/islands.

What I want to do, is using these data to be able to determine if between two points there is a polygon - the line between these two points intersects a polygon in the layer.

This is a calculation that I have to make thousands times - I am not interested in presenting these data or my calculations on a map.

Can I achieve something like that with SharpMap? As this is a long running task, performance is a concern.


Thanks,
George J.
Editor
May 29, 2013 at 11:13 AM
Is the point data static or dynamic?

If its dynamic, how often do the points update?

If you have no need to represent the output on a map, then you could probably just use NTS or probably PostGIS.

When you say thousands of points, are you doing the calculation per point (so each point looks at each other point) or per pair (so you have a thousand pairs of points that need the calculation on them?
May 29, 2013 at 11:29 AM
Hi, Robert.

The point data is dynamic - they represent vessels. I use these points in pairs; I do a buffer - unfortunately the version of SharpMap I am using Buffer is not implemented - to find other points inside the buffer, and for each point inside the buffer, I create a line, and do a spatial query if it intersects a polygon in the layer.

I have already tried to use a database; opening and dropping the database connection is killing performance; the query its self has no problems. There are large queues of statements waiting to execute. That's why I am trying to have the layer data in memory, and use a capable spatial query framework, that could speed up the operation.


Thanks,
George J.
Coordinator
May 29, 2013 at 11:49 AM
Edited May 29, 2013 at 12:13 PM
George: Do you want to upgrade SharpMap or not?
Editor
May 29, 2013 at 11:52 AM
Can't you keep the database connection open? or manage that better?

I believe SharpMap is more concerned about displaying maps. As far as I know it uses NTS for the geometry operations (which definitely includes a buffer function) so you might be better just using NTS, rather than going through sharpmap?
Editor
May 29, 2013 at 11:55 AM
Are you referring to me?
May 29, 2013 at 1:17 PM
FObermaier wrote:
George: Do you want to upgrade SharpMap or not?
Actually, I do. One of the reasons I made this post is that the object model I saw in v1.0RC3 is very different from the one I had used a year ago.


George J.
Coordinator
May 29, 2013 at 1:58 PM
In that case, Buffer is now implemented using NTS Buffer operation.

Assuming you have your vessels in a FeatureDataTable that won't change for the iteration sth along the following lines should work.
SharpMap.Data.FeatureDataTable vessels;
var vesselProvider = new SharpMap.Data.Providers.GeometryFeatureProvider(vessels);
SharpMap.Data.Providers.Shapefile landPolygonProvider; //= ...;

foreach(var vesselOfInterest in vessels.Select())
{
    var fds = new SharpMap.Data.FeatureDataSet();
    vesselProvider.ExecuteIntersectionQuery(vesselOfInterest.Geometry.Buffer(bufferValue), fds);
    if (fds.Tables.Count > 0 && fds.Tables[0].Rows.Count > 0)
    {
         foreach(var otherVessels in fds.Tables(0).Select())
         {
             var lineString = vesselOfInterest.Geometry.Factory.CreateLineString(new [] { vesselOfInterest.Geometry.Coordinate, 
                 ((SharpMap.Data.FeatureDataRow)otherVessel).Geometry.Coordinate } );
             var fdsLand = new SharpMap.Data.FeatureDataSet();
             landPolygonProvider.ExecuteIntersectionQuery(lineString, fdsLand))
             if (fdsLand.Tables.Count > 0 && fdsLand.Tables[0].Rows.Count > 0)
             {
                    //Here be dragons
             }
        }
    }
}
This code has not seen a compiler, there may be typos, but the idea should be clear.
May 29, 2013 at 6:52 PM
Hi,

Thank you for the code.

My issue is, I do not have a FeatureDataTable for my points - probably I could create one, but the points data is changing.

This is the code I have come with:
GeoAPI.Geometries.Coordinate p1 = new GeoAPI.Geometries.Coordinate(thisVessel.Lontitude, thisVessel.Lattitude);
GeoAPI.Geometries.Coordinate p2 = new GeoAPI.Geometries.Coordinate(otherVessel.Lontitude, otherVessel.Lattitude);

GeoAPI.GeometryServiceProvider.Instance = new NetTopologySuite.NtsGeometryServices();
GeoAPI.Geometries.IGeometryFactory geometryFactory = GeoAPI.GeometryServiceProvider.Instance.CreateGeometryFactory();

GeoAPI.Geometries.ILineString l = geometryFactory.CreateLineString(new[] { p1, p2 });
SharpMap.Data.FeatureDataSet rs1 = new SharpMap.Data.FeatureDataSet();
GeoDatabase.Instance.LandLayer.ExecuteIntersectionQuery(l, rs1);
int r1 = rs1.Tables.Count > 0 ? rs1.Tables[0].Rows.Count : 0;
but I get the exception:
   at SharpMap.Data.Providers.ShapeFile.GetObjectIDsInView(Envelope bbox) in C:\Projekt\opensource\sharpmap-tfs-v1.0\SharpMap\Data\Providers\ShapeFile.cs:line 686
   at SharpMap.Data.Providers.ShapeFile.ExecuteIntersectionQuery(IGeometry geom, FeatureDataSet ds) in C:\Projekt\opensource\sharpmap-tfs-v1.0\SharpMap\Data\Providers\ShapeFile.cs:line 767
   at Radar.RadarVisibilityExtensions.IsVisible(Vessel OwnShip, Vessel TargetTrack) in c:\Projects\CACS\Tracking\Module.Sensors.Simulator\Radar\RadarVisibilityExtensions.cs:line 92
Is the problem that I use Coordinates, instead of features? Or my GeometryFactory? LandLayer is a SharpMap.Layers.VectorLayer object populated with the following code.
_layer = new SharpMap.Layers.VectorLayer("Aegean");
_layer.DataSource = new SharpMap.Data.Providers.ShapeFile(@"Aegean.shp", true, true, 4326);
_layer.SRID = 4326;
Thanks in advance,
George J.
Coordinator
May 31, 2013 at 8:46 AM
You don't need a FeatureDataTable. The sample assumed you need to get the points within a particular buffer, too. But you already have them and that is fine.

What exactly is the exception you get?
It might be necessary to change
_layer.DataSource = new SharpMap.Data.Providers.ShapeFile(@"Aegean.shp", true, true, 4326);
to
var shp = new SharpMap.Data.Providers.ShapeFile(@"Aegean.shp", true, true, 4326);
shp.DoTrueIntersectionQuery = true;
_layer.DataSource = shp;
Jun 3, 2013 at 7:38 PM
Well, I got it working - not craching - but it doesn't see to return any results; I always get 'false' as results.

This is the full working code:
using SharpMap.Layers;
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Spatial.Los
{
    class Program
    {
        // Aegean.shp can be found in url: http://goo.gl/Yauws

        static void Main(string[] args)
        {
            Vessel[] Vessels = {
                new Vessel() { Height = 300, Lattitude = 37.0055555555556, Lontitude = 25.0513888888889 },
                new Vessel() { Height = 300, Lattitude = 37.3502777777778, Lontitude = 24.8961111111111 },
                new Vessel() { Height = 300, Lattitude = 37.4383333333333, Lontitude = 25.3805555555556 },
                new Vessel() { Height = 5000, Lattitude = 37.4347222222222, Lontitude = 26.2591666666667 },
                new Vessel() { Height = 5000, Lattitude = 37.39, Lontitude = 26.2605555555556 },
                new Vessel() { Height = 5000, Lattitude = 37.6611111111111, Lontitude = 25.3394444444444 },
                new Vessel() { Height = 5000, Lattitude = 37.6502777777778, Lontitude = 25.3888888888889 },
                new Vessel() { Height = 5000, Lattitude = 37.1319444444444, Lontitude = 24.5716666666667 },
                new Vessel() { Height = 0, Lattitude = 37.32, Lontitude = 25.6866666666667 },
                new Vessel() { Height = 0, Lattitude = 37.0258333333333, Lontitude = 24.8586111111111 },
                new Vessel() { Height = 0, Lattitude = 37.3297222222222, Lontitude = 25.1369444444444 },
                new Vessel() { Height = 0, Lattitude = 37.2961111111111, Lontitude = 25.5505555555556 },
                new Vessel() { Height = 0, Lattitude = 37.2405555555556, Lontitude = 25.1611111111111 },
                new Vessel() { Height = 0, Lattitude = 37.1030555555556, Lontitude = 24.8936111111111 },
                new Vessel() { Height = 0, Lattitude = 37.3611111111111, Lontitude = 25.5516666666667 },
                new Vessel() { Height = 0, Lattitude = 37.1066666666667, Lontitude = 24.8283333333333 },
                new Vessel() { Height = 0, Lattitude = 37.2369444444444, Lontitude = 25.0838888888889 },
                new Vessel() { Height = 0, Lattitude = 36.9511111111111, Lontitude = 24.8019444444444 },
                new Vessel() { Height = 0, Lattitude = 37.1719444444444, Lontitude = 25.3758333333333 },
                new Vessel() { Height = 0, Lattitude = 37.1911111111111, Lontitude = 25.5938888888889 },
                new Vessel() { Height = 0, Lattitude = 37.1452777777778, Lontitude = 25.1466666666667 },
                new Vessel() { Height = 0, Lattitude = 37.355, Lontitude = 25.2719444444444 },
                new Vessel() { Height = 0, Lattitude = 37.3538888888889, Lontitude = 25.0105555555556 },
                new Vessel() { Height = 0, Lattitude = 37.0972222222222, Lontitude = 24.7247222222222 },
                new Vessel() { Height = 0, Lattitude = 37.0644444444444, Lontitude = 25.0213888888889 },
                new Vessel() { Height = 0, Lattitude = 37.2647222222222, Lontitude = 24.8188888888889 },
                new Vessel() { Height = 0, Lattitude = 37.4430555555556, Lontitude = 25.7169444444444 },
                new Vessel() { Height = 0, Lattitude = 37.3322222222222, Lontitude = 25.8252777777778 },
                new Vessel() { Height = 0, Lattitude = 37.2466666666667, Lontitude = 25.8375 },
                new Vessel() { Height = 0, Lattitude = 37.1333333333333, Lontitude = 25.8722222222222 },
                new Vessel() { Height = 0, Lattitude = 37.4852777777778, Lontitude = 25.1661111111111 },
                new Vessel() { Height = 0, Lattitude = 37.0886111111111, Lontitude = 25.3058333333333 },
                new Vessel() { Height = 0, Lattitude = 37.3188888888889, Lontitude = 25.3686111111111},
                new Vessel() { Height = 0, Lattitude = 37.2211111111111, Lontitude = 24.6911111111111}
            };

            Vessel OwnShip = new Vessel() { Lattitude = 37.2961111111111, Lontitude = 25.5505555555556 };

            for (int i = 0; i < Vessels.Length; i++)
            {
                Stopwatch sw = Stopwatch.StartNew();

                Vessel OtherShip = Vessels[i];
                bool visible = GeoDatabase.IsVisible(OwnShip, OtherShip);

                sw.Stop();
                Console.WriteLine(string.Format("Vessel {0:00} - Visibility: {1} (Time elapsed: {2}ms)", i + 1, visible, sw.ElapsedMilliseconds));
            }
            Console.WriteLine();
            Console.WriteLine();


            //Stopwatch s1 = Stopwatch.StartNew();
            //Dictionary<int, bool> results = GeoDatabase.CreateVisibilityMatrixByVessel(OwnShip, Vessels);
            //s1.Stop();
            //Console.WriteLine(string.Format("CreateVisibilityMatrixByVessel() - Time elapsed: {0}ms", s1.ElapsedMilliseconds));

            //foreach (var item in results)
            //{
            //    Console.WriteLine(string.Format("Vessel {0} - Visibility: {1}", item.Key, item.Value));
            //}

            Console.Write("Press a key to continue");
            Console.Read();
        }
    }

    public class Vessel
    {
        public double Lontitude { get; set; }
        public double Lattitude { get; set; }
        public double Height { get; set; }
    }

    public class GeoDatabase
    {
        readonly static Lazy<GeoDatabase> _instance = new Lazy<GeoDatabase>(() => { return new GeoDatabase(); });
        private SharpMap.Layers.VectorLayer _layer;

        private GeoDatabase()
        {
            SharpMap.Data.Providers.ShapeFile _shapeFile = new SharpMap.Data.Providers.ShapeFile(@"App_Data\Aegean.shp", true, true, 4326);
            _shapeFile.DoTrueIntersectionQuery = true;
            _layer = new SharpMap.Layers.VectorLayer("Aegean");
            _layer.DataSource = _shapeFile;
            _layer.SRID = 4326;
        }

        public static GeoDatabase Instance
        {
            get { return _instance.Value; }
        }

        public VectorLayer LandLayer
        {
            get { return _layer; }
        }

        public static bool IsVisible(Vessel OwnShip, Vessel OtherShip)
        {
            GeoAPI.Geometries.Coordinate p1 = new GeoAPI.Geometries.Coordinate(OwnShip.Lontitude, OwnShip.Lattitude);
            GeoAPI.Geometries.Coordinate p2 = new GeoAPI.Geometries.Coordinate(OtherShip.Lontitude, OtherShip.Lattitude);

            GeoAPI.GeometryServiceProvider.Instance = new NetTopologySuite.NtsGeometryServices();
            GeoAPI.Geometries.IGeometryFactory geometryFactory = GeoAPI.GeometryServiceProvider.Instance.CreateGeometryFactory();

            GeoAPI.Geometries.ILineString l = geometryFactory.CreateLineString(new[] { p1, p2 });
            SharpMap.Data.FeatureDataSet rs1 = new SharpMap.Data.FeatureDataSet();
            GeoDatabase.Instance.LandLayer.ExecuteIntersectionQuery(l, rs1);
            int r1 = rs1.Tables.Count > 0 ? rs1.Tables[0].Rows.Count : 0;

            GeoAPI.Geometries.IGeometry b = geometryFactory.CreatePoint(p2).Buffer(300);
            SharpMap.Data.FeatureDataSet rs2 = new SharpMap.Data.FeatureDataSet();
            GeoDatabase.Instance.LandLayer.ExecuteIntersectionQuery(b, rs2);
            int r2 = rs2.Tables.Count > 0 ? rs2.Tables[0].Rows.Count : 0;

            return r1 + r2 == 0 ? true : false;
        }

        public static Dictionary<int, bool> CreateVisibilityMatrixByVessel(Vessel OwnShip, Vessel[] Vessels)
        {
            Dictionary<int, bool> results = new Dictionary<int, bool>();

            for (int i = 0; i < Vessels.Length; i++)
            {
                Vessel OtherShip = Vessels[i];
                bool visible = GeoDatabase.IsVisible(OwnShip, OtherShip);
                results.Add(i + 1, visible);
            }

            return results;
        }
    }
}
The data I am using is availiable here: http://goo.gl/Yauws

Any suggestions, pointers?


George J.
Coordinator
Jun 3, 2013 at 9:33 PM
Are you sure that you should get true?

One problem is the size of the buffer you chose. If your layer is lon/lat and you choose a buffer of 300, your buffer covers all of the earth.
I set the buffer value to 0.01 and added a new vessel just 0.0001 degrees off in each direction and got true.

Hth FObermaier
Jun 3, 2013 at 10:56 PM
My 300 value, supposed to be on meters - I have tested the same thing inside SQL Server, running the same spatial queries and using the same data.

What unit does the Buffer function uses?

Ok, the intersections with the testing part with Buffer is wrong - what about the intersections with the line? They seem always to return 0.


Regards,
George J.
Coordinator
Jun 4, 2013 at 7:57 AM
Edited Jun 4, 2013 at 8:27 AM
SharpMap/NTS Buffer operation is unit agnostic. The value should be in the same unit as your data, in case of lon/lat, degrees.

After a look on the provided Shapefile it seems to me that you exchanged longitude and latitude. That is why you don't get any matches.

Edit: That was it, if you exchange longitude and latitude you get both true and false as return value.
Although I must admit I don't understand the logic :)

Hth FObermaier
Jun 4, 2013 at 6:28 PM
Edited Jun 4, 2013 at 8:09 PM
You are right! I do not understand how that happened - all I did was to export data using ogr2ogr utility.

What I am trying to do:
  • Vessels move along islands. As they move they are not allowed to move near the coastline, for security reasons (mostly swimmers)
  • Some vessels are equipped with radars, some not
  • Vessels equipped with radars, they cannot detect other vessels that move near the coastline due to radar limitations - radars can identify tracks that are at least 300m from each other.
Greece generally has a large volume of marine traffic. Due to not conforming to marine rules, yachts mostly, there are marine accidents.

There are a lot of AIS data, that can point out the movement of vessels. Due the complexity of the island coastline and the fact that islands are really close to each other, there is always the excuse "I didn't pick it up in radar".

So, in order to have a "tool" to analyze data, to positively tell, if the radar-lines where clear between the vessels, and there was the possibility for the radar to track a vessel, we have to be able to tell if the line between vessels is clear, and at the same time, the second vessel was clear from the coast by 300m.

I guess now, it's clear the purpose of the above code. There are tons of data produced by the hour, and it's near to impossible to visualize them through a desktop GIS application. These data is stored in a database. It very easy to export a text file and pass it through a command line utility and have some evaluation.

Although, I am afraid the Buffer() function in NTS does not work just like SQL STBuffer function. I want to create a circle with a radius of 300m from a point. The result I get is a polygon, that when I visualize it, it's huge.

How can I implement, a geometry (circle) that it has a 300m radius from a given point as a polygon in NTS?


Grateful for your time,
George J.
Coordinator
Jun 4, 2013 at 8:03 PM
I always have to look up "Great circle math" myself: Here are my favourite sites:
Jun 4, 2013 at 8:58 PM
The second one is like a bible to me. The first one I have to read it.

So, if I understand what you mean; I have to calculate the circle, but I should not use the common circle formula in Euclid's geometry. Instead I should use the distance/bearing formula.

Correct?


George J.
Coordinator
Jun 4, 2013 at 9:20 PM
Yes
Jun 4, 2013 at 10:26 PM
So, the IsVisible() fuction becomes:
        public static bool IsVisible(Vessel OwnShip, Vessel OtherShip)
        {
            GeoAPI.Geometries.Coordinate p1 = new GeoAPI.Geometries.Coordinate(OwnShip.Lontitude, OwnShip.Lattitude, 0);
            GeoAPI.Geometries.Coordinate p2 = new GeoAPI.Geometries.Coordinate(OtherShip.Lontitude, OtherShip.Lattitude, 0);

            GeoAPI.GeometryServiceProvider.Instance = new NetTopologySuite.NtsGeometryServices();
            GeoAPI.Geometries.IGeometryFactory geometryFactory = GeoAPI.GeometryServiceProvider.Instance.CreateGeometryFactory();

            GeoAPI.Geometries.ILineString l = geometryFactory.CreateLineString(new[] { p1, p2 });
            SharpMap.Data.FeatureDataSet rs1 = new SharpMap.Data.FeatureDataSet();
            GeoDatabase.Instance.LandLayer.ExecuteIntersectionQuery(l, rs1);
            int r1 = rs1.Tables.Count > 0 ? rs1.Tables[0].Rows.Count : 0;

            GeoAPI.Geometries.IGeometry b = geometryFactory.CreateEllipse(p2, 0.3, 8);
            SharpMap.Data.FeatureDataSet rs2 = new SharpMap.Data.FeatureDataSet();
            GeoDatabase.Instance.LandLayer.ExecuteIntersectionQuery(b, rs2);
            int r2 = rs2.Tables.Count > 0 ? rs2.Tables[0].Rows.Count : 0;

            Console.WriteLine(string.Format("Feature Count: {0}-{1}", r1, r2));
            return r1 + r2 == 0 ? true : false;
        }
and some extension methods to handle the spherical geometry:
    public static class GeographyExtensions
    {
        public const double EarthRadius = 6371D; // Km

        /// <summary>
        /// Converts degrees to radians
        /// </summary>
        /// <param name="degrees">number in degress</param>
        /// <returns></returns>
        public static double ToRadians(this double degrees)
        {
            return degrees * Math.PI / 180.0D;
        }

        /// <summary>
        /// Converts radians to degrees
        /// </summary>
        /// <param name="radians">number in radians</param>
        /// <returns></returns>
        public static double ToDegrees(this double radians)
        {
            return radians * 180.0D / Math.PI;
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="p1">Start Point</param>
        /// <param name="d">Distance in Kilometers</param>
        /// <param name="bearing">Bearing in degrees</param>
        /// <returns></returns>
        public static GeoAPI.Geometries.Coordinate Destination(this GeoAPI.Geometries.Coordinate p1, double d, double bearing)
        {
            double lon1r = p1.X.ToRadians();
            double lat1r = p1.Y.ToRadians();
            double b = bearing.ToRadians();

            double lat2 = Math.Asin(Math.Sin(lat1r) * Math.Cos(d / EarthRadius) + Math.Cos(lat1r) * Math.Sin(d / EarthRadius) * Math.Cos(b));
            double lon2 = lon1r + Math.Atan2(Math.Sin(b) * Math.Sin(d / EarthRadius) * Math.Cos(lat1r), Math.Cos(d / EarthRadius) - Math.Sin(lat1r) * Math.Sin(lat2));

            return new GeoAPI.Geometries.Coordinate(lon2.ToDegrees(), lat2.ToDegrees());
        }

        public static GeoAPI.Geometries.ILineString CreateEllipse(this GeoAPI.Geometries.IGeometryFactory factory, GeoAPI.Geometries.Coordinate center, double radius, int segmentsPerQuadrant)
        {
            double step = 360 / (segmentsPerQuadrant * 4);
            List<GeoAPI.Geometries.Coordinate> pts = new List<GeoAPI.Geometries.Coordinate>();

            double angle = 0d;
            for (int i = 0; i < 4 * segmentsPerQuadrant; i++)
            {
                GeoAPI.Geometries.Coordinate pt = center.Destination(radius, angle);
                pts.Add(pt);
                angle += step;
            }
            pts.Add(pts[0]);

            return factory.CreateLineString(pts.ToArray());
        }
    }
Thank you, again
Marked as answer by gcapnias on 10/13/2013 at 9:39 AM
Coordinator
Jun 5, 2013 at 7:31 AM
Thanks, actually this is sth that has been requested over and over...
FObermaier