Find Intersection Point Between 2 LineStrings

1.7k views Asked by At

I created a formula to form a grid on the google earth. I want to get the intersection point between the lat/long. Please tell me how we can get the intersection.I am using SharpKML lib for generating KML

for (int x = 90; x >= 0; x = x - 15)
                {
                    Placemark placemark = new Placemark();
                    LineString line = new LineString();
                    CoordinateCollection co = new CoordinateCollection();
                    for (int i = 0; i <= 180; i = i + 15)
                    {
                        Vector cords = new Vector()
                            {
                                Latitude = x,
                                Longitude = i,
                                Altitude = 1000
                            };

                        co.Add(cords);
                    }
                    for (int i = -180; i <= 0; i = i + 15)
                    {
                        Vector cords = new Vector()
                        {
                            Latitude = x,
                            Longitude = i,
                            Altitude = 1000
                        };

                        co.Add(cords);
                    }
                    line.Coordinates = co;
                    placemark.Geometry = line;
                    document.AddFeature(placemark);
                }
            for (int x = -90; x <= 0; x = x + 15)
            {
                Placemark placemark = new Placemark();
                LineString line = new LineString();
                CoordinateCollection co = new CoordinateCollection();
                for (int i = 0; i <= 180; i = i + 15)
                {
                    Vector cords = new Vector()
                    {
                        Latitude = x,
                        Longitude = i,
                        Altitude = 1000
                    };

                    co.Add(cords);
                }
                for (int i = -180; i <= 0; i = i + 15)
                {
                    Vector cords = new Vector()
                    {
                        Latitude = x,
                        Longitude = i,
                        Altitude = 1000
                    };

                    co.Add(cords);
                }
                line.Coordinates = co;
                placemark.Geometry = line;
                document.AddFeature(placemark);
            }

            for (int i = 0; i <= 180; i = i + 15)
            {
                Placemark placemark = new Placemark();
                LineString line = new LineString();
                CoordinateCollection co = new CoordinateCollection();
                for (int x = 0; x <= 90; x = x + 15)
                {
                    Vector cords = new Vector()
                    {
                        Latitude = x,
                        Longitude = i,
                        Altitude = 1000
                    };

                    co.Add(cords);
                }
                for (int x = -90; x <= 0; x = x + 15)
                {
                    Vector cords = new Vector()
                    {
                        Latitude = x,
                        Longitude = i,
                        Altitude = 1000
                    };

                    co.Add(cords);
                }
                line.Coordinates = co;
                placemark.Geometry = line;
                document.AddFeature(placemark);
            }
            for (int i = -180; i <= 0; i = i + 15)
            {
                Placemark placemark = new Placemark();
                LineString line = new LineString();
                CoordinateCollection co = new CoordinateCollection();
                for (int x = 0; x <= 90; x = x + 15)
                {
                    Vector cords = new Vector()
                    {
                        Latitude = x,
                        Longitude = i,
                        Altitude = 1000
                    };

                    co.Add(cords);
                }
                for (int x = -90; x <= 0; x = x + 15)

                {
                    Vector cords = new Vector()
                    {
                        Latitude = x,
                        Longitude = i,
                        Altitude = 1000
                    };

                    co.Add(cords);
                }
                line.Coordinates = co;
                placemark.Geometry = line;
                document.AddFeature(placemark);
            }
2

There are 2 answers

0
Matthew On

I think the DotSpatial library should meet your needs, I have used this library in the past but not made use of the intersections function:

http://dotspatial.codeplex.com/wikipage?title=DotSpatial.Data.FeatureSetExt.Intersection

If you try and do your own line intersection analysis, know that a simplistic Cartesian plane approach will introduce errors (which [I think] become more obvious as you approach the poles).

See here: http://www.geog.ubc.ca/courses/klink/gis.notes/ncgia/u32.html

And here: Intersection between two geographic lines

0
Ted On

Matthew is correct if the question is how you can find the intersection point of an arbitrary LineString object with your grid using C#. In C++ you can use GEOS http://trac.osgeo.org/geos/ in Java it would be JTS http://www.vividsolutions.com/jts/JTSHome.htm.

