|
| 1 | +// Copyright 2015 |
| 2 | +// |
| 3 | +// This file is part of ProjNet. |
| 4 | +// ProjNet is free software; you can redistribute it and/or modify |
| 5 | +// it under the terms of the GNU Lesser General Public License as published by |
| 6 | +// the Free Software Foundation; either version 2 of the License, or |
| 7 | +// (at your option) any later version. |
| 8 | +// |
| 9 | +// ProjNet is distributed in the hope that it will be useful, |
| 10 | +// but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 11 | +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 12 | +// GNU Lesser General Public License for more details. |
| 13 | + |
| 14 | +// You should have received a copy of the GNU Lesser General Public License |
| 15 | +// along with ProjNet; if not, write to the Free Software |
| 16 | +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
| 17 | + |
| 18 | +/* |
| 19 | + * Copyright (C) 2002 Urban Science Applications, Inc. |
| 20 | + * |
| 21 | + * This library is free software; you can redistribute it and/or |
| 22 | + * modify it under the terms of the GNU Lesser General Public |
| 23 | + * License as published by the Free Software Foundation; either |
| 24 | + * version 2.1 of the License, or (at your option) any later version. |
| 25 | + * |
| 26 | + * This library is distributed in the hope that it will be useful, |
| 27 | + * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 28 | + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 29 | + * Lesser General Public License for more details. |
| 30 | + * |
| 31 | + * You should have received a copy of the GNU Lesser General Public |
| 32 | + * License along with this library; if not, write to the Free Software |
| 33 | + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
| 34 | + * |
| 35 | + */ |
| 36 | + |
| 37 | +using System; |
| 38 | +using System.Collections.Generic; |
| 39 | +using ProjNet.CoordinateSystems.Transformations; |
| 40 | + |
| 41 | +namespace ProjNet.CoordinateSystems.Projections |
| 42 | +{ |
| 43 | + /// <summary> |
| 44 | + /// Implements the Polar Stereographic Projection. |
| 45 | + /// </summary> |
| 46 | + [Serializable] |
| 47 | + internal class PolarStereographicProjection : MapProjection |
| 48 | + { |
| 49 | + private readonly double _globalScale; |
| 50 | + private readonly double _reciprocGlobalScale; |
| 51 | + |
| 52 | + private static int MAXIMUM_ITERATIONS = 15; |
| 53 | + private static double ITERATION_TOLERANCE = 1E-14; |
| 54 | + private static double EPS15 = 1E-15; |
| 55 | + private static double M_HALFPI = 0.5 * Math.PI; |
| 56 | + private double phits, akm1; |
| 57 | + private bool N_POLE; |
| 58 | + |
| 59 | + |
| 60 | + /// <summary> |
| 61 | + /// Initializes the PolarStereographicProjection object with the specified parameters. |
| 62 | + /// </summary> |
| 63 | + /// <param name="parameters">List of parameters to initialize the projection.</param> |
| 64 | + /// <remarks> |
| 65 | + /// <para>The parameters this projection expects are listed below.</para> |
| 66 | + /// <list type="table"> |
| 67 | + /// <listheader><term>Items</term><description>Descriptions</description></listheader> |
| 68 | + /// <item><term>central_meridian</term><description>The 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).</description></item> |
| 69 | + /// <item><term>latitude_of_origin</term><description>The 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).</description></item> |
| 70 | + /// <item><term>scale_factor</term><description>The factor by which the map grid is reduced or enlarged during the projection process, defined by its value at the natural origin.</description></item> |
| 71 | + /// <item><term>false_easting</term><description>Since 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).</description></item> |
| 72 | + /// <item><term>false_northing</term><description>Since 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.</description></item> |
| 73 | + /// </list> |
| 74 | + /// </remarks> |
| 75 | + public PolarStereographicProjection(IEnumerable<ProjectionParameter> parameters) |
| 76 | + : this(parameters, null) |
| 77 | + { |
| 78 | + } |
| 79 | + |
| 80 | + /// <summary> |
| 81 | + /// Initializes the PolarStereographicProjection object with the specified parameters. |
| 82 | + /// </summary> |
| 83 | + /// <param name="parameters">List of parameters to initialize the projection.</param> |
| 84 | + /// <param name="inverse">Inverse projection</param> |
| 85 | + /// <remarks> |
| 86 | + /// <para>The parameters this projection expects are listed below.</para> |
| 87 | + /// <list type="table"> |
| 88 | + /// <listheader><term>Items</term><description>Descriptions</description></listheader> |
| 89 | + /// <item><term>central_meridian</term><description>The 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).</description></item> |
| 90 | + /// <item><term>latitude_of_origin</term><description>The 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).</description></item> |
| 91 | + /// <item><term>scale_factor</term><description>The factor by which the map grid is reduced or enlarged during the projection process, defined by its value at the natural origin.</description></item> |
| 92 | + /// <item><term>false_easting</term><description>Since 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).</description></item> |
| 93 | + /// <item><term>false_northing</term><description>Since 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.</description></item> |
| 94 | + /// </list> |
| 95 | + /// </remarks> |
| 96 | + public PolarStereographicProjection(IEnumerable<ProjectionParameter> parameters, PolarStereographicProjection inverse) |
| 97 | + : base(parameters, inverse) |
| 98 | + { |
| 99 | + Name = "Polar_Stereographic"; |
| 100 | + |
| 101 | + _globalScale = scale_factor * _semiMajor; |
| 102 | + _reciprocGlobalScale = 1.0 / _globalScale; |
| 103 | + |
| 104 | + if (_e == 0.0) throw new Exception("Polar Stereographics: only ellipsoidal formulation"); |
| 105 | + N_POLE = (lat_origin > 0.0); // N or S hemisphere |
| 106 | + phits = Math.Abs(lat_origin); |
| 107 | + |
| 108 | + if (Math.Abs(phits - M_HALFPI) < EPS10) |
| 109 | + { |
| 110 | + double one_p_e = 1.0 + _e; |
| 111 | + double one_m_e = 1.0 - _e; |
| 112 | + double pow_p = Math.Pow(one_p_e, one_p_e); |
| 113 | + double pow_m = Math.Pow(one_m_e, one_m_e); |
| 114 | + akm1 = 2.0 / Math.Sqrt(pow_p * pow_m); |
| 115 | + } |
| 116 | + else |
| 117 | + { |
| 118 | + double sinphits = Math.Sin(phits); |
| 119 | + double cosphits = Math.Cos(phits); |
| 120 | + akm1 = cosphits / tsfn(cosphits, sinphits, _e); |
| 121 | + |
| 122 | + double t = _e * sinphits; |
| 123 | + akm1 /= Math.Sqrt(1.0 - t * t); |
| 124 | + } |
| 125 | + } |
| 126 | + |
| 127 | + /// <summary> |
| 128 | + /// Converts coordinates in projected meters to radians. |
| 129 | + /// </summary> |
| 130 | + /// <param name="x"></param> |
| 131 | + /// <param name="y"></param> |
| 132 | + protected override void MetersToRadians(ref double x, ref double y) |
| 133 | + { |
| 134 | + x *= _reciprocGlobalScale; |
| 135 | + y *= _reciprocGlobalScale; |
| 136 | + |
| 137 | + if (N_POLE) y = -y; |
| 138 | + double rho = Math.Sqrt(x * x + y * y); |
| 139 | + double tp = -rho / akm1; |
| 140 | + double phi_l = M_HALFPI - 2.0 * Math.Atan(tp); |
| 141 | + double halfe = -0.5 * _e; |
| 142 | + |
| 143 | + double lp_phi = 0.0; |
| 144 | + for (int iter = MAXIMUM_ITERATIONS; ;) |
| 145 | + { |
| 146 | + double sinphi = _e * Math.Sin(phi_l); |
| 147 | + double one_p_sinphi = 1.0 + sinphi; |
| 148 | + double one_m_sinphi = 1.0 - sinphi; |
| 149 | + lp_phi = 2.0 * Math.Atan(tp * Math.Pow(one_p_sinphi / one_m_sinphi, halfe)) + M_HALFPI; |
| 150 | + if (Math.Abs(phi_l - lp_phi) < ITERATION_TOLERANCE) |
| 151 | + { |
| 152 | + break; |
| 153 | + } |
| 154 | + |
| 155 | + phi_l = lp_phi; |
| 156 | + if (--iter < 0) |
| 157 | + { |
| 158 | + throw new Exception("Polar Stereographics doesn't converge"); |
| 159 | + } |
| 160 | + |
| 161 | + } |
| 162 | + |
| 163 | + if (!N_POLE) lp_phi = -lp_phi; |
| 164 | + double lp_lam = (x == 0.0 && y == 0.0) ? 0.0 : Math.Atan2(x, y); |
| 165 | + |
| 166 | + x = lp_lam + central_meridian; |
| 167 | + y = lp_phi; |
| 168 | + } |
| 169 | + |
| 170 | + /// <summary> |
| 171 | + /// Method to convert a point (lon, lat) in radians to (x, y) in meters |
| 172 | + /// </summary> |
| 173 | + /// <param name="lon">The longitude of the point in radians when entering, its x-ordinate in meters after exit.</param> |
| 174 | + /// <param name="lat">The latitude of the point in radians when entering, its y-ordinate in meters after exit.</param> |
| 175 | + protected override void RadiansToMeters(ref double lon, ref double lat) |
| 176 | + { |
| 177 | + double lp_lam = lon - central_meridian; |
| 178 | + double lp_phi = lat; |
| 179 | + |
| 180 | + double coslam = Math.Cos(lp_lam); |
| 181 | + double sinlam = Math.Sin(lp_lam); |
| 182 | + |
| 183 | + if (!N_POLE) |
| 184 | + { |
| 185 | + lp_phi = -lp_phi; |
| 186 | + coslam = -coslam; |
| 187 | + } |
| 188 | + |
| 189 | + double sinphi = Math.Sin(lp_phi); |
| 190 | + double cosphi = Math.Cos(lp_phi); |
| 191 | + |
| 192 | + double x = (Math.Abs(lp_phi - M_HALFPI) < EPS15) ? 0.0 : akm1 * tsfn(cosphi, sinphi, _e); |
| 193 | + lon = x * sinlam * _globalScale; |
| 194 | + lat = -x * coslam * _globalScale; |
| 195 | + } |
| 196 | + |
| 197 | + |
| 198 | + /// <summary> |
| 199 | + /// Returns the inverse of this projection. |
| 200 | + /// </summary> |
| 201 | + /// <returns>IMathTransform that is the reverse of the current projection.</returns> |
| 202 | + public override MathTransform Inverse() |
| 203 | + { |
| 204 | + if (_inverse == null) |
| 205 | + { |
| 206 | + _inverse = new PolarStereographicProjection(_Parameters.ToProjectionParameter(), this); |
| 207 | + } |
| 208 | + |
| 209 | + return _inverse; |
| 210 | + } |
| 211 | + |
| 212 | + private double tsfn(double cosphi, double sinphi, double e) |
| 213 | + { |
| 214 | + double t = (sinphi > 0.0) ? cosphi / (1.0 + sinphi) : (1.0 - sinphi) / cosphi; |
| 215 | + return Math.Exp(e * Math.Atanh(e * sinphi)) * t; |
| 216 | + } |
| 217 | + } |
| 218 | +} |
0 commit comments