diff --git a/src/NetTopologySuite.IO.SqlServerBytes/NetTopologySuite.IO.SqlServerBytes.csproj b/src/NetTopologySuite.IO.SqlServerBytes/NetTopologySuite.IO.SqlServerBytes.csproj index 7750838..e3295bc 100644 --- a/src/NetTopologySuite.IO.SqlServerBytes/NetTopologySuite.IO.SqlServerBytes.csproj +++ b/src/NetTopologySuite.IO.SqlServerBytes/NetTopologySuite.IO.SqlServerBytes.csproj @@ -22,7 +22,7 @@ - + @@ -40,4 +40,8 @@ + + + + diff --git a/src/NetTopologySuite.IO.SqlServerBytes/SqlServerBytesWriter.cs b/src/NetTopologySuite.IO.SqlServerBytes/SqlServerBytesWriter.cs index bb20c5c..46ede1a 100644 --- a/src/NetTopologySuite.IO.SqlServerBytes/SqlServerBytesWriter.cs +++ b/src/NetTopologySuite.IO.SqlServerBytes/SqlServerBytesWriter.cs @@ -4,7 +4,8 @@ using System.Linq; using NetTopologySuite.Geometries; using NetTopologySuite.IO.Properties; - +using NetTopologySuite.IO.Utility; +using NetTopologySuite.Operation.Valid; using Figure = NetTopologySuite.IO.Serialization.Figure; using FigureAttribute = NetTopologySuite.IO.Serialization.FigureAttribute; using Geography = NetTopologySuite.IO.Serialization.Geography; @@ -131,9 +132,13 @@ private Geography ToGeography(Geometry geometry) var geography = new Geography { SRID = Math.Max(0, geometry.SRID), - IsValid = geometry.IsValid + IsValid = CheckValid(geometry), + //IsLargerThanAHemisphere = }; + if (IsGeography) + geography.IsLargerThanAHemisphere = CheckLargerThanAHemisphere(geometry); + while (geometries.Count > 0) { var (currentGeometry, parentOffset) = geometries.Dequeue(); @@ -254,5 +259,23 @@ private OpenGisType ToOpenGisType(OgcGeometryType type) return (OpenGisType)type; } + + private static bool CheckValid(Geometry geometry) + { + // ToDo: Consider deriving special IsValidSqlServerOp class. + // #see issue gh issue #12 + var isValidOp = new IsValidOp(geometry) { + IsSelfTouchingRingFormingHoleValid = true + }; + + return isValidOp.IsValid; + } + + private static bool CheckLargerThanAHemisphere(Geometry geometry) + { + var ge = GeographyEnvelope.Compute(geometry); + return ge.Angle > 90d; + } + } } diff --git a/src/NetTopologySuite.IO.SqlServerBytes/Utility/GeocentricTransform.cs b/src/NetTopologySuite.IO.SqlServerBytes/Utility/GeocentricTransform.cs new file mode 100644 index 0000000..9764f88 --- /dev/null +++ b/src/NetTopologySuite.IO.SqlServerBytes/Utility/GeocentricTransform.cs @@ -0,0 +1,221 @@ +using System; +using NetTopologySuite.Algorithm; + +namespace NetTopologySuite.IO.Utility +{ + /// + /// Wenzel, H.-G.(1985): Hochauflösende Kugelfunktionsmodelle für + /// das Gravitationspotential der Erde. Wiss. Arb. Univ. Hannover + /// Nr. 137, p. 130-131. + /// + /// Programmed by + /// + /// GGA- Leibniz-Institue of Applied Geophysics + /// Stilleweg 2 + /// D-30655 Hannover + /// Federal Republic of Germany + /// + /// Internet: www.gga-hannover.de + /// + /// Hannover, March 1999, April 2004. + /// + ///see also: comments in statements + /// + /// + /// Mathematically exact and because of symmetry of rotation-ellipsoid, + /// each point(X, Y, Z) has at least two solutions(Latitude1, Longitude1, Height1) and + /// (Latitude2, Longitude2, Height2). Is point = (0., 0., Z)(P = 0.), so you get even + /// four solutions, every two symmetrical to the semi-minor axis. + /// Here Height1 and Height2 have at least a difference in order of + /// radius of curvature (e.g. (0, 0, b) => (90., 0., 0.) or(-90., 0., -2b); + /// (a+100.)*(sqrt(2.)/2., sqrt(2.)/2., 0.) => (0., 45., 100.) or + /// (0., 225., -(2a+100.))). + /// The algorithm always computes(Latitude, Longitude) with smallest |Height|. + /// For normal computations, that means |Height|<10000.m, algorithm normally + /// converges after to 2-3 steps!!! + /// But if |Height| has the amount of length of ellipsoid's axis + /// (e.g. -6300000.m), algorithm needs about 15 steps. + /// + public class GeocentricTransform + { + public static GeocentricTransform CreateFor(int srid) + { + switch (srid) + { + default: + case 4326: + return new GeocentricTransform(6378137, 6378137); + //return new GeocentricTransform(6378137, 6356752); + } + } + + /* local defintions and variables */ + /* end-criterium of loop, accuracy of sin(Latitude) */ + private const double GENAU = 1E-12; + private const double GENAU2 = GENAU * GENAU; + + private const int MAXITER = 30; + + // private const double COS_67P5 = 0.38268343236508977; /* cosine of 67.5 degrees */ + // private const double AD_C = 1.0026000; /* Toms region 1 constant */ + private const double PI = 3.14159265358979323e0; + private const double PI_OVER_2 = (PI / 2.0e0); + + private readonly double _a; + private readonly double _a2; + private readonly double _b; + private readonly double _b2; + private readonly double _e2; + + /// + /// Creates a new instance of GeocentricGeodetic + /// + public GeocentricTransform(double equatorialRadius, double polarRadius) + { + _a = equatorialRadius; + _b = polarRadius; + _a2 = _a * _a; + _b2 = _b * _b; + _e2 = (_a2 - _b2) / _a2; + //_ep2 = (_a2 - _b2)/_b2; + } + + /// + /// Converts lon, lat, height to x, y, z where lon and lat are in radians and everything else is meters + /// + /// + /// + /// + public (double x, double y, double z) GeodeticToGeocentric(double lon, double lat, double height) + { + lon = AngleUtility.ToRadians(lon); + lat = AngleUtility.ToRadians(lat); + + /* + * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates + * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z), + * according to the current ellipsoid parameters. + * + * Latitude : Geodetic latitude in radians (input) + * Longitude : Geodetic longitude in radians (input) + * Height : Geodetic height, in meters (input) + * X : Calculated Geocentric X coordinate, in meters (output) + * Y : Calculated Geocentric Y coordinate, in meters (output) + * Z : Calculated Geocentric Z coordinate, in meters (output) + * + */ + + /* + ** Don't blow up if Latitude is just a little out of the value + ** range as it may just be a rounding issue. Also removed longitude + ** test, it should be wrapped by cos() and sin(). NFW for PROJ.4, Sep/2001. + */ + if (lat < -PI_OVER_2 && lat > -1.001 * PI_OVER_2) + lat = -PI_OVER_2; + else if (lat > PI_OVER_2 && lat < 1.001 * PI_OVER_2) + lat = PI_OVER_2; + else if ((lat < -PI_OVER_2) || (lat > PI_OVER_2)) + { + /* lat out of range */ + return (double.NaN, double.NaN, double.NaN); + } + + if (lon > PI) + lon -= (2 * PI); + double sinLat = Math.Sin(lat); + double cosLat = Math.Cos(lat); + double sin2Lat = sinLat * sinLat; /* Square of sin(Latitude) */ + double rn = _a / (Math.Sqrt(1.0e0 - _e2 * sin2Lat)); /* Earth radius at location */ + double x = (rn + height) * cosLat * Math.Cos(lon); + double y = (rn + height) * cosLat * Math.Sin(lon); + double z = ((rn * (1 - _e2)) + height) * sinLat; + + return (x, y, z); + } + + /// + /// Converts x, y, z to lon, lat, height + /// + /// + /// + /// + public (double lon, double lat, double height) GeocentricToGeodetic(double x, double y, double z) + { + double lon; + double lat; + double height; + + double cosPhi; /* cos of searched geodetic latitude */ + double sinPhi; /* sin of searched geodetic latitude */ + double sinDiffPhi; /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */ + //bool At_Pole = false; /* indicates location is in polar region */ + + double p = Math.Sqrt(x * x + y * y); + double rr = Math.Sqrt(x * x + y * y + z * z); /* distance between center and location */ + + /* special cases for latitude and longitude */ + if (p / _a < GENAU) + { + /* special case, if P=0. (X=0., Y=0.) */ + //At_Pole = true; + lon = 0.0; + + /* if (X, Y, Z)=(0., 0., 0.) then Height becomes semi-minor axis + * of ellipsoid (=center of mass), Latitude becomes PI/2 */ + if (rr / _a < GENAU) + { + lat = PI_OVER_2; + height = -_b; + return (lon, AngleUtility.ToDegrees(lat), height); + } + } + else + { + /* ellipsoidal (geodetic) longitude + * interval: -PI < Longitude <= +PI */ + lon = Math.Atan2(y, x); + } + + /* -------------------------------------------------------------- + * Following iterative algorithm was developped by + * "Institut für Erdmessung", University of Hannover, July 1988. + * Internet: www.ife.uni-hannover.de + * Iterative computation of CPHI, SPHI and Height. + * Iteration of CPHI and SPHI to 10**-12 radian resp. + * 2*10**-7 arcsec. + * -------------------------------------------------------------- + */ + double ct = z / rr; // sin of geocentric latitude // looks like these two should be flipped (TD). + double st = p / rr; // cos of geocentric latitude + double rx = 1.0 / Math.Sqrt(1.0 - _e2 * (2.0 - _e2) * st * st); + double cosPhi0 = st * (1.0 - _e2) * rx; /* cos of start or old geodetic latitude in iterations */ + double sinPhi0 = ct * rx; /* sin of start or old geodetic latitude in iterations */ + int iter = 0; /* # of continous iteration, max. 30 is always enough (s.a.) */ + + /* loop to find sin(Latitude) resp. Latitude + * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */ + do + { + iter++; + double earthRadius = _a / Math.Sqrt(1.0 - _e2 * sinPhi0 * sinPhi0); + + /* ellipsoidal (geodetic) height */ + height = p * cosPhi0 + z * sinPhi0 - earthRadius * (1.0 - _e2 * sinPhi0 * sinPhi0); + + double rk = _e2 * earthRadius / (earthRadius + height); + rx = 1.0 / Math.Sqrt(1.0 - rk * (2.0 - rk) * st * st); + cosPhi = st * (1.0 - rk) * rx; + sinPhi = ct * rx; + sinDiffPhi = sinPhi * cosPhi0 - cosPhi * sinPhi0; + cosPhi0 = cosPhi; + sinPhi0 = sinPhi; + } while (sinDiffPhi * sinDiffPhi > GENAU2 && iter < MAXITER); + + /* ellipsoidal (geodetic) latitude */ + lat = Math.Atan(sinPhi / Math.Abs(cosPhi)); + + return (AngleUtility.ToDegrees(lon), AngleUtility.ToDegrees(lat), height); + } + } +} + diff --git a/src/NetTopologySuite.IO.SqlServerBytes/Utility/GeographyEnvelope.cs b/src/NetTopologySuite.IO.SqlServerBytes/Utility/GeographyEnvelope.cs new file mode 100644 index 0000000..2a1ade2 --- /dev/null +++ b/src/NetTopologySuite.IO.SqlServerBytes/Utility/GeographyEnvelope.cs @@ -0,0 +1,123 @@ +using System; +using System.Collections.Generic; +using NetTopologySuite.Algorithm; +using NetTopologySuite.Geometries; +using NetTopologySuite.Mathematics; + +namespace NetTopologySuite.IO.Utility +{ + /// + /// A geography envelope is made up of a and an + /// + internal struct GeographyEnvelope + { + /// + /// Computes the envelope of a geography + /// + /// The geometry + /// The envelope of a geography + public static GeographyEnvelope Compute(Geometry geometry) + { + if (geometry.IsEmpty) + return new GeographyEnvelope(); + + var flt = new GeographyEnvelopeFilter(GeocentricTransform.CreateFor(geometry.SRID)); + geometry.Apply(flt); + flt.Compute(out var coordinate, out double angle); + + return new GeographyEnvelope { Angle = angle, Center = coordinate}; + } + + /// + /// Gets the opening angle between and a vector to the + /// envelope circle. + /// + public double Angle { get; private set; } + + /// + /// Gets a value indicating the center of the envelope + /// + public Coordinate Center { get; private set; } + + private class GeographyEnvelopeFilter : IEntireCoordinateSequenceFilter + { + private readonly GeocentricTransform _transform; + private readonly List _vectors = new List(); + + /// + /// Creates an instance of this class + /// + /// + public GeographyEnvelopeFilter(GeocentricTransform transform) + { + _transform = transform; + } + + void IEntireCoordinateSequenceFilter.Filter(CoordinateSequence seq) + { + int count = seq.Count; + if (CoordinateSequences.IsRing(seq)) + count--; + + for (int i = 0; i < count; i++) + { + var xyz = _transform.GeodeticToGeocentric(seq.GetX(i), seq.GetY(i), 0d); + _vectors.Add(new Vector3D(xyz.x, xyz.y, xyz.z)); + } + } + + bool IEntireCoordinateSequenceFilter.Done => false; + + bool IEntireCoordinateSequenceFilter.GeometryChanged => false; + + public void Compute(out CoordinateZ center, out double angle) + { + // Build average center vector + var vCenterX = new DD(_vectors[0].X); + var vCenterY = new DD(_vectors[0].Y); + var vCenterZ = new DD(_vectors[0].Z); + var divisor = new DD(_vectors.Count); + for (int i = 1; i < _vectors.Count; i++) + { + vCenterX += new DD(_vectors[i].X); + vCenterY += new DD(_vectors[i].Y); + vCenterZ += new DD(_vectors[i].Z); + } + var vCenter = new Vector3D( + (vCenterX / divisor).ToDoubleValue(), + (vCenterY / divisor).ToDoubleValue(), + (vCenterZ / divisor).ToDoubleValue()); + + // look for max angle between center and vectors + angle = -double.MaxValue; + double vCenterLength = vCenter.Length(); + for (int i = 0; i < _vectors.Count; i++) + { + double dot = vCenter.Dot(_vectors[i]); + double testAngle; + if (dot == 0d) + { + testAngle = 90d; + } + else + { + double cosPhi = dot / (vCenterLength * _vectors[i].Length()); + testAngle = AngleUtility.ToDegrees(Math.Acos(cosPhi)); + // Is angle obtuse + if (dot < 0) + testAngle += 180d; + } + + // Update largest angle + if (testAngle > angle) + angle = testAngle; + } + + // Transform to geodetic + var centerLL = _transform.GeocentricToGeodetic(vCenter.X, vCenter.Y, vCenter.Z); + center = new CoordinateZ(centerLL.lon, centerLL.lat, centerLL.height); + } + + } + } +} diff --git a/test/NetTopologySuite.IO.SqlServerBytes.Test/NetTopologySuite.IO.SqlServerBytes.Test.csproj b/test/NetTopologySuite.IO.SqlServerBytes.Test/NetTopologySuite.IO.SqlServerBytes.Test.csproj index 437d86d..0443a5f 100644 --- a/test/NetTopologySuite.IO.SqlServerBytes.Test/NetTopologySuite.IO.SqlServerBytes.Test.csproj +++ b/test/NetTopologySuite.IO.SqlServerBytes.Test/NetTopologySuite.IO.SqlServerBytes.Test.csproj @@ -9,4 +9,8 @@ + + + + diff --git a/test/NetTopologySuite.IO.SqlServerBytes.Test/Utility/GeographyEnvelopeTest.cs b/test/NetTopologySuite.IO.SqlServerBytes.Test/Utility/GeographyEnvelopeTest.cs new file mode 100644 index 0000000..1cdd2a5 --- /dev/null +++ b/test/NetTopologySuite.IO.SqlServerBytes.Test/Utility/GeographyEnvelopeTest.cs @@ -0,0 +1,71 @@ +using System; +using System.Collections.Generic; +using NetTopologySuite.Geometries; +using NetTopologySuite.Geometries.Implementation; +using Xunit; + +namespace NetTopologySuite.IO.Utility +{ + public class GeographyEnvelopeTest + { + private NtsGeometryServices _ntsServices = + new NtsGeometryServices(CoordinateArraySequenceFactory.Instance, new PrecisionModel(), 4326); + + private readonly DoubleToleranceComparer _toleranceComparer = new DoubleToleranceComparer(2e-6); + + /* + * SQL statements for inline data: + * + + WITH T AS (SELECT geography::STGeomFromText('POINT (0 90)', 4326).STBuffer(500) AS g) + SELECT g.STIsValid(), g.STAsText(), g.EnvelopeCenter().STAsText(), g.EnvelopeAngle() FROM T; + + WITH T AS (SELECT geography::STGeomFromText('LINESTRING (7 52, 9 49)', 4326) AS g) + SELECT g.STIsValid(), g.STAsText(), g.EnvelopeCenter().STAsText(), g.EnvelopeAngle() FROM T; + + */ + + [Theory] + [InlineData("POINT (1 2)", 1, 2, 0)] + [InlineData("LINESTRING (7 52, 9 49)", 8.03176928564099, 50.504278988995246, 1.6291318048158)] + [InlineData("POLYGON ((0 89.995523482987139, 2.8124999999977969 89.995523482987139, 5.6249999999956 89.995523482987139, 8.43749999999341 89.995523482987139, 11.249999999991239 89.995523482987139, 14.062499999989088 89.995523482987139, 16.874999999986965 89.995523482987139, 19.687499999984873 89.995523482987139, 22.499999999982819 89.995523482987139, 25.312499999980805 89.995523482987139, 28.124999999978833 89.995523482987139, 30.937499999976911 89.995523482987139, 33.749999999975053 89.995523482987139, 36.562499999973248 89.995523482987139, 39.374999999971507 89.995523482987139, 42.187499999969845 89.995523482987139, 44.999999999968253 89.995523482987139, 47.812499999966732 89.995523482987139, 50.6249999999653 89.995523482987139, 53.43749999996394 89.995523482987139, 56.249999999962668 89.995523482987139, 59.062499999961496 89.995523482987139, 61.8749999999604 89.995523482987139, 64.687499999959414 89.995523482987139, 67.499999999958519 89.995523482987139, 70.312499999957723 89.995523482987139, 73.124999999957026 89.995523482987139, 75.937499999956444 89.995523482987139, 78.749999999955975 89.995523482987139, 81.562499999955591 89.995523482987139, 84.3749999999553 89.995523482987139, 87.187499999955151 89.995523482987139, 89.999999999955108 89.995523482987139, 92.812499999955151 89.995523482987139, 95.6249999999553 89.995523482987139, 98.437499999955577 89.995523482987139, 101.24999999995597 89.995523482987139, 104.06249999995644 89.995523482987139, 106.87499999995703 89.995523482987139, 109.68749999995772 89.995523482987139, 112.49999999995852 89.995523482987139, 115.3124999999594 89.995523482987139, 118.12499999996041 89.995523482987139, 120.93749999996147 89.995523482987139, 123.74999999996265 89.995523482987139, 126.56249999996392 89.995523482987139, 129.3749999999653 89.995523482987139, 132.18749999996672 89.995523482987139, 134.99999999996825 89.995523482987139, 137.81249999996987 89.995523482987139, 140.62499999997152 89.995523482987139, 143.43749999997326 89.995523482987139, 146.24999999997505 89.995523482987139, 149.06249999997692 89.995523482987139, 151.87499999997885 89.995523482987139, 154.68749999998079 89.995523482987139, 157.49999999998283 89.995523482987139, 160.31249999998488 89.995523482987139, 163.12499999998698 89.995523482987139, 165.93749999998911 89.995523482987139, 168.74999999999125 89.995523482987139, 171.56249999999341 89.995523482987139, 174.37499999999562 89.995523482987139, 177.18749999999778 89.995523482987139, 180 89.995523482987139, -177.18749999999778 89.995523482987139, -174.37499999999562 89.995523482987139, -171.56249999999341 89.995523482987139, -168.74999999999125 89.995523482987139, -165.93749999998911 89.995523482987139, -163.12499999998698 89.995523482987139, -160.31249999998488 89.995523482987139, -157.49999999998283 89.995523482987139, -154.68749999998079 89.995523482987139, -151.87499999997885 89.995523482987139, -149.06249999997692 89.995523482987139, -146.24999999997505 89.995523482987139, -143.43749999997326 89.995523482987139, -140.62499999997152 89.995523482987139, -137.81249999996987 89.995523482987139, -134.99999999996825 89.995523482987139, -132.18749999996672 89.995523482987139, -129.3749999999653 89.995523482987139, -126.56249999996392 89.995523482987139, -123.74999999996265 89.995523482987139, -120.93749999996147 89.995523482987139, -118.12499999996041 89.995523482987139, -115.3124999999594 89.995523482987139, -112.49999999995852 89.995523482987139, -109.68749999995772 89.995523482987139, -106.87499999995703 89.995523482987139, -104.06249999995644 89.995523482987139, -101.24999999995597 89.995523482987139, -98.437499999955577 89.995523482987139, -95.6249999999553 89.995523482987139, -92.812499999955151 89.995523482987139, -89.999999999955108 89.995523482987139, -87.187499999955151 89.995523482987139, -84.3749999999553 89.995523482987139, -81.562499999955591 89.995523482987139, -78.749999999955975 89.995523482987139, -75.937499999956444 89.995523482987139, -73.124999999957026 89.995523482987139, -70.312499999957723 89.995523482987139, -67.499999999958519 89.995523482987139, -64.687499999959414 89.995523482987139, -61.8749999999604 89.995523482987139, -59.062499999961496 89.995523482987139, -56.249999999962668 89.995523482987139, -53.43749999996394 89.995523482987139, -50.6249999999653 89.995523482987139, -47.812499999966732 89.995523482987139, -44.999999999968253 89.995523482987139, -42.187499999969845 89.995523482987139, -39.374999999971507 89.995523482987139, -36.562499999973248 89.995523482987139, -33.749999999975053 89.995523482987139, -30.937499999976911 89.995523482987139, -28.124999999978833 89.995523482987139, -25.312499999980805 89.995523482987139, -22.499999999982819 89.995523482987139, -19.687499999984873 89.995523482987139, -16.874999999986965 89.995523482987139, -14.062499999989088 89.995523482987139, -11.249999999991239 89.995523482987139, -8.43749999999341 89.995523482987139, -5.6249999999956 89.995523482987139, -2.8124999999977969 89.995523482987139, 0 89.995523482987139))", + 0.0011397786255458117, 90d, 0.00447651701285935, + Skip = "A buffer around the north pole should result in lon = 0d, shouldn't it?")] + [InlineData("POLYGON ((0 89.995523482987139, 2.8124999999977969 89.995523482987139, 5.6249999999956 89.995523482987139, 8.43749999999341 89.995523482987139, 11.249999999991239 89.995523482987139, 14.062499999989088 89.995523482987139, 16.874999999986965 89.995523482987139, 19.687499999984873 89.995523482987139, 22.499999999982819 89.995523482987139, 25.312499999980805 89.995523482987139, 28.124999999978833 89.995523482987139, 30.937499999976911 89.995523482987139, 33.749999999975053 89.995523482987139, 36.562499999973248 89.995523482987139, 39.374999999971507 89.995523482987139, 42.187499999969845 89.995523482987139, 44.999999999968253 89.995523482987139, 47.812499999966732 89.995523482987139, 50.6249999999653 89.995523482987139, 53.43749999996394 89.995523482987139, 56.249999999962668 89.995523482987139, 59.062499999961496 89.995523482987139, 61.8749999999604 89.995523482987139, 64.687499999959414 89.995523482987139, 67.499999999958519 89.995523482987139, 70.312499999957723 89.995523482987139, 73.124999999957026 89.995523482987139, 75.937499999956444 89.995523482987139, 78.749999999955975 89.995523482987139, 81.562499999955591 89.995523482987139, 84.3749999999553 89.995523482987139, 87.187499999955151 89.995523482987139, 89.999999999955108 89.995523482987139, 92.812499999955151 89.995523482987139, 95.6249999999553 89.995523482987139, 98.437499999955577 89.995523482987139, 101.24999999995597 89.995523482987139, 104.06249999995644 89.995523482987139, 106.87499999995703 89.995523482987139, 109.68749999995772 89.995523482987139, 112.49999999995852 89.995523482987139, 115.3124999999594 89.995523482987139, 118.12499999996041 89.995523482987139, 120.93749999996147 89.995523482987139, 123.74999999996265 89.995523482987139, 126.56249999996392 89.995523482987139, 129.3749999999653 89.995523482987139, 132.18749999996672 89.995523482987139, 134.99999999996825 89.995523482987139, 137.81249999996987 89.995523482987139, 140.62499999997152 89.995523482987139, 143.43749999997326 89.995523482987139, 146.24999999997505 89.995523482987139, 149.06249999997692 89.995523482987139, 151.87499999997885 89.995523482987139, 154.68749999998079 89.995523482987139, 157.49999999998283 89.995523482987139, 160.31249999998488 89.995523482987139, 163.12499999998698 89.995523482987139, 165.93749999998911 89.995523482987139, 168.74999999999125 89.995523482987139, 171.56249999999341 89.995523482987139, 174.37499999999562 89.995523482987139, 177.18749999999778 89.995523482987139, 180 89.995523482987139, -177.18749999999778 89.995523482987139, -174.37499999999562 89.995523482987139, -171.56249999999341 89.995523482987139, -168.74999999999125 89.995523482987139, -165.93749999998911 89.995523482987139, -163.12499999998698 89.995523482987139, -160.31249999998488 89.995523482987139, -157.49999999998283 89.995523482987139, -154.68749999998079 89.995523482987139, -151.87499999997885 89.995523482987139, -149.06249999997692 89.995523482987139, -146.24999999997505 89.995523482987139, -143.43749999997326 89.995523482987139, -140.62499999997152 89.995523482987139, -137.81249999996987 89.995523482987139, -134.99999999996825 89.995523482987139, -132.18749999996672 89.995523482987139, -129.3749999999653 89.995523482987139, -126.56249999996392 89.995523482987139, -123.74999999996265 89.995523482987139, -120.93749999996147 89.995523482987139, -118.12499999996041 89.995523482987139, -115.3124999999594 89.995523482987139, -112.49999999995852 89.995523482987139, -109.68749999995772 89.995523482987139, -106.87499999995703 89.995523482987139, -104.06249999995644 89.995523482987139, -101.24999999995597 89.995523482987139, -98.437499999955577 89.995523482987139, -95.6249999999553 89.995523482987139, -92.812499999955151 89.995523482987139, -89.999999999955108 89.995523482987139, -87.187499999955151 89.995523482987139, -84.3749999999553 89.995523482987139, -81.562499999955591 89.995523482987139, -78.749999999955975 89.995523482987139, -75.937499999956444 89.995523482987139, -73.124999999957026 89.995523482987139, -70.312499999957723 89.995523482987139, -67.499999999958519 89.995523482987139, -64.687499999959414 89.995523482987139, -61.8749999999604 89.995523482987139, -59.062499999961496 89.995523482987139, -56.249999999962668 89.995523482987139, -53.43749999996394 89.995523482987139, -50.6249999999653 89.995523482987139, -47.812499999966732 89.995523482987139, -44.999999999968253 89.995523482987139, -42.187499999969845 89.995523482987139, -39.374999999971507 89.995523482987139, -36.562499999973248 89.995523482987139, -33.749999999975053 89.995523482987139, -30.937499999976911 89.995523482987139, -28.124999999978833 89.995523482987139, -25.312499999980805 89.995523482987139, -22.499999999982819 89.995523482987139, -19.687499999984873 89.995523482987139, -16.874999999986965 89.995523482987139, -14.062499999989088 89.995523482987139, -11.249999999991239 89.995523482987139, -8.43749999999341 89.995523482987139, -5.6249999999956 89.995523482987139, -2.8124999999977969 89.995523482987139, 0 89.995523482987139))", + 0d, 90d, 0.00447651701285935)] + public void TestCenterAndAngle(string wkt, double lon, double lat, double angle) + { + // Arrange + var rdr = new WKTReader { DefaultSRID = 4326 }; + var geometry = rdr.Read(wkt); + + // Act + var ge = GeographyEnvelope.Compute(geometry); + + // Assert + Assert.Equal(lon, ge.Center.X, _toleranceComparer); + Assert.Equal(lat, ge.Center.Y, _toleranceComparer); + Assert.Equal(angle, ge.Angle, _toleranceComparer); + } + + private class DoubleToleranceComparer : EqualityComparer + { + private readonly double _tolerance; + + public DoubleToleranceComparer(double tolerance) + { + _tolerance = tolerance; + } + + public override bool Equals(double x, double y) + { + return Math.Abs(x - y) <= _tolerance; + } + + public override int GetHashCode(double obj) + { + return obj.GetHashCode(); + } + } + } +}