A simple utility class to handle dateline/anti-meridian issues with geometries

DatelineUtility.cs

using NetTopologySuite.Geometries;
using NetTopologySuite.IO;
using NUnit.Framework;
using System;
using System.Collections.Generic;

namespace NTS.Package.Tests
{
    /// <summary>
    /// Utility class for dateline handling (IDEA)
    /// </summary>
    public class DatelineUtility
    {
        /// <summary>
        /// Maximum value for <see cref="DateLineGap"/>.
        /// </summary>
        public const double MaxDateLineGap = 0.001;

        private static double _dateLineGap;

        /// <summary>
        /// Tests if a geometry might possibly cross the dateline
        /// </summary>
        /// <param name="geometry">The geometry to test</param>
        /// <returns></returns>
        /// <exception cref="ArgumentNullException">Thrown, if <paramref name="geometry"/> is <c>null</c></exception>
        /// <exception cref="ArgumentException">Thrown, if geometry is not geographic (see: <see cref="IsGeographic(Geometry)"/>)</exception>
        public static bool IsCrossingDateline(Geometry geometry)
        {
            // Argument checking
            if (geometry == null)
                throw new ArgumentNullException(nameof(geometry));

            // If we don't have a geographic coordinate system 
            if (!IsGeographic(geometry))
                throw new ArgumentException("Not a geography", nameof(geometry));

            // Can't do anything with empty geometries
            if (geometry.IsEmpty)
                return false;

            // If we have a geometry that is only consisting of points we don't need to do any
            // dateline handling
            if (geometry is IPuntal)
                return false;

            // if the bounds don't cross the dateline, the geometry doesn't either
            var env = geometry.EnvelopeInternal;
            if (env.MinX > -180d && env.MaxX < 180)
                return false;

            // Build filter that tests if we need to perform special dateline handling
            var cdf = new IsCrossingDatelineFilter();
            geometry.Apply(cdf);

            if (!cdf.IsCrossingDateline && geometry.EnvelopeInternal.MaxX <= 180)
                return false;

            return true;
        }

        /// <summary>
        /// Gets or sets a value 
        /// </summary>
        /// <remarks>
        /// The allowed range is [0, <see cref="MaxDateLineGap"/>]
        /// </remarks>
        public static double DateLineGap
        {
            get => _dateLineGap; 
            set
            {
                if (value < 0 || value > MaxDateLineGap)
                    throw new ArgumentOutOfRangeException(nameof(value), $"Not in allowed range [0 ... {MaxDateLineGap}]");

                _dateLineGap = value;
            }
        }

