diff --git a/.github/workflows/full-ci.yml b/.github/workflows/full-ci.yml index 134f301..529c6c6 100644 --- a/.github/workflows/full-ci.yml +++ b/.github/workflows/full-ci.yml @@ -15,10 +15,10 @@ jobs: - name: Get source uses: actions/checkout@v2 - - name: Setup .NET Core 3.1 - uses: actions/setup-dotnet@v1 + - name: Setup .NET Core 6 + uses: actions/setup-dotnet@v4 with: - dotnet-version: 3.1.406 + dotnet-version: '6.x' - name: Build run: dotnet build -c Release -v minimal -p:WarningLevel=3 @@ -44,10 +44,10 @@ jobs: if: github.event_name == 'push' && (github.ref == 'refs/heads/develop' || github.ref == 'refs/heads/master') steps: - - name: Setup .NET Core 3.1 - uses: actions/setup-dotnet@v1 + - name: Setup .NET Core 6 + uses: actions/setup-dotnet@v4 with: - dotnet-version: 3.1.406 + dotnet-version: '6.x' - name: Download Package Files uses: actions/download-artifact@v2 @@ -69,10 +69,10 @@ jobs: if: github.event_name == 'push' && github.ref == 'refs/heads/master' steps: - - name: Setup .NET Core 3.1 - uses: actions/setup-dotnet@v1 + - name: Setup .NET Core 6 + uses: actions/setup-dotnet@v4 with: - dotnet-version: 3.1.406 + dotnet-version: '6.x' - name: Download Package Files uses: actions/download-artifact@v2 diff --git a/src/ProjNet/CoordinateSystemServices.cs b/src/ProjNet/CoordinateSystemServices.cs index 84c992f..94fd8e3 100644 --- a/src/ProjNet/CoordinateSystemServices.cs +++ b/src/ProjNet/CoordinateSystemServices.cs @@ -305,6 +305,11 @@ public ICoordinateTransformation CreateTransformation(CoordinateSystem source, C return _ctFactory.CreateFromCoordinateSystems(source, target); } + /// + /// AddCoordinateSystem + /// + /// + /// protected void AddCoordinateSystem(int srid, CoordinateSystem coordinateSystem) { lock (((IDictionary) _csBySrid).SyncRoot) @@ -332,6 +337,11 @@ protected void AddCoordinateSystem(int srid, CoordinateSystem coordinateSystem) } } + /// + /// AddCoordinateSystem + /// + /// + /// protected virtual int AddCoordinateSystem(CoordinateSystem coordinateSystem) { int srid = (int) coordinateSystem.AuthorityCode; @@ -340,11 +350,17 @@ protected virtual int AddCoordinateSystem(CoordinateSystem coordinateSystem) return srid; } + /// + /// Clear + /// protected void Clear() { _csBySrid.Clear(); } + /// + /// Count + /// protected int Count { get @@ -354,11 +370,21 @@ protected int Count } } + /// + /// RemoveCoordinateSystem + /// + /// + /// + /// public bool RemoveCoordinateSystem(int srid) { throw new NotSupportedException(); } + /// + /// GetEnumerator + /// + /// public IEnumerator> GetEnumerator() { _initialization.WaitOne(); diff --git a/src/ProjNet/CoordinateSystems/CoordinateSystemFactory.cs b/src/ProjNet/CoordinateSystems/CoordinateSystemFactory.cs index 044bb6c..28d831d 100644 --- a/src/ProjNet/CoordinateSystems/CoordinateSystemFactory.cs +++ b/src/ProjNet/CoordinateSystems/CoordinateSystemFactory.cs @@ -27,11 +27,8 @@ namespace ProjNet.CoordinateSystems /// /// /// CoordinateSystemFactory allows applications to make coordinate systems that - /// cannot be created by a . This factory is very - /// flexible, whereas the authority factory is easier to use. - /// So can be used to make 'standard' coordinate - /// systems, and can be used to make 'special' - /// coordinate systems. + /// is very flexible, whereas the other factories are easier to use. + /// So this Factory can be used to make 'special' coordinate systems. /// For example, the EPSG authority has codes for USA state plane coordinate systems /// using the NAD83 datum, but these coordinate systems always use meters. EPSG does not /// have codes for NAD83 state plane coordinate systems that use feet units. This factory diff --git a/src/ProjNet/CoordinateSystems/Info.cs b/src/ProjNet/CoordinateSystems/Info.cs index edab9a5..56caa3e 100644 --- a/src/ProjNet/CoordinateSystems/Info.cs +++ b/src/ProjNet/CoordinateSystems/Info.cs @@ -35,7 +35,7 @@ public abstract class Info : IInfo /// were specified in the Simple Features interfaces, so they have been kept here. /// This specification does not dictate what the contents of these items /// should be. However, the following guidelines are suggested: - /// When is used to create an object, the ‘Authority’ + /// When is used to create an object, the ‘Authority’ /// and 'AuthorityCode' values should be set to the authority name of the factory object, and the authority /// code supplied by the client, respectively. The other values may or may not be set. (If the authority is /// EPSG, the implementer may consider using the corresponding metadata values in the EPSG tables.) diff --git a/src/ProjNet/CoordinateSystems/Projections/MapProjection.cs b/src/ProjNet/CoordinateSystems/Projections/MapProjection.cs index a369567..4dac204 100644 --- a/src/ProjNet/CoordinateSystems/Projections/MapProjection.cs +++ b/src/ProjNet/CoordinateSystems/Projections/MapProjection.cs @@ -49,10 +49,19 @@ namespace ProjNet.CoordinateSystems.Projections [Serializable] public abstract class MapProjection : MathTransform, IProjection { + /// + /// EPS10 => 1e-10. + /// protected const double EPS10 = 1e-10; + /// + /// EPS7 => 1e-7. + /// protected const double EPS7 = 1e-7; + /// + /// HUGE_VAL => double.NaN. + /// protected const double HUGE_VAL = double.NaN; // ReSharper disable InconsistentNaming @@ -1154,6 +1163,11 @@ protected static double LatitudeToRadians(double y, bool edge) private const double P20 = 0.01677689594356261023; /* 761 / 45360 */ + /// + /// authset + /// + /// + /// protected static double[] authset(double es) { double[] APA = new double[3]; @@ -1169,6 +1183,12 @@ protected static double[] authset(double es) return APA; } + /// + /// authlat + /// + /// + /// + /// protected static double authlat(double beta, double[] APA) { double t = beta + beta; diff --git a/src/ProjNet/CoordinateSystems/Projections/PolarStereographicProjection.cs b/src/ProjNet/CoordinateSystems/Projections/PolarStereographicProjection.cs new file mode 100644 index 0000000..42692a9 --- /dev/null +++ b/src/ProjNet/CoordinateSystems/Projections/PolarStereographicProjection.cs @@ -0,0 +1,229 @@ +// Copyright 2015 +// +// This file is part of ProjNet. +// ProjNet is free software; you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// ProjNet is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. + +// You should have received a copy of the GNU Lesser General Public License +// along with ProjNet; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +/* + * Copyright (C) 2002 Urban Science Applications, Inc. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + */ + +using System; +using System.Collections.Generic; +using ProjNet.CoordinateSystems.Transformations; + +namespace ProjNet.CoordinateSystems.Projections +{ + /// + /// Implements the Polar Stereographic Projection. + /// + [Serializable] + internal class PolarStereographicProjection : MapProjection + { + private readonly double _globalScale; + private readonly double _reciprocGlobalScale; + + private static int MAXIMUM_ITERATIONS = 15; + private static double ITERATION_TOLERANCE = 1E-14; + private static double EPS15 = 1E-15; + private static double M_HALFPI = 0.5 * Math.PI; + private double phits, akm1; + private bool N_POLE; + + + /// + /// Initializes the PolarStereographicProjection object with the specified parameters. + /// + /// List of parameters to initialize the projection. + /// + /// The parameters this projection expects are listed below. + /// + /// ItemsDescriptions + /// central_meridianThe longitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the longitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0). + /// latitude_of_originThe latitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the latitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0). + /// scale_factorThe factor by which the map grid is reduced or enlarged during the projection process, defined by its value at the natural origin. + /// false_eastingSince the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Easting, FE, is the easting value assigned to the abscissa (east). + /// false_northingSince the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Northing, FN, is the northing value assigned to the ordinate. + /// + /// + public PolarStereographicProjection(IEnumerable parameters) + : this(parameters, null) + { + } + + /// + /// Initializes the PolarStereographicProjection object with the specified parameters. + /// + /// List of parameters to initialize the projection. + /// Inverse projection + /// + /// The parameters this projection expects are listed below. + /// + /// ItemsDescriptions + /// central_meridianThe longitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the longitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0). + /// latitude_of_originThe latitude of the point from which the values of both the geographical coordinates on the ellipsoid and the grid coordinates on the projection are deemed to increment or decrement for computational purposes. Alternatively it may be considered as the latitude of the point which in the absence of application of false coordinates has grid coordinates of (0,0). + /// scale_factorThe factor by which the map grid is reduced or enlarged during the projection process, defined by its value at the natural origin. + /// false_eastingSince the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Easting, FE, is the easting value assigned to the abscissa (east). + /// false_northingSince the natural origin may be at or near the centre of the projection and under normal coordinate circumstances would thus give rise to negative coordinates over parts of the mapped area, this origin is usually given false coordinates which are large enough to avoid this inconvenience. The False Northing, FN, is the northing value assigned to the ordinate. + /// + /// + public PolarStereographicProjection(IEnumerable parameters, PolarStereographicProjection inverse) + : base(parameters, inverse) + { + Name = "Polar_Stereographic"; + + _globalScale = scale_factor * _semiMajor; + _reciprocGlobalScale = 1.0 / _globalScale; + + if (_e == 0.0) throw new Exception("Polar Stereographics: only ellipsoidal formulation"); + N_POLE = (lat_origin > 0.0); // N or S hemisphere + phits = Math.Abs(lat_origin); + + if (Math.Abs(phits - M_HALFPI) < EPS10) + { + double one_p_e = 1.0 + _e; + double one_m_e = 1.0 - _e; + double pow_p = Math.Pow(one_p_e, one_p_e); + double pow_m = Math.Pow(one_m_e, one_m_e); + akm1 = 2.0 / Math.Sqrt(pow_p * pow_m); + } + else + { + double sinphits = Math.Sin(phits); + double cosphits = Math.Cos(phits); + akm1 = cosphits / tsfn(cosphits, sinphits, _e); + + double t = _e * sinphits; + akm1 /= Math.Sqrt(1.0 - t * t); + } + } + + /// + /// Converts coordinates in projected meters to radians. + /// + /// + /// + protected override void MetersToRadians(ref double x, ref double y) + { + x *= _reciprocGlobalScale; + y *= _reciprocGlobalScale; + + if (N_POLE) y = -y; + double rho = Math.Sqrt(x * x + y * y); + double tp = -rho / akm1; + double phi_l = M_HALFPI - 2.0 * Math.Atan(tp); + double halfe = -0.5 * _e; + + double lp_phi = 0.0; + for (int iter = MAXIMUM_ITERATIONS; ;) + { + double sinphi = _e * Math.Sin(phi_l); + double one_p_sinphi = 1.0 + sinphi; + double one_m_sinphi = 1.0 - sinphi; + lp_phi = 2.0 * Math.Atan(tp * Math.Pow(one_p_sinphi / one_m_sinphi, halfe)) + M_HALFPI; + if (Math.Abs(phi_l - lp_phi) < ITERATION_TOLERANCE) + { + break; + } + + phi_l = lp_phi; + if (--iter < 0) + { + throw new Exception("Polar Stereographics doesn't converge"); + } + + } + + if (!N_POLE) lp_phi = -lp_phi; + double lp_lam = (x == 0.0 && y == 0.0) ? 0.0 : Math.Atan2(x, y); + + x = lp_lam + central_meridian; + y = lp_phi; + } + + /// + /// Method to convert a point (lon, lat) in radians to (x, y) in meters + /// + /// The longitude of the point in radians when entering, its x-ordinate in meters after exit. + /// The latitude of the point in radians when entering, its y-ordinate in meters after exit. + protected override void RadiansToMeters(ref double lon, ref double lat) + { + double lp_lam = lon - central_meridian; + double lp_phi = lat; + + double coslam = Math.Cos(lp_lam); + double sinlam = Math.Sin(lp_lam); + + if (!N_POLE) + { + lp_phi = -lp_phi; + coslam = -coslam; + } + + double sinphi = Math.Sin(lp_phi); + double cosphi = Math.Cos(lp_phi); + + double x = (Math.Abs(lp_phi - M_HALFPI) < EPS15) ? 0.0 : akm1 * tsfn(cosphi, sinphi, _e); + lon = x * sinlam * _globalScale; + lat = -x * coslam * _globalScale; + } + + + /// + /// Returns the inverse of this projection. + /// + /// IMathTransform that is the reverse of the current projection. + public override MathTransform Inverse() + { + if (_inverse == null) + { + _inverse = new PolarStereographicProjection(_Parameters.ToProjectionParameter(), this); + } + + return _inverse; + } + + private double tsfn(double cosphi, double sinphi, double e) + { + double t = (sinphi > 0.0) ? cosphi / (1.0 + sinphi) : (1.0 - sinphi) / cosphi; + return Math.Exp(e * Atanh(e * sinphi)) * t; + } + + + /// + /// Atanh - Inverse of Math.Tanh + /// + /// The Math.Atanh is not available for netstandard2.0. + /// + private static double Atanh(double x) + { + return Math.Log((1 + x) / (1 - x)) * 0.5; + } + } +} diff --git a/src/ProjNet/CoordinateSystems/Projections/ProjectionsRegistry.cs b/src/ProjNet/CoordinateSystems/Projections/ProjectionsRegistry.cs index 5b46de5..8c4bd1d 100644 --- a/src/ProjNet/CoordinateSystems/Projections/ProjectionsRegistry.cs +++ b/src/ProjNet/CoordinateSystems/Projections/ProjectionsRegistry.cs @@ -49,6 +49,7 @@ static ProjectionsRegistry() Register("oblique_mercator", typeof(ObliqueMercatorProjection)); Register("oblique_stereographic", typeof(ObliqueStereographicProjection)); Register("orthographic", typeof(OrthographicProjection)); + Register("polar_stereographic", typeof(PolarStereographicProjection)); } /// diff --git a/src/ProjNet/CoordinateSystems/Transformations/GeographicTransform.cs b/src/ProjNet/CoordinateSystems/Transformations/GeographicTransform.cs index 88766d0..bc419a8 100644 --- a/src/ProjNet/CoordinateSystems/Transformations/GeographicTransform.cs +++ b/src/ProjNet/CoordinateSystems/Transformations/GeographicTransform.cs @@ -65,11 +65,17 @@ public override string XML } } + /// + /// DimSource + /// public override int DimSource { get { return SourceGCS.Dimension; } } + /// + /// DimTarget + /// public override int DimTarget { get { return TargetGCS.Dimension; } diff --git a/src/ProjNet/CoordinateSystems/VerticalDatum.cs b/src/ProjNet/CoordinateSystems/VerticalDatum.cs index f2634f4..8705bdb 100644 --- a/src/ProjNet/CoordinateSystems/VerticalDatum.cs +++ b/src/ProjNet/CoordinateSystems/VerticalDatum.cs @@ -24,6 +24,9 @@ public VerticalDatum(DatumType type, string name, string authority, long code, s { } + /// + /// ODN - VerticalDatum + /// public static VerticalDatum ODN { get diff --git a/test/ProjNet.Tests/CoordinateTransformTests.cs b/test/ProjNet.Tests/CoordinateTransformTests.cs index 4ae6adf..cbfaea6 100644 --- a/test/ProjNet.Tests/CoordinateTransformTests.cs +++ b/test/ProjNet.Tests/CoordinateTransformTests.cs @@ -644,6 +644,117 @@ public void TestObliqueStereographicProjection() Assert.AreEqual(Coord2171[1], transformedCoord2171[1], 1); } + [Test] + public void TestUniversalPolarStereographicProjection() + { + //test data from http://epsg.io/transform + double[] Coord4326 = new double[] { 15.00, 73.00 }; + double[] Coord32661 = new double[] { 2491967.01029204, 163954.12194234435 }; + + string wkt4326 = "" + + "GEOGCS[\"WGS 84\"," + + "DATUM[\"WGS_1984\"," + + "SPHEROID[\"WGS 84\",6378137,298.257223563," + + "AUTHORITY[\"EPSG\",\"7030\"]]," + + "AUTHORITY[\"EPSG\",\"6326\"]]," + + "PRIMEM[\"Greenwich\",0," + + "AUTHORITY[\"EPSG\",\"8901\"]]," + + "UNIT[\"degree\",0.01745329251994328," + + "AUTHORITY[\"EPSG\",\"9122\"]]," + + "AUTHORITY[\"EPSG\",\"4326\"]]"; + + string wkt32661 = "" + + "PROJCS[\"WGS 84 / UPS North (N,E)\"," + + "GEOGCS[\"WGS 84\"," + + "DATUM[\"WGS_1984\"," + + "SPHEROID[\"WGS 84\",6378137,298.257223563," + + "AUTHORITY[\"EPSG\",\"7030\"]]," + + "AUTHORITY[\"EPSG\",\"6326\"]]," + + "PRIMEM[\"Greenwich\",0," + + "AUTHORITY[\"EPSG\",\"8901\"]]," + + "UNIT[\"degree\",0.0174532925199433," + + "AUTHORITY[\"EPSG\",\"9122\"]]," + + "AUTHORITY[\"EPSG\",\"4326\"]]," + + "PROJECTION[\"Polar_Stereographic\"]," + + "PARAMETER[\"latitude_of_origin\",90]," + + "PARAMETER[\"central_meridian\",0]," + + "PARAMETER[\"scale_factor\",0.994]," + + "PARAMETER[\"false_easting\",2000000]," + + "PARAMETER[\"false_northing\",2000000]," + + "UNIT[\"metre\",1," + + "AUTHORITY[\"EPSG\",\"9001\"]]," + + "AUTHORITY[\"EPSG\",\"32661\"]]"; + + var cs1 = CoordinateSystemFactory.CreateFromWkt(wkt4326); + var cs2 = CoordinateSystemFactory.CreateFromWkt(wkt32661); + var ctf = new CoordinateTransformationFactory(); + + var ict = ctf.CreateFromCoordinateSystems(cs2, cs1); + var ict2 = ctf.CreateFromCoordinateSystems(cs1, cs2); + double[] transformedCoord4326 = ict.MathTransform.Transform(Coord32661); + double[] transformedCoord32661 = ict2.MathTransform.Transform(Coord4326); + + Assert.AreEqual(Coord4326[0], transformedCoord4326[0], 0.01); + Assert.AreEqual(Coord4326[1], transformedCoord4326[1], 0.01); + Assert.AreEqual(Coord32661[0], transformedCoord32661[0], 1); + Assert.AreEqual(Coord32661[1], transformedCoord32661[1], 1); + } + + [Test] + public void TestAustralianAntarcticPolarStereographicProjection() + { + //test data from http://epsg.io/transform + double[] Coord4326 = new double[] { 15.00, -73.00 }; + double[] Coord3032 = new double[] { 4476201.247377692, 7066975.373300694 }; + + string wkt4326 = "" + + "GEOGCS[\"WGS 84\"," + + "DATUM[\"WGS_1984\"," + + "SPHEROID[\"WGS 84\",6378137,298.257223563," + + "AUTHORITY[\"EPSG\",\"7030\"]]," + + "AUTHORITY[\"EPSG\",\"6326\"]]," + + "PRIMEM[\"Greenwich\",0," + + "AUTHORITY[\"EPSG\",\"8901\"]]," + + "UNIT[\"degree\",0.01745329251994328," + + "AUTHORITY[\"EPSG\",\"9122\"]]," + + "AUTHORITY[\"EPSG\",\"4326\"]]"; + + string wkt3032 = "" + + "PROJCS[\"WGS 84 / Australian Antarctic Polar Stereographic\"," + + "GEOGCS[\"WGS 84\"," + + "DATUM[\"WGS_1984\"," + + "SPHEROID[\"WGS 84\",6378137,298.257223563," + + "AUTHORITY[\"EPSG\",\"7030\"]]," + + "AUTHORITY[\"EPSG\",\"6326\"]]," + + "PRIMEM[\"Greenwich\",0," + + "AUTHORITY[\"EPSG\",\"8901\"]]," + + "UNIT[\"degree\",0.0174532925199433," + + "AUTHORITY[\"EPSG\",\"9122\"]]," + + "AUTHORITY[\"EPSG\",\"4326\"]]," + + "PROJECTION[\"Polar_Stereographic\"]," + + "PARAMETER[\"latitude_of_origin\",-71]," + + "PARAMETER[\"central_meridian\",70]," + + "PARAMETER[\"false_easting\",6000000]," + + "PARAMETER[\"false_northing\",6000000]," + + "UNIT[\"metre\",1," + + "AUTHORITY[\"EPSG\",\"9001\"]]," + + "AUTHORITY[\"EPSG\",\"3032\"]]"; + + var cs1 = CoordinateSystemFactory.CreateFromWkt(wkt4326); + var cs2 = CoordinateSystemFactory.CreateFromWkt(wkt3032); + var ctf = new CoordinateTransformationFactory(); + + var ict = ctf.CreateFromCoordinateSystems(cs2, cs1); + var ict2 = ctf.CreateFromCoordinateSystems(cs1, cs2); + double[] transformedCoord4326 = ict.MathTransform.Transform(Coord3032); + double[] transformedCoord3032 = ict2.MathTransform.Transform(Coord4326); + + Assert.AreEqual(Coord4326[0], transformedCoord4326[0], 0.01); + Assert.AreEqual(Coord4326[1], transformedCoord4326[1], 0.01); + Assert.AreEqual(Coord3032[0], transformedCoord3032[0], 1); + Assert.AreEqual(Coord3032[1], transformedCoord3032[1], 1); + } + [Test] public void TestUnitTransforms() { diff --git a/test/ProjNet.Tests/ProjNET.Tests.csproj b/test/ProjNet.Tests/ProjNET.Tests.csproj index 3ade108..610cb67 100644 --- a/test/ProjNet.Tests/ProjNET.Tests.csproj +++ b/test/ProjNet.Tests/ProjNET.Tests.csproj @@ -2,14 +2,14 @@ - netcoreapp3.1 + net6 true 1701;1702;1591 - +