< Summary

Class:Itinero.Geo.GeoExtensions
Assembly:Itinero
File(s):/home/runner/work/routing2/routing2/src/Itinero/Geo/GeoExtensions.cs
Covered lines:138
Uncovered lines:117
Coverable lines:255
Total lines:528
Line coverage:54.1% (138 of 255)
Covered branches:32
Total branches:88
Branch coverage:36.3% (32 of 88)
Tag:263_26948838820

Metrics

MethodBranch coverage Cyclomatic complexity Line coverage
DistanceEstimateInMeter(...)100%1100%
DistanceEstimateInMeterShape(...)100%4100%
DistanceEstimateInMeter(...)0%20%
OffsetWithDistanceX(...)100%1100%
OffsetWithDistanceY(...)100%1100%
PositionAlongLineInMeters(...)100%10%
PositionAlongLine(...)100%10%
PositionAlongLineInMeters(...)0%100%
PositionAlongLine(...)50%476.92%
ProjectOn(...)83.33%6100%
Center(...)50%470%
Expand(...)0%120%
Intersect(...)27.77%3650%
A(...)100%1100%
B(...)100%1100%
C(...)100%1100%
BoxAround(...)100%1100%
Overlaps(...)100%6100%
AngleWithMeridian(...)75%483.33%

File(s)

/home/runner/work/routing2/routing2/src/Itinero/Geo/GeoExtensions.cs

