From b3e306cf2a24c5868832e3d311108eb728bf3a17 Mon Sep 17 00:00:00 2001 From: Christoph Perger Date: Wed, 22 Jun 2016 13:10:12 +0200 Subject: [PATCH] Fixed LAEA projection y coordinate that resulted in NaN. Fixes #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 | 1 + Contributors | 1 + .../Projected/Europe.cs | 24 +++++++++++++++++++ .../Transforms/LambertAzimuthalEqualArea.cs | 5 ++-- 4 files changed, 29 insertions(+), 2 deletions(-) diff --git a/Changelog.md b/Changelog.md index 2bffe4926..d80adf9f8 100644 --- a/Changelog.md +++ b/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) diff --git a/Contributors b/Contributors index c36f4ae74..d88a72c59 100644 --- a/Contributors +++ b/Contributors @@ -17,3 +17,4 @@ Maxim Miroshnikov Contributors: Martin Karing +Christoph Perger diff --git a/Source/DotSpatial.Projections.Tests/Projected/Europe.cs b/Source/DotSpatial.Projections.Tests/Projected/Europe.cs index 24acb0f28..5dacf0e6f 100644 --- a/Source/DotSpatial.Projections.Tests/Projected/Europe.cs +++ b/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 GetProjections() { return ProjectionInfoDesc.GetForCoordinateSystemCategory(KnownCoordinateSystems.Projected.Europe); diff --git a/Source/DotSpatial.Projections/Transforms/LambertAzimuthalEqualArea.cs b/Source/DotSpatial.Projections/Transforms/LambertAzimuthalEqualArea.cs index 8da64d9f5..30c781984 100644 --- a/Source/DotSpatial.Projections/Transforms/LambertAzimuthalEqualArea.cs +++ b/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; }