Tronald / CoordinateSharp

A library designed to ease geographic coordinate format conversions, and determine sun/moon information in C#
Other
364 stars 59 forks source link

Determine bounding rectangle for two coordinate (a line on map) #130

Closed mohammad4x closed 4 years ago

mohammad4x commented 4 years ago

Hi. I have A and B locations in "lat, lon" format and the distance between A and A1 in "meter, kilometer,...". How can i calculate A1, A2, B1, B2 locations in "lat, lon" format? (image) i know it's not a question about using this library but i hope you could help me based on your knowledge of Geometry and Geo-Locations, also i used CoordinateSharp and Move() method for this problem but i hit a dead end, my spherical geometry petered out :)

Tronald commented 4 years ago

Hi Mohammad

As I guess you have seen, moving a coordinate on an initial bearing (using Haversine or Vincenty) won't give the best results when trying to create solid geofences/shapes, especially if you are covering great distances due to the shape of the earth.

Lastly, I will say the moving operations in CoordinateSharp are low precision and best used for estimating. You may need a stronger gis tool depending on what you are trying to achieve, but we may be able to make it work depending on what you need.

mohammad4x commented 4 years ago

1- not a fully perfect rectangle, but on a globe (I'm using Open Street Map) 2- distances vary too much, between hundreds of meters to maybe couple thousands kilometers and i need under 50m (it's better to be under 10m) error in drawing that rectangle 3- issues were about my lack of geometry knowledge, i found a solution at the end and wrote this code:

public class BoxFinderTest
    {
        private readonly GeoDistance distance;
        private readonly Coordinate start;
        private readonly Coordinate end;
        private readonly Distance routeDistance;

        public BoxFinderTest(GeoDistance distance, GeoPoint start, GeoPoint end)
        {
            this.distance = distance;
            var offEagerLoad = new EagerLoad(false);

            if (!Coordinate.TryParse(start.ToString(), offEagerLoad, out this.start))
                throw new ArgumentException("invalid start point!");

            if (!Coordinate.TryParse(end.ToString(), offEagerLoad, out this.end))
                throw new ArgumentException("invalid end point!");

            routeDistance = new Distance(this.start,this.end);
        }

        public IEnumerable<GeoPoint> GetBoundingBox()
        {
            return new List<GeoPoint>() { GetTopLeft(), GetTopRight(), GetBottomRight(), GetBottomLeft(), GetTopLeft() };
        }

        private GeoPoint GetBottomLeft()
        {
            Coordinate newPoint;
            if (this.routeDistance.Bearing == 0)
                newPoint = GetNewPoint(270, this.start);
            else if (this.routeDistance.Bearing > 180)
                newPoint = GetNewPoint(routeDistance.Bearing - 90, this.end);
            else if (this.routeDistance.Bearing < 180)
                newPoint = GetNewPoint(routeDistance.Bearing + 90, this.start);
            else
                newPoint = GetNewPoint(270, this.end);
            return new GeoPoint
            {
                Lat = newPoint.Latitude.ToDouble(),
                Lon = newPoint.Longitude.ToDouble()
            };
        }

        private GeoPoint GetBottomRight()
        {
            Coordinate newPoint;
            if (this.routeDistance.Bearing == 0)
                newPoint = GetNewPoint(90, this.start);
            else if (this.routeDistance.Bearing > 180)
                newPoint = GetNewPoint(routeDistance.Bearing - 90, this.start);
            else if (this.routeDistance.Bearing < 180)
                newPoint = GetNewPoint(routeDistance.Bearing + 90, this.end);
            else
                newPoint = GetNewPoint(90, this.end);
            return new GeoPoint
            {
                Lat = newPoint.Latitude.ToDouble(),
                Lon = newPoint.Longitude.ToDouble()
            };
        }

        private GeoPoint GetTopRight()
        {
            Coordinate newPoint;
            if (this.routeDistance.Bearing == 0)
                newPoint = GetNewPoint(90, this.end);
            else if (this.routeDistance.Bearing > 180)
                newPoint = GetNewPoint((routeDistance.Bearing + 90) % 360, this.start);
            else if (this.routeDistance.Bearing < 180)
                newPoint = GetNewPoint((routeDistance.Bearing - 90) % 360, this.end);
            else
                newPoint = GetNewPoint(90, this.start);
            return new GeoPoint
            {
                Lat = newPoint.Latitude.ToDouble(),
                Lon = newPoint.Longitude.ToDouble()
            };
        }

        private GeoPoint GetTopLeft()
        {
            Coordinate newPoint;
            if (this.routeDistance.Bearing == 0)
                newPoint = GetNewPoint(270, this.end);
            else if (this.routeDistance.Bearing > 180)
                newPoint = GetNewPoint((routeDistance.Bearing + 90) % 360, this.end);
            else if (this.routeDistance.Bearing < 180)
                newPoint = GetNewPoint((routeDistance.Bearing - 90) % 360, this.start);
            else
                newPoint = GetNewPoint(270, this.start);
            return new GeoPoint
            {
                Lat = newPoint.Latitude.ToDouble(),
                Lon = newPoint.Longitude.ToDouble()
            };
        }

        private Coordinate GetNewPoint(double bearing, Coordinate startPoint)
        {
            Coordinate newPoint = new Coordinate(startPoint.Latitude.ToDouble(),startPoint.Longitude.ToDouble());
            newPoint.Move(distance.Distance, bearing, Shape.Ellipsoid);
            return newPoint;
        }
    }

