Fixed LAEA projection y coordinate that resulted in NaN. Fixes DotSpa…

…tial#813

Fix was applied according to Proj.4 C library: used the proj.4 c code in the master branch of the project under: https://github.com/OSGeo/proj.4/blob/master/src/PJ_laea.c
line 130 and line 143

Changelog.md

@@ -18,6 +18,7 @@ All notable changes to this project will be documented in this file.
 - Removed obsolete methods\properties (#797)
 
 ### Fixed
+- Fixed LAEA reprojected y coordinate that resulted in n.def (#813)
 - ShapeReader skipping one entry when switching the page (#774)
 - DotSpatial.Projections dll file is very big (#27)
 - ParseEsriString leaves datums.xml open (#713)

Contributors

@@ -17,3 +17,4 @@ Maxim Miroshnikov <[email protected]>
 
 Contributors:
 Martin Karing <[email protected]>
+Christoph Perger <[email protected]>

Source/DotSpatial.Projections.Tests/Projected/Europe.cs

@@ -18,6 +18,30 @@ public void EuropeProjectedTests(ProjectionInfoDesc pInfo)
             Assert.AreEqual(false, pInfo.ProjectionInfo.IsLatLon);
         }
 
+        [Test]
+        public void ETRS1989LAEA()
+        {
+            ProjectionInfo pStart = KnownCoordinateSystems.Geographic.World.WGS1984;
+            ProjectionInfo pEnd = KnownCoordinateSystems.Projected.Europe.ETRS1989LAEA;
+
+            // Vienna, Austria
+            var lon = 16.4;
+            var lat = 48.2;
+
+            double[] xy = new double[] { lon, lat };
+            double[] z = new double[] { 0 };
+
+            Reproject.ReprojectPoints(xy, z, pStart, pEnd, 0, 1);
+
+            Reproject.ReprojectPoints(xy, z, pEnd, pStart, 0, 1);
+
+            // Test X
+            Assert.AreEqual(16.4, xy[0], 0.01);
+
+            // Test Y
+            Assert.AreEqual(48.2, xy[1], 0.01);
+        }
+
         private static IEnumerable<ProjectionInfoDesc> GetProjections()
         {
             return ProjectionInfoDesc.GetForCoordinateSystemCategory(KnownCoordinateSystems.Projected.Europe);

Source/DotSpatial.Projections/Transforms/LambertAzimuthalEqualArea.cs

@@ -22,6 +22,7 @@
 //        Name         |    Date    |        Comment
 // --------------------|------------|------------------------------------------------------------
 // Ted Dunsford        |   5/3/2010 |  Updated project to DotSpatial.Projection and license to LGPL
+// Christoph Perger    |  22/6/2016 |  Fixed some issues with projection from LAEA to WGS84
 // ********************************************************************************************************
 
 using System;
@@ -189,7 +190,7 @@ protected override void EllipticalInverse(double[] xy, double[] lp, int startInd
                 double cy = xy[y];
                 if (_mode == Modes.Equitorial || _mode == Modes.Oblique)
                 {
-                    double rho = Proj.Hypot(cx /= _dd, cy * _dd);
+                    double rho = Proj.Hypot(cx /= _dd, cy *= _dd);
                     if (rho < EPS10)
                     {
                         lp[lam] = 0;
@@ -201,7 +202,7 @@ protected override void EllipticalInverse(double[] xy, double[] lp, int startInd
                     cx *= (sCe = Math.Sin(sCe));
                     if (_mode == Modes.Oblique)
                     {
-                        ab = cCe * _sinb1 + y * sCe * _cosb1 / rho;
+                        ab = cCe * _sinb1 + cy * sCe * _cosb1 / rho;
                         //q = _qp*(ab);
                         cy = rho * _cosb1 * cCe - cy * _sinb1 * sCe;
                     }