| | | 1 | | using System; |
| | | 2 | | using System.Collections.Generic; |
| | | 3 | | using Itinero.Geo.Elevation; |
| | | 4 | | |
| | | 5 | | namespace Itinero.Geo; |
| | | 6 | | |
| | | 7 | | /// <summary> |
| | | 8 | | /// Contains extension methods to work with coordinates, lines, bounding boxes and basic spatial operations. |
| | | 9 | | /// </summary> |
| | | 10 | | public 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) |
| | 1926535 | 20 | | { |
| | 1926535 | 21 | | var lat1Rad = coordinate1.latitude / 180d * Math.PI; |
| | 1926535 | 22 | | var lon1Rad = coordinate1.longitude / 180d * Math.PI; |
| | 1926535 | 23 | | var lat2Rad = coordinate2.latitude / 180d * Math.PI; |
| | 1926535 | 24 | | var lon2Rad = coordinate2.longitude / 180d * Math.PI; |
| | | 25 | | |
| | 1926535 | 26 | | var x = (lon2Rad - lon1Rad) * Math.Cos((lat1Rad + lat2Rad) / 2.0); |
| | 1926535 | 27 | | var y = lat2Rad - lat1Rad; |
| | | 28 | | |
| | 1926535 | 29 | | var m = Math.Sqrt((x * x) + (y * y)) * Constants.RadiusOfEarth; |
| | | 30 | | |
| | 1926535 | 31 | | return m; |
| | 1926535 | 32 | | } |
| | | 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) |
| | 50670 | 38 | | { |
| | 50670 | 39 | | if (shape == null) |
| | 426 | 40 | | { |
| | 426 | 41 | | return coordinate1.DistanceEstimateInMeter(coordinate2); |
| | | 42 | | } |
| | | 43 | | |
| | 50244 | 44 | | var distance = 0.0; |
| | | 45 | | |
| | 50244 | 46 | | using var shapeEnumerator = shape.GetEnumerator(); |
| | 50244 | 47 | | var previous = coordinate1; |
| | | 48 | | |
| | 116624 | 49 | | while (shapeEnumerator.MoveNext()) |
| | 66380 | 50 | | { |
| | 66380 | 51 | | var current = shapeEnumerator.Current; |
| | 66380 | 52 | | distance += previous.DistanceEstimateInMeter(current); |
| | 66380 | 53 | | previous = current; |
| | 66380 | 54 | | } |
| | | 55 | | |
| | 50244 | 56 | | distance += previous.DistanceEstimateInMeter(coordinate2); |
| | | 57 | | |
| | 50244 | 58 | | return distance; |
| | 50670 | 59 | | } |
| | | 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) |
| | 0 | 68 | | { |
| | 0 | 69 | | var distance = 0.0; |
| | | 70 | | |
| | 0 | 71 | | using var shapeEnumerator = lineString.GetEnumerator(); |
| | 0 | 72 | | shapeEnumerator.MoveNext(); |
| | 0 | 73 | | var previous = shapeEnumerator.Current; |
| | | 74 | | |
| | 0 | 75 | | while (shapeEnumerator.MoveNext()) |
| | 0 | 76 | | { |
| | 0 | 77 | | var current = shapeEnumerator.Current; |
| | 0 | 78 | | distance += previous.DistanceEstimateInMeter(current); |
| | 0 | 79 | | previous = current; |
| | 0 | 80 | | } |
| | | 81 | | |
| | 0 | 82 | | return distance; |
| | 0 | 83 | | } |
| | | 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) |
| | 35 | 93 | | { |
| | | 94 | | const double offset = 0.001; |
| | 35 | 95 | | var offsetLon = (coordinate.longitude + offset, coordinate.latitude).AddElevation(coordinate.e); |
| | 35 | 96 | | var lonDistance = offsetLon.DistanceEstimateInMeter(coordinate); |
| | | 97 | | |
| | 35 | 98 | | return (coordinate.longitude + (meter / lonDistance * offset), |
| | 35 | 99 | | coordinate.latitude).AddElevation(coordinate.e); |
| | 35 | 100 | | } |
| | | 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) |
| | 26 | 111 | | { |
| | | 112 | | const double offset = 0.001; |
| | 26 | 113 | | var offsetLat = (coordinate.longitude, coordinate.latitude + offset).AddElevation(coordinate.e); |
| | 26 | 114 | | var latDistance = offsetLat.DistanceEstimateInMeter(coordinate); |
| | | 115 | | |
| | 26 | 116 | | return (coordinate.longitude, |
| | 26 | 117 | | coordinate.latitude + (meter / latDistance * offset)).AddElevation(coordinate.e); |
| | 26 | 118 | | } |
| | | 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) |
| | 0 | 128 | | { |
| | | 129 | | // ReSharper disable once PossibleMultipleEnumeration |
| | 0 | 130 | | var length = line.DistanceEstimateInMeter(); |
| | | 131 | | |
| | | 132 | | // ReSharper disable once PossibleMultipleEnumeration |
| | 0 | 133 | | return line.PositionAlongLineInMeters(length, position); |
| | 0 | 134 | | } |
| | | 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) |
| | 0 | 144 | | { |
| | | 145 | | // ReSharper disable once PossibleMultipleEnumeration |
| | 0 | 146 | | var length = line.DistanceEstimateInMeter(); |
| | 0 | 147 | | var targetLength = length * (offset / (double)ushort.MaxValue); |
| | | 148 | | |
| | | 149 | | // ReSharper disable once PossibleMultipleEnumeration |
| | 0 | 150 | | return line.PositionAlongLineInMeters(length, targetLength); |
| | 0 | 151 | | } |
| | | 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) |
| | 0 | 155 | | { |
| | 0 | 156 | | var currentLength = 0.0; |
| | | 157 | | |
| | | 158 | | // ReSharper disable once PossibleMultipleEnumeration |
| | 0 | 159 | | using var enumerator = line.GetEnumerator(); |
| | 0 | 160 | | if (!enumerator.MoveNext()) throw new Exception("Line doesn't have 2 locations"); |
| | 0 | 161 | | var previous = enumerator.Current; |
| | 0 | 162 | | while (enumerator.MoveNext()) |
| | 0 | 163 | | { |
| | 0 | 164 | | var current = enumerator.Current; |
| | 0 | 165 | | var segmentLength = current.DistanceEstimateInMeter(previous); |
| | | 166 | | |
| | | 167 | | // check if the target is in this segment or not. |
| | 0 | 168 | | if (segmentLength + currentLength < targetLength) |
| | 0 | 169 | | { |
| | 0 | 170 | | currentLength += segmentLength; |
| | 0 | 171 | | previous = current; |
| | 0 | 172 | | continue; |
| | | 173 | | } |
| | | 174 | | |
| | 0 | 175 | | var segmentOffsetLength = segmentLength + currentLength - targetLength; |
| | 0 | 176 | | var segmentOffset = 1 - (segmentOffsetLength / segmentLength); |
| | | 177 | | |
| | 0 | 178 | | float? e = null; |
| | 0 | 179 | | if (previous.e.HasValue && |
| | 0 | 180 | | current.e.HasValue) |
| | 0 | 181 | | { |
| | 0 | 182 | | e = (float)(previous.e.Value + (segmentOffset * (current.e.Value - previous.e.Value))); |
| | 0 | 183 | | } |
| | | 184 | | |
| | 0 | 185 | | return (previous.longitude + (segmentOffset * (current.longitude - previous.longitude)), |
| | 0 | 186 | | previous.latitude + (segmentOffset * (current.latitude - previous.latitude)), e); |
| | | 187 | | } |
| | | 188 | | |
| | 0 | 189 | | return previous; |
| | 0 | 190 | | } |
| | | 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) |
| | 40 | 201 | | { |
| | 40 | 202 | | var coordinate1 = line.coordinate1; |
| | 40 | 203 | | var coordinate2 = line.coordinate2; |
| | | 204 | | |
| | 40 | 205 | | var latitude = coordinate1.latitude + ((coordinate2.latitude - coordinate1.latitude) * offset); |
| | 40 | 206 | | var longitude = coordinate1.longitude + ((coordinate2.longitude - coordinate1.longitude) * offset); |
| | 40 | 207 | | float? e = null; |
| | 40 | 208 | | if (coordinate1.e.HasValue && |
| | 40 | 209 | | coordinate2.e.HasValue) |
| | 0 | 210 | | { |
| | 0 | 211 | | e = (float)(coordinate1.e.Value - ((coordinate2.e.Value - coordinate1.e.Value) * offset)); |
| | 0 | 212 | | } |
| | | 213 | | |
| | 40 | 214 | | return (longitude, latitude).AddElevation(e); |
| | 40 | 215 | | } |
| | | 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) |
| | 109214 | 250 | | { |
| | 109214 | 251 | | var a = line.coordinate1; |
| | 109214 | 252 | | var b = line.coordinate2; |
| | | 253 | | |
| | | 254 | | const double metersPerDegLat = 111320.0; |
| | 109214 | 255 | | var refLatRad = (a.latitude + b.latitude) * 0.5 * Math.PI / 180.0; |
| | 109214 | 256 | | var metersPerDegLon = metersPerDegLat * Math.Cos(refLatRad); |
| | | 257 | | |
| | 109214 | 258 | | var abDx = (b.longitude - a.longitude) * metersPerDegLon; |
| | 109214 | 259 | | var abDy = (b.latitude - a.latitude) * metersPerDegLat; |
| | 109214 | 260 | | var ab2 = abDx * abDx + abDy * abDy; |
| | 109214 | 261 | | if (ab2 < E) return null; |
| | | 262 | | |
| | 109214 | 263 | | var aqDx = (coordinate.longitude - a.longitude) * metersPerDegLon; |
| | 109214 | 264 | | var aqDy = (coordinate.latitude - a.latitude) * metersPerDegLat; |
| | 109214 | 265 | | var t = (abDx * aqDx + abDy * aqDy) / ab2; |
| | 210225 | 266 | | if (t < 0.0 || t > 1.0) return null; |
| | | 267 | | |
| | 8203 | 268 | | return (a.longitude + t * (b.longitude - a.longitude), |
| | 8203 | 269 | | a.latitude + t * (b.latitude - a.latitude), |
| | 8203 | 270 | | (float?)null); |
| | 109214 | 271 | | } |
| | | 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) |
| | 2210 | 281 | | { |
| | 2210 | 282 | | float? e = null; |
| | 2210 | 283 | | if (box.topLeft.e.HasValue && |
| | 2210 | 284 | | box.bottomRight.e.HasValue) |
| | 0 | 285 | | { |
| | 0 | 286 | | e = box.topLeft.e.Value + box.bottomRight.e.Value; |
| | 0 | 287 | | } |
| | | 288 | | |
| | 2210 | 289 | | return ((box.topLeft.longitude + box.bottomRight.longitude) / 2, |
| | 2210 | 290 | | (box.topLeft.latitude + box.bottomRight.latitude) / 2).AddElevation(e); |
| | 2210 | 291 | | } |
| | | 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) |
| | 0 | 306 | | { |
| | 0 | 307 | | if (!box.Overlaps(other.topLeft)) |
| | 0 | 308 | | { |
| | 0 | 309 | | var center = box.Center(); |
| | | 310 | | |
| | | 311 | | // handle left. |
| | 0 | 312 | | var left = box.topLeft.longitude; |
| | 0 | 313 | | if (!box.Overlaps((other.topLeft.longitude, center.latitude, null))) |
| | 0 | 314 | | { |
| | 0 | 315 | | left = other.topLeft.longitude; |
| | 0 | 316 | | } |
| | | 317 | | |
| | | 318 | | // handle top. |
| | 0 | 319 | | var top = box.topLeft.latitude; |
| | 0 | 320 | | if (!box.Overlaps((center.longitude, other.topLeft.latitude, null))) |
| | 0 | 321 | | { |
| | 0 | 322 | | top = other.topLeft.latitude; |
| | 0 | 323 | | } |
| | | 324 | | |
| | 0 | 325 | | box = ((left, top, null), box.bottomRight); |
| | 0 | 326 | | } |
| | | 327 | | |
| | 0 | 328 | | if (!box.Overlaps(other.bottomRight)) |
| | 0 | 329 | | { |
| | 0 | 330 | | var center = box.Center(); |
| | | 331 | | |
| | | 332 | | // handle right. |
| | 0 | 333 | | var right = box.bottomRight.longitude; |
| | 0 | 334 | | if (!box.Overlaps((other.bottomRight.longitude, center.latitude, null))) |
| | 0 | 335 | | { |
| | 0 | 336 | | right = other.bottomRight.longitude; |
| | 0 | 337 | | } |
| | | 338 | | |
| | | 339 | | // handle bottom. |
| | 0 | 340 | | var bottom = box.bottomRight.latitude; |
| | 0 | 341 | | if (!box.Overlaps((center.longitude, other.bottomRight.latitude, null))) |
| | 0 | 342 | | { |
| | 0 | 343 | | bottom = other.bottomRight.latitude; |
| | 0 | 344 | | } |
| | | 345 | | |
| | 0 | 346 | | box = (box.topLeft, (right, bottom, null)); |
| | 0 | 347 | | } |
| | | 348 | | |
| | 0 | 349 | | return box; |
| | 0 | 350 | | } |
| | | 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) |
| | 4 | 364 | | { |
| | 4 | 365 | | var det = (double)((line.A() * thisLine.B()) - (thisLine.A() * line.B())); |
| | 4 | 366 | | if (Math.Abs(det) <= E) |
| | 1 | 367 | | { |
| | | 368 | | // lines are parallel; no intersections. |
| | 1 | 369 | | return null; |
| | | 370 | | } |
| | | 371 | | else |
| | 3 | 372 | | { |
| | | 373 | | // lines are not the same and not parallel so they will intersect. |
| | 3 | 374 | | var x = ((thisLine.B() * line.C()) - (line.B() * thisLine.C())) / det; |
| | 3 | 375 | | var y = ((line.A() * thisLine.C()) - (thisLine.A() * line.C())) / det; |
| | | 376 | | |
| | 3 | 377 | | (double latitude, double longitude, float? e) coordinate = (x, y, (float?)null); |
| | | 378 | | |
| | | 379 | | // check if the coordinate is on this line. |
| | 3 | 380 | | if (checkSegment) |
| | 3 | 381 | | { |
| | 3 | 382 | | var dist = (thisLine.A() * thisLine.A()) + (thisLine.B() * thisLine.B()); |
| | 3 | 383 | | var line1 = (coordinate, thisLine.coordinate1); |
| | 3 | 384 | | var distTo1 = (line1.A() * line1.A()) + (line1.B() * line1.B()); |
| | 3 | 385 | | if (distTo1 > dist) |
| | 1 | 386 | | { |
| | 1 | 387 | | return null; |
| | | 388 | | } |
| | | 389 | | |
| | 2 | 390 | | var line2 = (coordinate, thisLine.coordinate2); |
| | 2 | 391 | | var distTo2 = (line2.A() * line2.A()) + (line2.B() * line2.B()); |
| | 2 | 392 | | if (distTo2 > dist) |
| | 1 | 393 | | { |
| | 1 | 394 | | return null; |
| | | 395 | | } |
| | 1 | 396 | | } |
| | | 397 | | |
| | 1 | 398 | | if (thisLine.coordinate1.e == null || thisLine.coordinate2.e == null) |
| | 1 | 399 | | { |
| | 1 | 400 | | return coordinate; |
| | | 401 | | } |
| | | 402 | | |
| | 0 | 403 | | float? e = null; |
| | 0 | 404 | | if (Math.Abs(thisLine.coordinate1.e.Value - thisLine.coordinate2.e.Value) < E) |
| | 0 | 405 | | { |
| | | 406 | | // don't calculate anything, elevation is identical. |
| | 0 | 407 | | e = thisLine.coordinate1.e; |
| | 0 | 408 | | } |
| | 0 | 409 | | else if (Math.Abs(thisLine.A()) < E && Math.Abs(thisLine.B()) < E) |
| | 0 | 410 | | { |
| | | 411 | | // tiny segment, not stable to calculate offset |
| | 0 | 412 | | e = thisLine.coordinate1.e; |
| | 0 | 413 | | } |
| | | 414 | | else |
| | 0 | 415 | | { |
| | | 416 | | // calculate offset and calculate an estimate of the elevation. |
| | 0 | 417 | | if (Math.Abs(thisLine.A()) > Math.Abs(thisLine.B())) |
| | 0 | 418 | | { |
| | 0 | 419 | | var diffLat = Math.Abs(thisLine.A()); |
| | 0 | 420 | | var diffLatIntersection = Math.Abs(coordinate.latitude - thisLine.coordinate1.latitude); |
| | | 421 | | |
| | 0 | 422 | | e = (float)(((thisLine.coordinate2.e - thisLine.coordinate1.e) * |
| | 0 | 423 | | (diffLatIntersection / diffLat)) + |
| | 0 | 424 | | thisLine.coordinate1.e); |
| | 0 | 425 | | } |
| | | 426 | | else |
| | 0 | 427 | | { |
| | 0 | 428 | | var diffLon = Math.Abs(thisLine.B()); |
| | 0 | 429 | | var diffLonIntersection = Math.Abs(coordinate.longitude - thisLine.coordinate1.longitude); |
| | | 430 | | |
| | 0 | 431 | | e = (float)(((thisLine.coordinate2.e - thisLine.coordinate1.e) * |
| | 0 | 432 | | (diffLonIntersection / diffLon)) + |
| | 0 | 433 | | thisLine.coordinate1.e); |
| | 0 | 434 | | } |
| | 0 | 435 | | } |
| | | 436 | | |
| | 0 | 437 | | return coordinate.AddElevation(e); |
| | | 438 | | } |
| | 4 | 439 | | } |
| | | 440 | | |
| | | 441 | | private static double A(this ((double longitude, double latitude, float? e) coordinate1, |
| | | 442 | | (double longitude, double latitude, float? e) coordinate2) line) |
| | 42 | 443 | | { |
| | 42 | 444 | | return line.coordinate2.latitude - line.coordinate1.latitude; |
| | 42 | 445 | | } |
| | | 446 | | |
| | | 447 | | private static double B(this ((double longitude, double latitude, float? e) coordinate1, |
| | | 448 | | (double longitude, double latitude, float? e) coordinate2) line) |
| | 42 | 449 | | { |
| | 42 | 450 | | return line.coordinate1.longitude - line.coordinate2.longitude; |
| | 42 | 451 | | } |
| | | 452 | | |
| | | 453 | | private static double C(this ((double longitude, double latitude, float? e) coordinate1, |
| | | 454 | | (double longitude, double latitude, float? e) coordinate2) line) |
| | 12 | 455 | | { |
| | 12 | 456 | | return (line.A() * line.coordinate1.longitude) + |
| | 12 | 457 | | (line.B() * line.coordinate1.latitude); |
| | 12 | 458 | | } |
| | | 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) |
| | 13733 | 469 | | { |
| | 13733 | 470 | | var offsetLat = (coordinate.longitude, coordinate.latitude + 0.1, (float?)null); |
| | 13733 | 471 | | var offsetLon = (coordinate.longitude + 0.1, coordinate.latitude, (float?)null); |
| | 13733 | 472 | | var latDistance = offsetLat.DistanceEstimateInMeter(coordinate); |
| | 13733 | 473 | | var lonDistance = offsetLon.DistanceEstimateInMeter(coordinate); |
| | | 474 | | |
| | 13733 | 475 | | return ((coordinate.longitude - (sizeInMeters / lonDistance * 0.1), |
| | 13733 | 476 | | coordinate.latitude + (sizeInMeters / latDistance * 0.1), null), |
| | 13733 | 477 | | (coordinate.longitude + (sizeInMeters / lonDistance * 0.1), |
| | 13733 | 478 | | coordinate.latitude - (sizeInMeters / latDistance * 0.1), null)); |
| | 13733 | 479 | | } |
| | | 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) |
| | 2998145 | 491 | | { |
| | 2998145 | 492 | | return box.bottomRight.latitude < coordinate.latitude && coordinate.latitude <= box.topLeft.latitude && |
| | 2998145 | 493 | | box.topLeft.longitude < coordinate.longitude && coordinate.longitude <= box.bottomRight.longitude; |
| | 2998145 | 494 | | } |
| | | 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) |
| | 212 | 506 | | { |
| | 212 | 507 | | var dy = c2.latitude - c1.latitude; |
| | 212 | 508 | | 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 |
| | 212 | 510 | | var phi = Math.Atan2(dy, dx); |
| | 212 | 511 | | var angle = |
| | 212 | 512 | | (phi - (Math.PI / 2)) // Rotate 90° to have the north up |
| | 212 | 513 | | * 180 / Math.PI; // Convert to degrees |
| | 212 | 514 | | angle = -angle; |
| | | 515 | | // A bit of normalization below: |
| | 212 | 516 | | if (angle < -180) |
| | 0 | 517 | | { |
| | 0 | 518 | | angle += 360; |
| | 0 | 519 | | } |
| | | 520 | | |
| | 212 | 521 | | if (angle > 180) |
| | 92 | 522 | | { |
| | 92 | 523 | | angle -= 360; |
| | 92 | 524 | | } |
| | | 525 | | |
| | 212 | 526 | | return angle; |
| | 212 | 527 | | } |
| | | 528 | | } |