If, however, you are creating the grid yourself, and want an answer to the far simpler question of how do I find the intersection points between the horizontal and vertical lines of the grid I just created, the answer would be to use the same exact latitude, longitude values that you used for the LineStrings in a nested loop:

Document document = new Document();
for(y = -90; y < 0; y += 15){
  for(x = -180; x < 0; x+= 15){
     Point point = new Point();
     point.Coordinate = new Vector(x, y);
     Placemark placemark = new Placemark();
     placemark.Geometry = point;
     document.AddFeature(placemark);
  }
}

    .. repeat for the other 4 quadrants

// It's conventional for the root element to be Kml,
// but you could use document instead.
Kml root = new Kml();
root.Feature = document;
XmlFile kml = KmlFile.Create(root, false);

Here is some source code if you wanted to use DotSpatial, for instance, to find the intersection between the grids and a Shapefile. In this case, the shapefile has river lines and only produced one intersection point. Be advised that the topology intersection code is kind of slow, so you will want to use extent checking to speed things up. In your case you may want to build new Features by using KMLSharp to read a the linestring coordinates in a kml source file, rather than opening a shapefile, but the intersection code would be similar.

As a side note, I don't think the seemingly easy to use FeatureSet.Intersection method is smart enough to handle the case where line intersections produce point features as intersections. It only works for points or polygons where the output is likely to be the same feature type as the input.

    using DotSpatial.Controls;
    using DotSpatial.Data;
    using DotSpatial.Topology;
    using DotSpatial.Symbology;


    private FeatureSet gridLines;

    private void buttonAddGrid_Click(object sender, EventArgs e)
    {
        gridLines = new FeatureSet(FeatureType.Line);
        for (int x = -180; x < 0; x += 15)
        {
            List<Coordinate> coords = new List<Coordinate>();
            coords.Add(new Coordinate(x, -90));
            coords.Add(new Coordinate(x, 90));
            LineString ls = new LineString(coords);
            gridLines.AddFeature(ls);
        }
        for (int y = -90; y < 0; y += 15)
        {
            List<Coordinate> coords = new List<Coordinate>();
            coords.Add(new Coordinate(-180, y));
            coords.Add(new Coordinate(180, y));
            LineString ls = new LineString(coords);
            gridLines.AddFeature(ls);
        }

        map1.Layers.Add(new MapLineLayer(gridLines));
    }

    private void buttonIntersect_Click(object sender, EventArgs e)
    {
        if (gridLines == null)
        {
            MessageBox.Show("First add the grid.");
        }

        IFeatureSet river = FeatureSet.Open(@"C:\Data\Rivers\River.shp");
        MapLineLayer riverLayer = new MapLineLayer(river);
        map1.Layers.Add(river);




        List<DotSpatial.Topology.Point> allResultPoints = new List<DotSpatial.Topology.Point>();
        foreach (Feature polygon in river.Features)
        {
            Geometry lineString = polygon.BasicGeometry as Geometry;
            foreach (Feature lineFeature in gridLines.Features)
            {
                // Speed up calculation with extent testing.
                if(!lineFeature.Envelope.Intersects(lineString.Envelope)){
                    continue;
                }
                IFeature intersectFeature = lineFeature.Intersection(lineString);
                if (intersectFeature == null)
                {
                    continue;
                }


                MultiPoint multi = intersectFeature.BasicGeometry as MultiPoint;
                if (multi != null)
                {
                    for(int i = 0; i < multi.NumGeometries; i++)
                    {
                        allResultPoints.Add(intersectFeature.GetBasicGeometryN(i) as DotSpatial.Topology.Point);
                    }
                }
                DotSpatial.Topology.Point single = intersectFeature.BasicGeometry as DotSpatial.Topology.Point;
                {
                    allResultPoints.Add(single);
                }
            }
        }

        FeatureSet finalPoints = new FeatureSet(FeatureType.Point);
        foreach(DotSpatial.Topology.Point pt in allResultPoints){
            finalPoints.AddFeature(pt);
        }
        map1.Layers.Add(new MapPointLayer(finalPoints));
    }