and i've tested it on some 2km lines, it seems almost accurate (under 50 meter error). if you have any suggestion to make this more accurate , i would be happy to take them :)

Tronald commented 4 years ago

Your code looks good to me. CoordinateSharp's precision in this area is frankly poor, as the requirement to dial it in further as not been expressed.

It admittedly bugs me that a 2km rectangle can't hit < 1m of precision though. This would be a good opportunity to re-examine this area of the library however to see if we can improve the formula's without sacrificing performance.

I'm going to research this a bit and get back to you.

Tronald commented 4 years ago

Ok, so reading through this a bit more and going through old notes, I think I have an idea as to what is occurring.

The move operation you are utilizing makes use of the Vincenty "direct" formula. I tested the formula against other tools and it is working correctly. Where the accuracy loss comes in is during the bearing change that occurs during the move operation I believe. Haversine (spherical earth) also has bearing changes that need to be considered if creating shapes.

For example:

If we start at 45° 01′ 04.7878″ N, 45 E and move 2 km with an initial bearing of 90°  we finish at 45° 01′ 04.7777″ N, 045° 01′ 31.3449″ E with a final bearing NOT facing 90° , but rather 090° 01′ 04.61″

I believe if we are to utilize drawing geofences we must take into consideration the final bearing of the line and offset the next bearing with the change. This should dial in the accuracy.

CoordinateSharp currently does not provide a final bearing property within the Distance class, but I think it should be added now that you have brought this scenario up.

An example of a Vincenty calculator can be found here

I will research and test this for the next release. My apologies however, that I cannot provide a more immediate solution within the scope of the current version of the library.

Keeping this open for tracking purposes.

mohammad4x commented 4 years ago

Thank you very much for your time :) i can live with the current CoordinateSharp precision for a while ( at least one month), after that either should find an accurately developed formula and calculate the distance by myself or use your new release (which will be pretty accurate i hope 💃 )

Tronald commented 4 years ago

You may see an example of adjusting for final bearing using CoordinateSharp's current version in this fiddle.

Accuracy is greatly increased with 100 km square being accurate withing 25 meters. A 2 km square is now accurate within millimeters vs meters. Precision loss still occurs and grows with distance, but it should help greatly at first glance.

mohammad4x commented 4 years ago

I ran your code and got some under millimeters error, https://github.com/Tronald/CoordinateSharp/issues/130#issuecomment-596065351 code is working with finding every four corner of the rectangle based on A or B image. From what i understand from your code, you are suggesting "i should move from A to A1 (actually draw a line) then from A1 to B1 then to B2 then A2 then A to complete the rectangle (find four corners) and it's precision will be in the millimeter territory." Did i get what you meant right?

