From 246879e5c79467d271b5bac21542f72b88a1ca5e Mon Sep 17 00:00:00 2001 From: Stig Rune Jensen Date: Fri, 2 Dec 2022 12:05:27 +0100 Subject: [PATCH] Add polar stereographic projection --- .../PolarStereographicProjection.cs | 218 ++++++++++++++++++ .../Projections/ProjectionsRegistry.cs | 1 + .../ProjNet.Tests/CoordinateTransformTests.cs | 111 +++++++++ 3 files changed, 330 insertions(+) create mode 100644 src/ProjNet/CoordinateSystems/Projections/PolarStereographicProjection.cs diff --git a/src/ProjNet/CoordinateSystems/Projections/PolarStereographicProjection.cs b/src/ProjNet/CoordinateSystems/Projections/PolarStereographicProjection.cs new file mode 100644 index 0000000..703f13a --- /dev/null +++ b/src/ProjNet/CoordinateSystems/Projections/PolarStereographicProjection.cs @@ -0,0 +1,218 @@ +// 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 * Math.Atanh(e * sinphi)) * t; + } + } +} 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/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() {