Fixes inconsistent use of affine coefficients. Closes issue DotSpatia…

…l#822

Changelog.md

@@ -25,6 +25,7 @@ All notable changes to this project will be documented in this file.
 - Removed obsolete methods\properties (#797)
 
 ### Fixed
+- inconsistent use of affine coefficients (#822)
 - Fixed the shift in x-coordinate when reprojecting from WGS84 to LAEA (#815)
 - Fixed LAEA reprojected y coordinate that resulted in n.def (#813)
 - ShapeReader skipping one entry when switching the page (#774)

Source/DotSpatial.Data.Rasters.GdalExtension/GdalImage.cs

@@ -124,6 +124,7 @@ private void ReadHeader()
             NumBands = _dataset.RasterCount;
             var test = new double[6];
             _dataset.GetGeoTransform(test);
+            test = (new AffineTransform(test)).TransfromToCorner(0.5, 0.5);//shift origin by half a cell
             ProjectionString = _dataset.GetProjection();
             Bounds = new RasterBounds(Height, Width, test);
             WorldFile.Affine = test;

Source/DotSpatial.Data.Rasters.GdalExtension/GdalImageSource.cs

@@ -184,6 +184,7 @@ private void ReadHeader()
 
             double[] test = new double[6];
             _dataset.GetGeoTransform(test);
+            test = (new AffineTransform(test)).TransfromToCorner(0.5, 0.5);//shift origin by half a cell
             Bounds = new RasterBounds(_dataset.RasterYSize, _dataset.RasterXSize, test);
             _red = _dataset.GetRasterBand(1);
             _numOverviews = _red.GetOverviewCount();

Source/DotSpatial.Data.Rasters.GdalExtension/GdalRaster.cs

@@ -582,6 +582,9 @@ private void ReadHeader()
             }
             double[] affine = new double[6];
             _dataset.GetGeoTransform(affine);
+            // in gdal (row,col) coordinates are defined relative to the top-left corner of the top-left cell
+            // shift them by half a cell to give coordinates relative to the center of the top-left cell
+            affine = (new AffineTransform(affine)).TransfromToCorner(0.5, 0.5);
             ProjectionString = projString;
             Bounds = new RasterBounds(base.NumRows, base.NumColumns, affine);
             PixelSpace = Marshal.SizeOf(typeof(T));

Source/DotSpatial.Data.Rasters.GdalExtension/GdalTiledImage.cs

@@ -68,6 +68,7 @@ private void ReadHeader()
             WorldFile = new WorldFile { Affine = new double[6] };
             double[] test = new double[6];
             _dataset.GetGeoTransform(test);
+            test = (new AffineTransform(test)).TransfromToCorner(0.5, 0.5);//shift origin by half a cell
             Bounds = new RasterBounds(Height, Width, test);
             WorldFile.Affine = test;
             Close();

Source/DotSpatial.Data.Tests/AffineCoefficientsTests.cs

@@ -0,0 +1,71 @@
+using System;
+using System.IO;
+using GeoAPI.Geometries;
+using DotSpatial.Data.Rasters.GdalExtension;
+using DotSpatial.Tests.Common;
+using NUnit.Framework;
+using TestClass = NUnit.Framework.TestFixtureAttribute;
+using TestMethod = NUnit.Framework.TestAttribute;
+using TestCleanup = NUnit.Framework.TearDownAttribute;
+using TestInitialize = NUnit.Framework.SetUpAttribute;
+using ClassCleanup = NUnit.Framework.TestFixtureTearDownAttribute;
+using ClassInitialize = NUnit.Framework.TestFixtureSetUpAttribute;
+
+namespace DotSpatial.Data.Tests
+{
+    /// <summary>
+    ///This is a class testing the implementation of affine coefficients.
+    ///</summary>
+    [TestClass()]
+    public class AffineCoefficientsTests
+    {
+        /// <summary>
+        ///A test for affine coefficients in AffineTransform: generate random points with 
+        ///known (row, column) and check whether AffineTransform.ProjToCell returns the correct values.
+        ///This test fails when using Math.Floor in AffineTransform.ProjToCell.
+        ///</summary>
+        [TestMethod]
+        public void AffineTransformTest()
+        {
+            var c = new double[6];
+            c[0] = -179.75; // x-coordinate of the *center* of the top-left cell
+            c[1] = 0.5;     // column width
+            c[2] = 0.0;     // rotation/skew term
+            c[3] = 89.75;   // y-coordinate of the *center* of the top-left cell
+            c[4] = 0.0;     // rotation/skew term
+            c[5] = -0.5;    // negative row height
+            var at = new AffineTransform(c);
+
+            var rand = new Random();
+
+            for (int row = 0; row < 9; row++)
+            {
+                for (int col = 0; col < 9; col++)
+                {
+                    Coordinate corner = at.CellTopLeft_ToProj(row, col);
+                    double frac_col = rand.NextDouble();
+                    double frac_row = rand.NextDouble();
+                    double dX = c[1] * frac_col + c[2] * frac_row;
+                    double dY = c[4] * frac_col + c[5] * frac_row;
+                    Coordinate point = new Coordinate(corner.X + dX, corner.Y + dY);
+                    Assert.AreEqual(at.ProjToCell(point), new RcIndex(row, col));
+                }
+            }
+        }
+
+        /// <summary>
+        ///A test for affine coefficients in GdalRaster: use gdal to load a raster file with
+        ///known geolocation and test the location of the center of the first grid cell.
+        ///Without a half-cell shift applied to the origin of the grid in GdalRaster.ReadHeader,
+        ///this test fails because affine coefficients are defined differently in gdal and dotspatial.
+        ///</summary>
+        [TestMethod]
+        public void GdalRasterTest()
+        {
+            var rp = new GdalRasterProvider();
+            var raster = rp.Open(@"Data\Grids\sample_geotiff.tif");
+            var at = new AffineTransform(raster.Bounds.AffineCoefficients);
+            Assert.AreEqual(at.CellCenter_ToProj(0, 0), new Coordinate(-179.9499969, 89.9499969));// correct location from sample_geotiff.tfw
+        }
+    }
+}

Source/DotSpatial.Data.Tests/Data/Grids/sample_geotiff.tfw

@@ -0,0 +1,6 @@
+      0.1000000
+      0.0000000
+      0.0000000
+     -0.1000000
+   -179.9499969
+     89.9499969

Source/DotSpatial.Data.Tests/DotSpatial.Data.Tests.csproj

@@ -86,6 +86,7 @@
     </CodeAnalysisDependentAssemblyPaths>
   </ItemGroup>
   <ItemGroup>
+    <Compile Include="AffineCoefficientsTests.cs" />
     <Compile Include="AttributeTableTests.cs" />
     <Compile Include="BinaryRasterProviderTests.cs" />
     <Compile Include="DataManagerTests.cs" />

Source/DotSpatial.Data/AffineTransform.cs

@@ -119,10 +119,10 @@ public Coordinate CellBottomRight_ToProj(int row, int column)
         /// <summary>
         ///  Given the row and column, this returns new affine coefficients transformed to that cell
         /// </summary>
-        /// <param name="startRow">The integer row index from 0 to numRows - 1</param>
-        /// <param name="startColumn">The integer column index from 0 to numColumns - 1</param>
-        /// <returns>Transoformed affine coefficients</returns>
-        public double[] TransfromToCorner(int startColumn, int startRow)
+        /// <param name="startRow">The row index from 0 to numRows - 1</param>
+        /// <param name="startColumn">The column index from 0 to numColumns - 1</param>
+        /// <returns>Transformed affine coefficients</returns>
+        public double[] TransfromToCorner(double startColumn, double startRow)
         {
             // X = [0] + [1] * column + [2] * row;
             // Y = [3] + [4] * column + [5] * row;
@@ -165,8 +165,8 @@ public RcIndex ProjToCell(Coordinate location)
                 rw = (c[3] + c[4] * location.X / c[1] - c[4] * c[0] / c[1] - location.Y) / div;
                 cl = (location.X - c[2] * rw - c[0]) / c[1];
             }
-            var iRow = (int)Math.Floor(rw);
-            var iCol = (int)Math.Floor(cl);
+            var iRow = (int)Math.Round(rw);
+            var iCol = (int)Math.Round(cl);
             return new RcIndex(iRow, iCol);
         }
 

Source/DotSpatial.Plugins.MapWindowProjectFileCompatibility.Tests/DotSpatial.Plugins.MapWindowProjectFileCompatibility.Tests.csproj

@@ -73,6 +73,9 @@
   <ItemGroup>
     <None Include="packages.config" />
   </ItemGroup>
+  <ItemGroup>
+    <Service Include="{82A7F48D-3B50-4B1E-B82E-3ADA8210C358}" />
+  </ItemGroup>
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
   <Import Project="$(SolutionDir)\.nuget\NuGet.targets" Condition="Exists('$(SolutionDir)\.nuget\NuGet.targets')" />
   <PropertyGroup>

Source/DotSpatial.Positioning.Tests/DotSpatial.Positioning.Tests.csproj

@@ -62,6 +62,9 @@
       <Name>DotSpatial.Positioning</Name>
     </ProjectReference>
   </ItemGroup>
+  <ItemGroup>
+    <Service Include="{82A7F48D-3B50-4B1E-B82E-3ADA8210C358}" />
+  </ItemGroup>
   <Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />
   <Import Project="$(SolutionDir)\.nuget\NuGet.targets" Condition="Exists('$(SolutionDir)\.nuget\NuGet.targets')" />
   <!-- To modify your build process, add your task inside one of the targets below and uncomment it.