Tronald commented 4 years ago

Exactly right (Yes)! As long as you are drawing in a continuous line, you can offset the bearing changes that occur in the formulas. This isn't perfect, but it should increase accuracy down to a few millimeters if working with small areas. Precision loss will increase with size still.

I haven't tested this heavily with actual map drawings though, so let me know how it behaves.

I am working on a new Draw class for the library that will run this formula for you and store each point in a continuous list for reference.

mohammad4x commented 4 years ago

Well, first i should confess that i was misreading the result the whole time and i wasted my time and yours, i apologize for that, then i'm gonna express the real results that i got from both codes.

I have tested them on actual map drawing. Me saying "and i've tested it on some 2km lines, it seems almost accurate (under 50 meter error)." was completely garbage.

maybe i translated your code badly, you can check these two codes out here in this fiddle. I compared two types of distances in my and your code. your suggested code behave oddly when the distance grows! 65km error. I don't know, i say it again mayne i translated it wrong.

Tronald commented 4 years ago

Lol @ your string interpolation comment.

A couple of things. I didn't realize your rectangle was being drawn with 6 points (A1-A-A2 / B1-B-B2), my bad. Since this is the case your would need to draw a 6 point rectangle like so. A -> A1 -> B1 -> B -> B2 -> A2 -> A if using my example. My formula requires a continuous drawing.

You could also do A2->A1->B1->B2->A2 and then grab A and B by traveling half the distance between A2->A1 and B2->B1.

Regarding the code, the test results from my method are showing poor, because the bearing's aren't being calculated correctly on top of the non-continuous drawing.

In short, you are using

//This shows start to end
var lastBearing = new Distance(start, end).Bearing;
var newBearing = (lastBearing + 90) % 360;

This is simply going to adjust on the initial bearing, not the final bearing of the last movement. It should look like this when the steps are broken down. My code showed this in 2 lines, but I'm going to expand for reference as to what is occurring.

//Distance must be grabbed in reverse (end to start) to calculate bearing shift that occurs slowly during travel over distance.
//Also ensure you are specifying 'Ellipsoid' in all distance calculations as the default is `Sphere`
var lastDistanceReversed = new Distance(end, start, Shape.Ellipsoid);

//Flip the bearing direction (180 degrees) and normalize to obtain final bearing of the last movement
var finalBearing = (lastDistanceReversed.Bearing + 180) % 360;

//Determine how much the final bearing has shifted from the initial bearing
var shift = finalBearing-bearing;

//Factor in the shift to the initial bearing (bearing), add/subtract `X` degrees (90 in this example) of rotation and normalize
var newBearing = (shift + bearing + 90) % 360;

the newBearing is what you would determine your next point with while drawing a continuous rectangle. I tested the precision of your coordinate and long precision (5 km) and it's within mm's using my formula if drawing A2->A1->B1->B2->A2. Precision is calculated by finding the difference between the first A2 point and the last A2 point.

With all that said, you definitely did not waste my time. I'm glad you found a working solution that meets your precision requirements, but regardless it was good to revisit this part of the library and improve it. I'm excited that we will have a solution for this situation in the next release that should make shape drawing really easy and accurate.

This example draws a square by rotating the bearing 90 degrees with every line that will be available next release. The points are then stored in the GeoFence_Drawer.Points property for reference/iteration.

Coordinate c = new Coordinate(22.5, -72.3, new EagerLoad(false));

GeoFence.Drawer gd = new GeoFence.Drawer(c, Shape.Ellipsoid, 0);
gd.Draw(new Distance(5), 0); //First line no bearing change
gd.Draw(new Distance(5), 90); //turn right 90
gd.Draw(new Distance(5), 90); //turn right 90
gd.Close(); //draw line to starting point

gd.Points; //reference all drawn points. Iterate this.
mohammad4x commented 4 years ago

Thank you for your effort and this great library, i don't have anything new to say, you can close this issue :)