        /// <summary>
        /// Unwraps a geometry at the dateline
        /// </summary>
        /// <param name="geometry">The geometry to unwrap</param>
        /// <returns>The unwrapped geometry</returns>
        /// <exception cref="ArgumentException">Thrown, if geometry is not geographic (see: <see cref="IsGeographic(Geometry)"/>)</exception>
        public static Geometry UnwrapAtDateline(Geometry geometry, bool checkIfCrossing = true, bool unionResult = true, double dateLineGap = double.NaN, PrecisionModel pm = null)
        {
            // Does client want us to check if geometry crosses the dateline
            if (checkIfCrossing)
            {
                // Test if the geometry has a chance of crossing the dateline
                if (!IsCrossingDateline(geometry))
                    return geometry;
            }
            else
            {
                // Check if the geometry is at least a geographic one!
                if (!IsGeographic(geometry))
                    throw new ArgumentException(nameof(geometry));
            }

            // Get number of shifts
            double numberOfFullShifts = -Math.Floor((geometry.EnvelopeInternal.MinX + 180) / 360);
            if (numberOfFullShifts > 0)
            {
                // Move whole geometry by numberOfFullShifts * 360°
                geometry.Apply(new TranslateXFilter(numberOfFullShifts * 360d));
            }

            // Build a list of geometries
            var parts = new List<Geometry>();

            // Adjust gap at dateline
            if (double.IsNaN(dateLineGap))
                dateLineGap = DateLineGap;

            // Test dateline gap value
            if (dateLineGap < 0 || dateLineGap > MaxDateLineGap)
                throw new ArgumentOutOfRangeException(nameof(dateLineGap), $"Not in allowed range [0 ... {MaxDateLineGap}]");

            // While input geometry is not empty compute intersection with world bounds and add to list
            // Renove covered part from input geometry
            while (!geometry.IsEmpty)
            {
                // Intersect geometry with default world bounds
                var box = geometry.Factory.ToGeometry(new Envelope(-180, 180 - dateLineGap, -90, 90));
                var geom0 = geometry.Intersection(box);
                // Add part
                parts.Add(geom0);
                // Get remaining geometry and translate into world bounds
                geometry = geometry.Difference(box);
                geometry.Apply(new TranslateXFilter(-360d));
            }

            // check if the precision model is scaled

            // Union the geometries and return result
            if (unionResult)
            {
                // Check that the precision model is fixed
                if (pm != null && pm.IsFloating)
                    throw new ArgumentException("Not a fixed precision model", nameof(pm));

                return NetTopologySuite.Operation.OverlayNG.UnaryUnionNG.Union(parts, pm ?? geometry.Factory.PrecisionModel);
            }

            // Build a geometry
            return geometry.Factory.BuildGeometry(parts);
        }

        /// <summary>
        /// Tests if the geometry was build for geographic geometries
        /// </summary>
        /// <param name="geometry">The geometry</param>
        /// <returns><c>true</c> if it does.</returns>
        private static bool IsGeographic(Geometry geometry)
         => _geographicCrs.Value.Contains(geometry.SRID);        

        private static Lazy<HashSet<int>> _geographicCrs = new Lazy<HashSet<int>>(ReadGeographicCrs());

        private static HashSet<int> ReadGeographicCrs()
        {
            var hashSet = new HashSet<int>();
            hashSet.Add(4326);

            var asm = System.Reflection.Assembly.GetExecutingAssembly();
            using var strm = asm.GetManifestResourceStream($"{typeof(DatelineUtility).Namespace}.GeographicCrs.txt");
            using var sr = new System.IO.StreamReader(strm);

            while (!sr.EndOfStream)
            {
                string line = sr.ReadLine();
                if (string.IsNullOrWhiteSpace(line)) continue;
                if (line.StartsWith("#")) continue;

                string[] parts = line.Split(',');
                foreach(var part in parts)
                {
                    if (string.IsNullOrWhiteSpace(part))
                        continue;
                    if (int.TryParse(part, out int srid))
                        hashSet.Add(srid);
                }
            }
            return hashSet;
        }

        private class IsCrossingDatelineFilter : IEntireCoordinateSequenceFilter
        {
            public bool Done => IsCrossingDateline;

            public bool GeometryChanged => false;

            public bool IsCrossingDateline { get; private set; }

            public void Filter(CoordinateSequence seq)
            {
                // Test if we have a segement crossing the dateline
                double currX = seq.GetX(0);
                bool wasLT = currX < -180d;
                bool wasGT = currX > -180d;
                for (int i = 1; i < seq.Count; i++)
                {
                    double lastX = currX;

                    // Prevent unmotivated sign change on dateline
                    currX = seq.GetX(i);
                    if ((currX == -180d && lastX == 180d) ||
                        (currX == 180d && lastX == -180d))
                        currX = lastX;

                    // Check if we cross the dateline
                    if (currX < -180d && wasGT ||
                        currX > -180d && wasLT)
                    {
                        IsCrossingDateline = true;
                        break;
                    }

                    // Update flag variables
                    wasLT = currX < -180d;
                    wasGT = currX > -180d;
                }
            }
        }