#LineLine coverage
 1using System;
 2using System.Collections.Generic;
 3using Itinero.Geo.Elevation;
 4
 5namespace Itinero.Geo;
 6
 7/// <summary>
 8/// Contains extension methods to work with coordinates, lines, bounding boxes and basic spatial operations.
 9/// </summary>
 10public static class GeoExtensions
 11{
 12    /// <summary>
 13    /// Returns an estimate of the distance between the two given coordinates.
 14    /// </summary>
 15    /// <param name="coordinate1">The first coordinate.</param>
 16    /// <param name="coordinate2">The second coordinate.</param>
 17    /// <remarks>Accuracy decreases with distance.</remarks>
 18    public static double DistanceEstimateInMeter(this (double longitude, double latitude, float? e) coordinate1,
 19        (double longitude, double latitude, float? e) coordinate2)
 192653520    {
 192653521        var lat1Rad = coordinate1.latitude / 180d * Math.PI;
 192653522        var lon1Rad = coordinate1.longitude / 180d * Math.PI;
 192653523        var lat2Rad = coordinate2.latitude / 180d * Math.PI;
 192653524        var lon2Rad = coordinate2.longitude / 180d * Math.PI;
 25
 192653526        var x = (lon2Rad - lon1Rad) * Math.Cos((lat1Rad + lat2Rad) / 2.0);
 192653527        var y = lat2Rad - lat1Rad;
 28
 192653529        var m = Math.Sqrt((x * x) + (y * y)) * Constants.RadiusOfEarth;
 30
 192653531        return m;
 192653532    }
 33
 34    internal static double DistanceEstimateInMeterShape(
 35        this (double longitude, double latitude, float? e) coordinate1,
 36        (double longitude, double latitude, float? e) coordinate2,
 37        IEnumerable<(double longitude, double latitude, float? e)>? shape = null)
 5067038    {
 5067039        if (shape == null)
 42640        {
 42641            return coordinate1.DistanceEstimateInMeter(coordinate2);
 42        }
 43
 5024444        var distance = 0.0;
 45
 5024446        using var shapeEnumerator = shape.GetEnumerator();
 5024447        var previous = coordinate1;
 48
 11662449        while (shapeEnumerator.MoveNext())
 6638050        {
 6638051            var current = shapeEnumerator.Current;
 6638052            distance += previous.DistanceEstimateInMeter(current);
 6638053            previous = current;
 6638054        }
 55
 5024456        distance += previous.DistanceEstimateInMeter(coordinate2);
 57
 5024458        return distance;
 5067059    }
 60
 61    /// <summary>
 62    /// Returns an estimate of the length of the given linestring.
 63    /// </summary>
 64    /// <param name="lineString">The linestring.</param>
 65    /// <remarks>Accuracy decreases with distance.</remarks>
 66    public static double DistanceEstimateInMeter(
 67        this IEnumerable<(double longitude, double latitude, float? e)> lineString)
 068    {
 069        var distance = 0.0;
 70
 071        using var shapeEnumerator = lineString.GetEnumerator();
 072        shapeEnumerator.MoveNext();
 073        var previous = shapeEnumerator.Current;
 74
 075        while (shapeEnumerator.MoveNext())
 076        {
 077            var current = shapeEnumerator.Current;
 078            distance += previous.DistanceEstimateInMeter(current);
 079            previous = current;
 080        }
 81
 082        return distance;
 083    }
 84
 85    /// <summary>
 86    /// Returns a coordinate offset with a given distance.
 87    /// </summary>
 88    /// <param name="coordinate">The coordinate.</param>
 89    /// <param name="meter">The distance.</param>
 90    /// <returns>An offset coordinate.</returns>
 91    public static (double longitude, double latitude, float? e) OffsetWithDistanceX(
 92        this (double longitude, double latitude, float? e) coordinate, double meter)
 3593    {
 94        const double offset = 0.001;
 3595        var offsetLon = (coordinate.longitude + offset, coordinate.latitude).AddElevation(coordinate.e);
 3596        var lonDistance = offsetLon.DistanceEstimateInMeter(coordinate);
 97
 3598        return (coordinate.longitude + (meter / lonDistance * offset),
 3599            coordinate.latitude).AddElevation(coordinate.e);
 35100    }
 101
 102    /// <summary>
 103    /// Returns a coordinate offset with a given distance.
 104    /// </summary>
 105    /// <param name="coordinate">The coordinate.</param>
 106    /// <param name="meter">The distance.</param>
 107    /// <returns>An offset coordinate.</returns>
 108    public static (double longitude, double latitude, float? e) OffsetWithDistanceY(
 109        this (double longitude, double latitude, float? e) coordinate,
 110        double meter)
 26111    {
 112        const double offset = 0.001;
 26113        var offsetLat = (coordinate.longitude, coordinate.latitude + offset).AddElevation(coordinate.e);
 26114        var latDistance = offsetLat.DistanceEstimateInMeter(coordinate);
 115
 26116        return (coordinate.longitude,
 26117            coordinate.latitude + (meter / latDistance * offset)).AddElevation(coordinate.e);
 26118    }
 119
 120    /// <summary>
 121    /// Calculates an offset position along the line segment.
 122    /// </summary>
 123    /// <param name="line">The line segment.</param>
 124    /// <param name="position">The position in meters relative to the start point.</param>
 125    /// <returns>The offset coordinate.</returns>
 126    public static (double longitude, double latitude, float? e) PositionAlongLineInMeters(
 127        this IEnumerable<(double longitude, double latitude, float? e)> line, double position)
 0128    {
 129        // ReSharper disable once PossibleMultipleEnumeration
 0130        var length = line.DistanceEstimateInMeter();
 131
 132        // ReSharper disable once PossibleMultipleEnumeration
 0133        return line.PositionAlongLineInMeters(length, position);
 0134    }
 135
 136    /// <summary>
 137    /// Calculates an offset position along the line segment.
 138    /// </summary>
 139    /// <param name="line">The line segment.</param>
 140    /// <param name="offset">The offset [0,1].</param>
 141    /// <returns>The offset coordinate.</returns>
 142    public static (double longitude, double latitude, float? e) PositionAlongLine(
 143        this IEnumerable<(double longitude, double latitude, float? e)> line, double offset)
 0144    {
 145        // ReSharper disable once PossibleMultipleEnumeration
 0146        var length = line.DistanceEstimateInMeter();
 0147        var targetLength = length * (offset / (double)ushort.MaxValue);
 148
 149        // ReSharper disable once PossibleMultipleEnumeration
 0150        return line.PositionAlongLineInMeters(length, targetLength);
 0151    }
 152
 153    private static (double longitude, double latitude, float? e) PositionAlongLineInMeters(
 154        this IEnumerable<(double longitude, double latitude, float? e)> line, double length, double targetLength)
 0155    {
 0156        var currentLength = 0.0;
 157
 158        // ReSharper disable once PossibleMultipleEnumeration
 0159        using var enumerator = line.GetEnumerator();
 0160        if (!enumerator.MoveNext()) throw new Exception("Line doesn't have 2 locations");
 0161        var previous = enumerator.Current;
 0162        while (enumerator.MoveNext())
 0163        {
 0164            var current = enumerator.Current;
 0165            var segmentLength = current.DistanceEstimateInMeter(previous);
 166
 167            // check if the target is in this segment or not.
 0168            if (segmentLength + currentLength < targetLength)
 0169            {
 0170                currentLength += segmentLength;
 0171                previous = current;
 0172                continue;
 173            }
 174
 0175            var segmentOffsetLength = segmentLength + currentLength - targetLength;
 0176            var segmentOffset = 1 - (segmentOffsetLength / segmentLength);
 177
 0178            float? e = null;
 0179            if (previous.e.HasValue &&
 0180                current.e.HasValue)
 0181            {
 0182                e = (float)(previous.e.Value + (segmentOffset * (current.e.Value - previous.e.Value)));
 0183            }
 184
 0185            return (previous.longitude + (segmentOffset * (current.longitude - previous.longitude)),
 0186                previous.latitude + (segmentOffset * (current.latitude - previous.latitude)), e);
 187        }
 188
 0189        return previous;
 0190    }
 191
 192    /// <summary>
 193    /// Calculates an offset position along the line segment.
 194    /// </summary>
 195    /// <param name="line">The line segment.</param>
 196    /// <param name="offset">The offset [0,1].</param>
 197    /// <returns>The offset coordinate.</returns>
 198    public static (double longitude, double latitude, float? e) PositionAlongLine(
 199        this ((double longitude, double latitude, float? e) coordinate1,
 200            (double longitude, double latitude, float? e) coordinate2) line, double offset)
 40201    {
 40202        var coordinate1 = line.coordinate1;
 40203        var coordinate2 = line.coordinate2;
 204
 40205        var latitude = coordinate1.latitude + ((coordinate2.latitude - coordinate1.latitude) * offset);
 40206        var longitude = coordinate1.longitude + ((coordinate2.longitude - coordinate1.longitude) * offset);
 40207        float? e = null;
 40208        if (coordinate1.e.HasValue &&
 40209            coordinate2.e.HasValue)
 0210        {
 0211            e = (float)(coordinate1.e.Value - ((coordinate2.e.Value - coordinate1.e.Value) * offset));
 0212        }
 213
 40214        return (longitude, latitude).AddElevation(e);
 40215    }
 216
 217    private const double E = 0.0000000001;
 218
 219    /// <summary>
 220    /// Projects <paramref name="coordinate"/> perpendicularly onto the segment
 221    /// <paramref name="line"/>. Returns the foot of the perpendicular if it
 222    /// falls on the segment; null otherwise (or if the segment is degenerate).
 223    ///
 224    /// <para>
 225    /// The result is direction-independent â€” calling with the line endpoints
 226    /// swapped produces the same projected coordinate (to within FP).
 227    /// </para>
 228    ///
 229    /// <para>
 230    /// We work in a local Cartesian frame in meters around the segment:
 231    /// latitude is scaled by a constant <c>~111320 m/°</c>, longitude by
 232    /// <c>cos(midpointLat) × 111320</c>. Using the segment midpoint as the
 233    /// reference latitude (rather than coordinate1) keeps the scaling
 234    /// symmetric under endpoint swap. The projection itself is the standard
 235    /// dot-product clamp:
 236    /// </para>
 237    /// <code>
 238    /// t = ((Q âˆ’ A) · (B âˆ’ A)) / ((B âˆ’ A) · (B âˆ’ A))
 239    /// foot = A + t × (B âˆ’ A)   when 0 â‰¤ t â‰¤ 1, else null
 240    /// </code>
 241    /// </summary>
 242    /// <param name="line">The segment.</param>
 243    /// <param name="coordinate">The point to project.</param>
 244    /// <returns>The foot of the perpendicular on the segment, or null if
 245    /// it falls outside the segment endpoints.</returns>
 246    public static (double longitude, double latitude, float? e)? ProjectOn(
 247        this ((double longitude, double latitude, float? e) coordinate1,
 248            (double longitude, double latitude, float? e) coordinate2) line,
 249        (double longitude, double latitude, float? e) coordinate)
 109214250    {
 109214251        var a = line.coordinate1;
 109214252        var b = line.coordinate2;
 253
 254        const double metersPerDegLat = 111320.0;
 109214255        var refLatRad = (a.latitude + b.latitude) * 0.5 * Math.PI / 180.0;
 109214256        var metersPerDegLon = metersPerDegLat * Math.Cos(refLatRad);
 257
 109214258        var abDx = (b.longitude - a.longitude) * metersPerDegLon;
 109214259        var abDy = (b.latitude - a.latitude) * metersPerDegLat;
 109214260        var ab2 = abDx * abDx + abDy * abDy;
 109214261        if (ab2 < E) return null;
 262
 109214263        var aqDx = (coordinate.longitude - a.longitude) * metersPerDegLon;
 109214264        var aqDy = (coordinate.latitude - a.latitude) * metersPerDegLat;
 109214265        var t = (abDx * aqDx + abDy * aqDy) / ab2;
 210225266        if (t < 0.0 || t > 1.0) return null;
 267
 8203268        return (a.longitude + t * (b.longitude - a.longitude),
 8203269                a.latitude + t * (b.latitude - a.latitude),
 8203270                (float?)null);
 109214271    }
 272
 273    /// <summary>
 274    /// Returns the center of the box.
 275    /// </summary>
 276    /// <param name="box">The box.</param>
 277    /// <returns>The center.</returns>
 278    public static (double longitude, double latitude, float? e) Center(
 279        this ((double longitude, double latitude, float? e) topLeft, (double longitude, double latitude, float? e)
 280            bottomRight) box)
 2210281    {
 2210282        float? e = null;
 2210283        if (box.topLeft.e.HasValue &&
 2210284            box.bottomRight.e.HasValue)
 0285        {
 0286            e = box.topLeft.e.Value + box.bottomRight.e.Value;
 0287        }
 288
 2210289        return ((box.topLeft.longitude + box.bottomRight.longitude) / 2,
 2210290            (box.topLeft.latitude + box.bottomRight.latitude) / 2).AddElevation(e);
 2210291    }
 292
 293    /// <summary>
 294    /// Expands the given box with the other box to encompass both.
 295    /// </summary>
 296    /// <param name="box">The original box.</param>
 297    /// <param name="other">The other box.</param>
 298    /// <returns>The expand box or the original box if the other was already contained.</returns>
 299    public static ((double longitude, double latitude, float? e) topLeft, (double longitude, double latitude, float?
 300        e) bottomRight)
 301        Expand(
 302            this ((double longitude, double latitude, float? e) topLeft, (double longitude, double latitude, float?
 303                e) bottomRight) box,
 304            ((double longitude, double latitude, float? e) topLeft, (double longitude, double latitude, float? e)
 305                bottomRight) other)
 0306    {
 0307        if (!box.Overlaps(other.topLeft))
 0308        {
 0309            var center = box.Center();
 310
 311            // handle left.
 0312            var left = box.topLeft.longitude;
 0313            if (!box.Overlaps((other.topLeft.longitude, center.latitude, null)))
 0314            {
 0315                left = other.topLeft.longitude;
 0316            }
 317
 318            // handle top.
 0319            var top = box.topLeft.latitude;
 0320            if (!box.Overlaps((center.longitude, other.topLeft.latitude, null)))
 0321            {
 0322                top = other.topLeft.latitude;
 0323            }
 324
 0325            box = ((left, top, null), box.bottomRight);
 0326        }
 327
 0328        if (!box.Overlaps(other.bottomRight))
 0329        {
 0330            var center = box.Center();
 331
 332            // handle right.
 0333            var right = box.bottomRight.longitude;
 0334            if (!box.Overlaps((other.bottomRight.longitude, center.latitude, null)))
 0335            {
 0336                right = other.bottomRight.longitude;
 0337            }
 338
 339            // handle bottom.
 0340            var bottom = box.bottomRight.latitude;
 0341            if (!box.Overlaps((center.longitude, other.bottomRight.latitude, null)))
 0342            {
 0343                bottom = other.bottomRight.latitude;
 0344            }
 345
 0346            box = (box.topLeft, (right, bottom, null));
 0347        }
 348
 0349        return box;
 0350    }
 351
 352    /// <summary>
 353    /// Calculates the intersection point of the given line with this line.
 354    ///
 355    /// Returns null if the lines have the same direction or don't intersect.
 356    ///
 357    /// Assumes the given line is not a segment and this line is a segment.
 358    /// </summary>
 359    public static (double longitude, double latitude, float? e)? Intersect(
 360        this ((double longitude, double latitude, float? e) coordinate1,
 361            (double longitude, double latitude, float? e) coordinate2) thisLine,
 362        ((double longitude, double latitude, float? e) coordinate1,
 363            (double longitude, double latitude, float? e) coordinate2) line, bool checkSegment = true)
 4364    {
 4365        var det = (double)((line.A() * thisLine.B()) - (thisLine.A() * line.B()));
 4366        if (Math.Abs(det) <= E)
 1367        {
 368            // lines are parallel; no intersections.
 1369            return null;
 370        }
 371        else
 3372        {
 373            // lines are not the same and not parallel so they will intersect.
 3374            var x = ((thisLine.B() * line.C()) - (line.B() * thisLine.C())) / det;
 3375            var y = ((line.A() * thisLine.C()) - (thisLine.A() * line.C())) / det;
 376
 3377            (double latitude, double longitude, float? e) coordinate = (x, y, (float?)null);
 378
 379            // check if the coordinate is on this line.
 3380            if (checkSegment)
 3381            {
 3382                var dist = (thisLine.A() * thisLine.A()) + (thisLine.B() * thisLine.B());
 3383                var line1 = (coordinate, thisLine.coordinate1);
 3384                var distTo1 = (line1.A() * line1.A()) + (line1.B() * line1.B());
 3385                if (distTo1 > dist)
 1386                {
 1387                    return null;
 388                }
 389
 2390                var line2 = (coordinate, thisLine.coordinate2);
 2391                var distTo2 = (line2.A() * line2.A()) + (line2.B() * line2.B());
 2392                if (distTo2 > dist)
 1393                {
 1394                    return null;
 395                }
 1396            }
 397
 1398            if (thisLine.coordinate1.e == null || thisLine.coordinate2.e == null)
 1399            {
 1400                return coordinate;
 401            }
 402
 0403            float? e = null;
 0404            if (Math.Abs(thisLine.coordinate1.e.Value - thisLine.coordinate2.e.Value) < E)
 0405            {
 406                // don't calculate anything, elevation is identical.
 0407                e = thisLine.coordinate1.e;
 0408            }
 0409            else if (Math.Abs(thisLine.A()) < E && Math.Abs(thisLine.B()) < E)
 0410            {
 411                // tiny segment, not stable to calculate offset
 0412                e = thisLine.coordinate1.e;
 0413            }
 414            else
 0415            {
 416                // calculate offset and calculate an estimate of the elevation.
 0417                if (Math.Abs(thisLine.A()) > Math.Abs(thisLine.B()))
 0418                {
 0419                    var diffLat = Math.Abs(thisLine.A());
 0420                    var diffLatIntersection = Math.Abs(coordinate.latitude - thisLine.coordinate1.latitude);
 421
 0422                    e = (float)(((thisLine.coordinate2.e - thisLine.coordinate1.e) *
 0423                                 (diffLatIntersection / diffLat)) +
 0424                                thisLine.coordinate1.e);
 0425                }
 426                else
 0427                {
 0428                    var diffLon = Math.Abs(thisLine.B());
 0429                    var diffLonIntersection = Math.Abs(coordinate.longitude - thisLine.coordinate1.longitude);
 430
 0431                    e = (float)(((thisLine.coordinate2.e - thisLine.coordinate1.e) *
 0432                                 (diffLonIntersection / diffLon)) +
 0433                                thisLine.coordinate1.e);
 0434                }
 0435            }
 436
 0437            return coordinate.AddElevation(e);
 438        }
 4439    }
 440
 441    private static double A(this ((double longitude, double latitude, float? e) coordinate1,
 442        (double longitude, double latitude, float? e) coordinate2) line)
 42443    {
 42444        return line.coordinate2.latitude - line.coordinate1.latitude;
 42445    }
 446
 447    private static double B(this ((double longitude, double latitude, float? e) coordinate1,
 448        (double longitude, double latitude, float? e) coordinate2) line)
 42449    {
 42450        return line.coordinate1.longitude - line.coordinate2.longitude;
 42451    }
 452
 453    private static double C(this ((double longitude, double latitude, float? e) coordinate1,
 454        (double longitude, double latitude, float? e) coordinate2) line)
 12455    {
 12456        return (line.A() * line.coordinate1.longitude) +
 12457               (line.B() * line.coordinate1.latitude);
 12458    }
 459
 460    /// <summary>
 461    /// Creates a box around this coordinate with width/height approximately the given size in meter.
 462    /// </summary>
 463    /// <param name="coordinate">The coordinate.</param>
 464    /// <param name="sizeInMeters">The size in meter.</param>
 465    /// <returns>The size in meter.</returns>
 466    public static ((double longitude, double latitude, float? e) topLeft, (double longitude, double latitude, float?
 467        e) bottomRight)
 468        BoxAround(this (double longitude, double latitude, float? e) coordinate, double sizeInMeters)
 13733469    {
 13733470        var offsetLat = (coordinate.longitude, coordinate.latitude + 0.1, (float?)null);
 13733471        var offsetLon = (coordinate.longitude + 0.1, coordinate.latitude, (float?)null);
 13733472        var latDistance = offsetLat.DistanceEstimateInMeter(coordinate);
 13733473        var lonDistance = offsetLon.DistanceEstimateInMeter(coordinate);
 474
 13733475        return ((coordinate.longitude - (sizeInMeters / lonDistance * 0.1),
 13733476                coordinate.latitude + (sizeInMeters / latDistance * 0.1), null),
 13733477            (coordinate.longitude + (sizeInMeters / lonDistance * 0.1),
 13733478                coordinate.latitude - (sizeInMeters / latDistance * 0.1), null));
 13733479    }
 480
 481    /// <summary>
 482    /// Returns true if the box overlaps the given coordinate.
 483    /// </summary>
 484    /// <param name="box">The box.</param>
 485    /// <param name="coordinate">The coordinate.</param>
 486    /// <returns>True of the coordinate is inside the bounding box.</returns>
 487    public static bool Overlaps(
 488        this ((double longitude, double latitude, float? e) topLeft, (double longitude, double latitude, float? e)
 489            bottomRight) box,
 490        (double longitude, double latitude, float? e) coordinate)
 2998145491    {
 2998145492        return box.bottomRight.latitude < coordinate.latitude && coordinate.latitude <= box.topLeft.latitude &&
 2998145493               box.topLeft.longitude < coordinate.longitude && coordinate.longitude <= box.bottomRight.longitude;
 2998145494    }
 495
 496    /// <summary>
 497    /// Given two WGS84 coordinates, if walking from c1 to c2, it gives the angle that one would be following.
 498    ///
 499    /// 0° is north, 90° is east, -90° is west, both 180 and -180 are south
 500    /// </summary>
 501    /// <param name="c1">The first coordinate.</param>
 502    /// <param name="c2">The second coordinate.</param>
 503    /// <returns>The angle with the meridian in Northern direction.</returns>
 504    public static double AngleWithMeridian(this (double longitude, double latitude, float? e) c1,
 505        (double longitude, double latitude, float? e) c2)
 212506    {
 212507        var dy = c2.latitude - c1.latitude;
 212508        var dx = Math.Cos(Math.PI / 180 * c1.latitude) * (c2.longitude - c1.longitude);
 509        // phi is the angle we search, but with 0 pointing eastwards and in radians
 212510        var phi = Math.Atan2(dy, dx);
 212511        var angle =
 212512            (phi - (Math.PI / 2)) // Rotate 90° to have the north up
 212513            * 180 / Math.PI; // Convert to degrees
 212514        angle = -angle;
 515        // A bit of normalization below:
 212516        if (angle < -180)
 0517        {
 0518            angle += 360;
 0519        }
 520
 212521        if (angle > 180)
 92522        {
 92523            angle -= 360;
 92524        }
 525
 212526        return angle;
 212527    }
 528}