C# IPosition.GetCoordinate方法代码示例

		public static SpatialCriteria WithinRadiusOf(this SpatialCriteriaFactory @this,
													double radius,
													IPosition position)
		{
			var coordinate = position.GetCoordinate();
			return @this.WithinRadiusOf(radius, coordinate.Longitude, coordinate.Latitude);
		}

 private double[] WriteCoordinate(IPosition position)
 {
     var coordinate = position.GetCoordinate();
     var pointZM = coordinate as CoordinateZM;
     if (pointZM != null)
         return new[] { pointZM.Longitude, pointZM.Latitude, pointZM.Elevation, pointZM.Measure };
     var pointZ = coordinate as CoordinateZ;
     if (pointZ != null)
         return new[] { pointZ.Longitude, pointZ.Latitude, pointZ.Elevation };
     //CoordinateM is not supported by GeoJSON
     return new[] { coordinate.Longitude, coordinate.Latitude };
 }

        public GeodeticLine CalculateLoxodromicLine(IPosition position1, IPosition position2)
        {
            var point1 = position1.GetCoordinate();
            var point2 = position2.GetCoordinate();
            var lat1 = point1.Latitude;
            var lon1 = point1.Longitude;
            var lat2 = point2.Latitude;
            var lon2 = point2.Longitude;

            if (Math.Abs(lat1 - lat2) < double.Epsilon && Math.Abs(lon1 - lon2) < double.Epsilon)
                return null;

            double distance;
            var latDeltaRad = (lat2 - lat1).ToRadians();
            var meridionalDistance = CalculateMeridionalDistance(lat2) - CalculateMeridionalDistance(lat1);
            var course = LoxodromicLineCourse(lat1, lon1, lat2, lon2);

            if (Math.Abs(latDeltaRad) < 0.0008)
            {
                // near parallel sailing

                var lonDelta = lon2 - lon1;
                if (lonDelta > 180)
                    lonDelta = lonDelta - 360;
                if (lonDelta < -180)
                    lonDelta = lonDelta + 360;
                var lonDeltaRad = lonDelta.ToRadians();

                var midLatRad = (0.5 * (lat1 + lat2)).ToRadians();
                // expand merid_dist/dmp about lat_mid_rad to order e2*dlat_rad^2
                var e2 = Math.Pow(Spheroid.Eccentricity, 2);
                var ratio = Math.Cos(midLatRad) /
                        Math.Sqrt(1 - e2 * Math.Pow(Math.Sin(midLatRad), 2)) *
                    (1.0 + (e2 * Math.Cos(2 * midLatRad) / 8 -
                        (1 + 2 * Math.Pow(Math.Tan(midLatRad), 2)) / 24 -
                        e2 / 12) * latDeltaRad * latDeltaRad);

                distance = Math.Sqrt(Math.Pow(meridionalDistance, 2) + Math.Pow(Spheroid.EquatorialAxis * ratio * lonDeltaRad, 2));
            }
            else
            {
                distance = Math.Abs(meridionalDistance / Math.Cos(course.ToRadians()));
            }
            return new GeodeticLine(new Coordinate(lat1, lon1), new Coordinate(lat2, lon2), distance, course, course > 180 ? course - 180 : course + 180);
        }

        private double[] CalculateOrthodromicLineInternal(IPosition position1, IPosition position2)
        {
            var point1 = position1.GetCoordinate();
            var point2 = position2.GetCoordinate();

            if (Math.Abs(point1.Latitude - point2.Latitude) < double.Epsilon && Math.Abs(point1.Longitude - point2.Longitude) < double.Epsilon)
                return null;

            var lon1 = point1.Longitude.ToRadians();
            var lat1 = point1.Latitude.ToRadians();
            var lon2 = point2.Longitude.ToRadians();
            var lat2 = point2.Latitude.ToRadians();
            /*
             * Solution of the geodetic inverse problem after T.Vincenty.
             * Modified Rainsford's method with Helmert's elliptical terms.
             * Effective in any azimuth and at any distance short of antipodal.
             *
             * Latitudes and longitudes in radians positive North and East.
             * Forward azimuths at both points returned in radians from North.
             *
             * Programmed for CDC-6600 by LCDR L.Pfeifer NGS ROCKVILLE MD 18FEB75
             * Modified for IBM SYSTEM 360 by John G.Gergen NGS ROCKVILLE MD 7507
             * Ported from Fortran to Java by Martin Desruisseaux.
             *
             * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/inverse.for
             *         subroutine INVER1
             */
            const int maxIterations = 100;
            const double eps = 0.5E-13;
            double R = 1 - Spheroid.Flattening;

            double tu1 = R * Math.Sin(lat1) / Math.Cos(lat1);
            double tu2 = R * Math.Sin(lat2) / Math.Cos(lat2);
            double cu1 = 1 / Math.Sqrt(tu1 * tu1 + 1);
            double cu2 = 1 / Math.Sqrt(tu2 * tu2 + 1);
            double su1 = cu1 * tu1;
            double s = cu1 * cu2;
            double baz = s * tu2;
            double faz = baz * tu1;
            double x = lon2 - lon1;
            for (int i = 0; i < maxIterations; i++)
            {
                double sx = Math.Sin(x);
                double cx = Math.Cos(x);
                tu1 = cu2 * sx;
                tu2 = baz - su1 * cu2 * cx;
                double sy = Math.Sqrt(Math.Pow(tu1, 2) + Math.Pow(tu2, 2));
                double cy = s * cx + faz;
                double y = Math.Atan2(sy, cy);
                double SA = s * sx / sy;
                double c2a = 1 - SA * SA;
                double cz = faz + faz;
                if (c2a > 0)
                {
                    cz = -cz / c2a + cy;
                }
                double e = cz * cz * 2 - 1;
                double c = ((-3 * c2a + 4) * Spheroid.Flattening + 4) * c2a * Spheroid.Flattening / 16;
                double d = x;
                x = ((e * cy * c + cz) * sy * c + y) * SA;
                x = (1 - c) * x * Spheroid.Flattening + lon2 - lon1;

                if (Math.Abs(d - x) <= eps)
                {
                    x = Math.Sqrt((1 / (R * R) - 1) * c2a + 1) + 1;
                    x = (x - 2) / x;
                    c = 1 - x;
                    c = (x * x / 4 + 1) / c;
                    d = (0.375 * x * x - 1) * x;
                    x = e * cy;
                    s = 1 - 2 * e;
                    s = ((((sy * sy * 4 - 3) * s * cz * d / 6 - x) * d / 4 + cz) * sy * d + y) * c * R * Spheroid.EquatorialAxis;
                    // 'faz' and 'baz' are forward azimuths at both points.
                    faz = Math.Atan2(tu1, tu2);
                    baz = Math.Atan2(cu1 * sx, baz * cx - su1 * cu2) + Math.PI;
                    return new[] { s, faz.ToDegrees(), baz.ToDegrees() };
                }
            }
            // No convergence. It may be because coordinate points
            // are equals or because they are at antipodes.
            const double leps = 1E-10;
            if (Math.Abs(lon1 - lon2) <= leps && Math.Abs(lat1 - lat2) <= leps)
            {
                // Coordinate points are equals
                return null;
            }
            if (Math.Abs(lat1) <= leps && Math.Abs(lat2) <= leps)
            {
                // Points are on the equator.
                return new[] { Math.Abs(lon1 - lon2) * Spheroid.EquatorialAxis, faz.ToDegrees(), baz.ToDegrees() };
            }
            // Other cases: no solution for this algorithm.
            throw new ArithmeticException();
        }

 public GeodeticLine CalculateOrthodromicLine(IPosition position1, IPosition position2)
 {
     var result = CalculateOrthodromicLineInternal(position1, position2);
     if (result == null)
         return null;
     return new GeodeticLine(position1.GetCoordinate(), position2.GetCoordinate(), result[0], result[1], result[2]);
 }

        public bool TryCalculate(IPosition position, DateTime utcDate, out GeomagnetismResult result)
        {
            var coordinate = position.GetCoordinate();
            var coordinateZ = coordinate as CoordinateZ ?? new CoordinateZ(coordinate.Latitude, coordinate.Longitude, 0);

            double lat = coordinateZ.Latitude.ToRadians(),
                lon = coordinateZ.Longitude.ToRadians(),
                ele = coordinateZ.Elevation / 1000,
                dat = JulianDate.JD(utcDate);

            var model = _models.SingleOrDefault(mod => mod.ValidFrom <= utcDate && mod.ValidTo > utcDate);

            if (model == null)
            {
                result = default(GeomagnetismResult);
                return false;
            }

            var bound = 1 + model.MainCoefficientsG.GetUpperBound(0);

            var sinLat = Math.Sin(lat);
            var cosLat = Math.Cos(lat);
            var a = _spheroid.EquatorialAxis / 1000;
            var f = _spheroid.Flattening;
            var b = a * (1.0 - f);

            var sinLat2 = sinLat * sinLat;
            var cosLat2 = cosLat * cosLat;
            var a2 = a * a;
            var a4 = a2 * a2;
            var b2 = b * b;
            var b4 = b2 * b2;

            var sr = Math.Sqrt(a2 * cosLat2 + b2 * sinLat2);
            var theta = Math.Atan2(cosLat * (ele * sr + a2), sinLat * (ele * sr + b2));
            var r = ele * ele + 2.0 * ele * sr + (a4 - (a4 - b4) * sinLat2) / (a2 - (a2 - b2) * sinLat2);

            r = Math.Sqrt(r);

            var c = Math.Cos(theta);
            var s = Math.Sin(theta);

            double invS;
            if (Math.Abs(s - 0) < double.Epsilon)
                invS = 1.0 / (s + 1E-08);
            else
                invS = 1.0 / (s + 0.0);

            var p = new double[bound, bound];
            var dp = new double[bound, bound];

            p[0, 0] = 1;
            p[1, 1] = s;
            dp[0, 0] = 0;
            dp[1, 1] = c;
            p[1, 0] = c;
            dp[1, 0] = -s;

            for (var i = 2; i < bound; i++)
            {
                var root = Math.Sqrt((2.0 * i - 1) / (2.0 * i));
                p[i, i] = p[i - 1, i - 1] * s * root;
                dp[i, i] = (dp[i - 1, i - 1] * s + p[i - 1, i - 1] * c) * root;
            }

            for (var i = 0; i < bound; i++)
            {
                double i2 = i*i;
                for (var j = Math.Max(i + 1, 2); j < bound; j++)
                {
                    var root1 = Math.Sqrt((j - 1) * (j - 1) - i2);
                    var root2 = 1.0 / Math.Sqrt(j * j - i2);
                    p[j, i] = (p[j - 1, i] * c * (2.0 * j - 1) - p[j - 2, i] * root1) * root2;
                    dp[j, i] = ((dp[j - 1, i] * c - p[j - 1, i] * s) * (2.0 * j - 1) - dp[j - 2, i] * root1) * root2;
                }
            }

            double[,] g = new double[bound, bound], h = new double[bound, bound];
            double bRadial = 0.0, bTheta = 0.0, bPhi = 0.0;
            var fn0 = _spheroid.MeanRadius / 1000 / r;
            var fn = fn0 * fn0;

            double[] sm = new double[bound], cm = new double[bound];
            sm[0] = Math.Sin(0);
            cm[0] = Math.Cos(0);

            var yearfrac = (dat - JulianDate.JD(model.ValidFrom)) / 365.25;
            for (var i = 1; i < bound; i++)
            {
                sm[i] = Math.Sin(i * lon);
                cm[i] = Math.Cos(i * lon);

                for (var j = 0; j < bound; j++)
                {
                    g[i, j] = model.MainCoefficientsG[i, j] + yearfrac * model.SecularCoefficientsG[i, j];
                    h[i, j] = model.MainCoefficientsH[i, j] + yearfrac * model.SecularCoefficientsH[i, j];
                }

                double c1 = 0, c2 = 0, c3 = 0;
                for (var j = 0; j <= i; j++)
//.........这里部分代码省略.........