        /// <summary>
        /// A filter class to translate every x-ordinate by a predefined value
        /// </summary>
        /// <remarks>Sets <see cref="GeometryChanged"/> to <c>true</c> if the translation value is <c>!= 0d</c></remarks>
        private class TranslateXFilter : IEntireCoordinateSequenceFilter
        {
            private readonly double _translate;

            public TranslateXFilter(double translate)
            {
                _translate = translate;
            }

            public bool Done => false;

            public bool GeometryChanged => _translate != 0d;

            public void Filter(CoordinateSequence seq)
            {
                // translate whole ge
                for (int i = 0; i < seq.Count; i++)
                    seq.SetX(i, seq.GetX(i) + _translate);
            }
        }
    }
    public class DatelineUtilityTest
    {
        NetTopologySuite.NtsGeometryServices _instance = new NetTopologySuite.NtsGeometryServices(new PrecisionModel(1000));

        [TestCase("SRID=4326;LINESTRING EMPTY", "LINESTRING EMPTY")]
        [TestCase("SRID=4326;LINESTRING (-179 0, -181 0)", "MULTILINESTRING ((180 0, 179 0), (-179 0, -180 0))")]
        [TestCase("SRID=4326;LINESTRING (-541 -20, 0 0, 541 20)",
            "MULTILINESTRING ((179 -20, 180 -19.963), (-180 -19.963, 180 -6.654), (-180 -6.654, 0 0, 180 6.654), (-180 6.654, 180 19.963), (-180 19.963, -179 20))")]
        [TestCase("SRID=4326;LINESTRING (-179 0, -180 0, -180 2, -179 2)", "LINESTRING(-179 0, -180 0, -180 2, -179 2)")]
        [TestCase("SRID=4326;LINESTRING (-179 0, 181 0)", "LINESTRING(-180 0, 180 0)")]
        [TestCase("SRID=4326;LINESTRING (-185 20, -175 15, -185 10, -175 5, -180 2.5, -180 -2.5, -175 -5, -185 -10, -175 -15, -185 -20)",
            "MULTILINESTRING ((175 20, 180 17.5), (180 12.5, 175 10, 180 7.5), (180 2.5, 180 -2.5), (180 -7.5, 175 -10, 180 -12.5), " +
                             "(180 -17.5, 175 -20), (-180 17.5, -175 15, -180 12.5), (-180 7.5, -175 5, -180 2.5), (-180 -2.5, -175 -5, -180 -7.5), " +
                             "(-180 -12.5, -175 -15, -180 -17.5))")]
        [TestCase("SRID=4326;POLYGON ((-179 1, -181 1, -181 -1, -179 -1, -179 1))",
            "MULTIPOLYGON (((-179 1, -179 -1, -180 -1, -180 1, -179 1)), ((179 -1, 179 1, 180 1, 180 -1, 179 -1)))")]
        [TestCase("SRID=4326;POLYGON ((-185 0, -185 10, -175 10, -175 0, -185 0),   (-182 3, -178 3, -178 7, -182 7, -182 3)))",
            "MULTIPOLYGON (((-175 10, -175 0, -180 0, -180 3, -178 3, -178 7, -180 7, -180 10, -175 10)), ((175 10, 180 10, 180 7, 178 7, 178 3, 180 3, 180 0, 175 0, 175 10)))")]
        [TestCase("SRID=4326;LINESTRING (-179 0, -181 0)", "MULTILINESTRING ((180 0, 179 0), (-179 0, -180 0))")]
        public void Test(string wktInput, string wktExpected)
        {
            var rdr = new WKTReader(_instance);
            var input = rdr.Read(wktInput);
            var actual = DatelineUtility.UnwrapAtDateline(input);
            var expected = rdr.Read(wktExpected);

            if (actual.IsEmpty && expected.IsEmpty)
                Assert.That(actual.OgcGeometryType, Is.EqualTo(expected.OgcGeometryType));
            else
                Assert.That(actual.Equals(expected), Is.True,
                    $"Expected: {expected}\nActual: {actual}");
        }
    }
}