From aafb55c2354a73274c322e3d827d53c7eac0dbc7 Mon Sep 17 00:00:00 2001
From: Anton Kudryashov <qubabox@mail.ru>
Date: Mon, 21 Nov 2016 17:29:49 +0300
Subject: [PATCH] - removed coordinates from cells

---
 MES_Wind/frmMain.cs                     | 32 +++++-----
 WindStressPRM/Objects/ClimateCell.cs    |  8 +--
 WindStressPRM/Objects/PrognosticCell.cs |  8 +--
 WindStressPRM/Utils/Matrix.cs           | 20 ++++++-
 WindStressPRM/WindStressPRM.cs          | 80 +++++++------------------
 5 files changed, 61 insertions(+), 87 deletions(-)

diff --git a/MES_Wind/frmMain.cs b/MES_Wind/frmMain.cs
index 21b2fb8..1d3ff9d 100644
--- a/MES_Wind/frmMain.cs
+++ b/MES_Wind/frmMain.cs
@@ -91,7 +91,7 @@ namespace MES_Wind
  
         private void btnCalcStress_Click(object sender, EventArgs e)
         {
-            const double eps = 0.0001;
+            const double eps = 2;
             const double RasterMissingValue = -9999;
             //extract prognostic u layer
             IMapRasterLayer uRasterLayer = default(IMapRasterLayer);
@@ -130,7 +130,7 @@ namespace MES_Wind
             int rcountPrognostic = uRasterLayer.DataSet.NumRows;
             int ccountPrognostic = uRasterLayer.DataSet.NumColumns;
             WindStressPRM.Matrix<WindStressPRM.PrognosticCell> prognosticMatrix = new WindStressPRM.Matrix<WindStressPRM.PrognosticCell>();
-            prognosticMatrix.Cells = new WindStressPRM.PrognosticCell[rcountPrognostic, ccountPrognostic];
+            prognosticMatrix.Cells = new WindStressPRM.PrognosticCell[ccountPrognostic, rcountPrognostic];
             prognosticMatrix.Size = new WindStressPRM.CellSize(uRasterLayer.DataSet.CellWidth, uRasterLayer.DataSet.CellHeight);
             prognosticMatrix.AffineCoefficients = uRasterLayer.Bounds.AffineCoefficients;
             // fill cells of prognostic matrix
@@ -138,8 +138,11 @@ namespace MES_Wind
             {
                 for (int j = 0; j < ccountPrognostic; j++)
                 {
-                    Coordinate dummyRCoords = uRasterLayer.Bounds.CellCenter_ToProj(i,j);
-                    WindStressPRM.Coordinate cellCoords =new WindStressPRM.Coordinate(dummyRCoords.X,dummyRCoords.Y);
+                    if (i == 0 && j == 0)
+                    {
+                        Coordinate dummyRCoords = uRasterLayer.Bounds.CellCenter_ToProj(i, j);
+                        prognosticMatrix.Origin = new WindStressPRM.Coordinate(dummyRCoords.X, dummyRCoords.Y);
+                    }
                     double uValue = uRasterLayer.DataSet.Value[i, j];
                     double vValue = vRasterLayer.DataSet.Value[i, j];
                     if (Math.Abs(uValue - RasterMissingValue) < eps) {
@@ -148,15 +151,15 @@ namespace MES_Wind
                     if (Math.Abs(vValue - RasterMissingValue) < eps) {
                         vValue = Double.NaN;
                     }
-                    WindStressPRM.PrognosticCell dummyPrognosticCell = new WindStressPRM.PrognosticCell(cellCoords, uValue, vValue);
-                    prognosticMatrix.Cells[i, j] = dummyPrognosticCell;                        
+                    WindStressPRM.PrognosticCell dummyPrognosticCell = new WindStressPRM.PrognosticCell(uValue, vValue);
+                    prognosticMatrix.Cells[j, i] = dummyPrognosticCell;                        
                 }
             }
             //Now we create climate raster class
             int rowCountClim = clim5RasterLayer.DataSet.NumRows;
             int columnCountClim = clim5RasterLayer.DataSet.NumColumns;
             WindStressPRM.Matrix<WindStressPRM.ClimateCell> climateMatrix = new WindStressPRM.Matrix<WindStressPRM.ClimateCell>();
-            climateMatrix.Cells = new WindStressPRM.ClimateCell[rowCountClim, columnCountClim];
+            climateMatrix.Cells = new WindStressPRM.ClimateCell[columnCountClim, rowCountClim];
             climateMatrix.Size = new WindStressPRM.CellSize(clim5RasterLayer.DataSet.CellWidth, clim5RasterLayer.DataSet.CellHeight);
             climateMatrix.AffineCoefficients = clim5RasterLayer.Bounds.AffineCoefficients;
             // fill cells of climate matrix
@@ -164,8 +167,10 @@ namespace MES_Wind
             {
                 for (int j = 0; j < columnCountClim; j++)
                 {
-                    Coordinate dummyCellCoords = clim15RasterLayer.Bounds.CellCenter_ToProj(i,j);
-                    WindStressPRM.Coordinate dummyClimCoords = new WindStressPRM.Coordinate(dummyCellCoords.X,dummyCellCoords.Y);
+                    if (i == 0 && j == 0) {
+                        Coordinate dummyCellCoords = clim15RasterLayer.Bounds.CellCenter_ToProj(i, j);
+                        climateMatrix.Origin = new WindStressPRM.Coordinate(dummyCellCoords.X, dummyCellCoords.Y);;
+                    }
                     //add or substruct in range 0 - 27 to change what lines will be broken 
                     double clim5 = clim5RasterLayer.DataSet.Value[i, j];
                     double clim10 = clim10RasterLayer.DataSet.Value[i, j] +1;
@@ -182,8 +187,8 @@ namespace MES_Wind
                     {
                         clim15 = Double.NaN;
                     }
-                    WindStressPRM.ClimateCell dummyClim = new WindStressPRM.ClimateCell(dummyClimCoords, clim5, clim10, clim15);
-                    climateMatrix.Cells[i, j] = dummyClim;
+                    WindStressPRM.ClimateCell dummyClim = new WindStressPRM.ClimateCell(clim5, clim10, clim15);
+                    climateMatrix.Cells[j, i] = dummyClim;
                 }
             }
             // create PRM_line list to pass to PRM_wind from loaded line layer
@@ -265,9 +270,8 @@ namespace MES_Wind
             {
                 Coordinate coords = new Coordinate(disabledStation.Coordinate.X, disabledStation.Coordinate.Y);
                 DotSpatial.Topology.Point disabledPoint = new DotSpatial.Topology.Point(coords);
-                Feature disabledPointFeature = new Feature(disabledPoint);
-                disabledPointFeature.DataRow[stypecolumn] = disabledStation.Stationtype;
-                disabledPointSet.AddFeature(disabledPointFeature);
+                IFeature disabledPointFeature = disabledPointSet.AddFeature(disabledPoint);
+                disabledPointFeature.DataRow.SetField(stypecolumn, disabledStation.Stationtype);
             }
             //add result layer to the map
             IMapPointLayer resultPointLayer = (MapPointLayer)map1.Layers.Add(disabledPointSet); 
diff --git a/WindStressPRM/Objects/ClimateCell.cs b/WindStressPRM/Objects/ClimateCell.cs
index 446828d..ecd01d7 100644
--- a/WindStressPRM/Objects/ClimateCell.cs
+++ b/WindStressPRM/Objects/ClimateCell.cs
@@ -27,20 +27,14 @@ namespace WindStressPRM
         /// </summary>
         public double Wind15 { get; private set; }
         /// <summary>
-        /// Cell center coordinate pair 
-        /// Координаты центра ячейки
-        /// </summary>
-        public Coordinate Coordinate { get; private set; }
-        /// <summary>
         /// designated constructor
         /// </summary>
         /// <param name="coord"></param>
         /// <param name="w5"></param>
         /// <param name="w10"></param>
         /// <param name="w15"></param>
-        public ClimateCell(Coordinate coord, double w5, double w10, double w15)
+        public ClimateCell(double w5, double w10, double w15)
         {
-            this.Coordinate = coord;
             this.Wind5 = w5;
             this.Wind10 = w10;
             this.Wind15 = w15;
diff --git a/WindStressPRM/Objects/PrognosticCell.cs b/WindStressPRM/Objects/PrognosticCell.cs
index b1fb992..320d794 100644
--- a/WindStressPRM/Objects/PrognosticCell.cs
+++ b/WindStressPRM/Objects/PrognosticCell.cs
@@ -22,19 +22,13 @@ namespace WindStressPRM
         /// </summary>
         public double VelocityY { get; private set; }
         /// <summary>
-        /// Cell center coordinates
-        /// Координаты центра ячейки
-        /// </summary>
-        public Coordinate Coordinate { get; private set; }
-        /// <summary>
         /// designated constructor
         /// </summary>
         /// <param name="coord"> Coordinate pair</param>
         /// <param name="vX"> U component</param>
         /// <param name="vY"> V component</param>
-        public PrognosticCell(Coordinate coord, double vX, double vY)
+        public PrognosticCell(double vX, double vY)
         {
-            this.Coordinate = coord;
             this.VelocityX = vX;
             this.VelocityY = vY;
             if (!(this.CheckValue()))
diff --git a/WindStressPRM/Utils/Matrix.cs b/WindStressPRM/Utils/Matrix.cs
index 1019fd7..343ab9f 100644
--- a/WindStressPRM/Utils/Matrix.cs
+++ b/WindStressPRM/Utils/Matrix.cs
@@ -18,10 +18,28 @@ namespace WindStressPRM
         /// </summary>
         /// <param name="index"></param>
         /// <returns></returns>
-        public Coordinate CoordinateAt(Index index) {
+        public Coordinate CellToProjection(Index index) {
             return new Coordinate(Origin.X + index.Row*Size.Width, Origin.Y - index.Col*Size.Height);
         }
         /// <summary>
+        /// get index of cell for coordinate
+        /// </summary>
+        /// <param name="coordinate"></param>
+        /// <returns></returns>
+        public Index ProjectionToCell(Coordinate coordinate)
+        {
+            double rwDiff = (coordinate.X - Origin.X) / Size.Width;
+            double colDiff = (coordinate.Y - Origin.Y) / -Size.Height;
+            int iRow = (int)Math.Round(rwDiff);
+            int iCol = (int)Math.Round(colDiff);
+
+            if (iRow < 0 || iCol < 0 || iRow >= RowsCount() || iCol >= ColumnCount())
+            {
+                throw new Exception("projectionToCell method trying to find uncorrect index");
+            }
+            return new Index(iRow, iCol);            
+        }
+        /// <summary>
         /// AffineCoefficients for matrix from raster projection
         /// </summary>
         public double[] AffineCoefficients { get; set; }
diff --git a/WindStressPRM/WindStressPRM.cs b/WindStressPRM/WindStressPRM.cs
index cc8a808..fa46b99 100644
--- a/WindStressPRM/WindStressPRM.cs
+++ b/WindStressPRM/WindStressPRM.cs
@@ -278,62 +278,6 @@ namespace WindStressPRM
             return result;
         }
 
-        private double[] affineCoefficients(FunctionType functionType)
-        {
-            switch (functionType)
-            {
-                case FunctionType.FunctionVelocityX:
-                case FunctionType.FunctionVelocityY:
-                    {
-                        return Input.PrognosticCells.AffineCoefficients;
-                    }
-                case FunctionType.FunctionClimate5:
-                case FunctionType.FunctionClimate10:
-                case FunctionType.FunctionClimate15:
-                    {
-                        return Input.ClimateCells.AffineCoefficients;
-                    }
-                default:
-                    break;
-            }
-            return null;
-        }
-
-        // Taken from DotSpatial https://github.com/ViceIce/DotSpatial/blob/22c156c7646b1595d88d2523c066a9c6ab4d3a53/DotSpatial.Data/RasterBoundsExt.cs
-        // RasterBoundsExt.cs. Define AffineCoefficients like this.
-        private Index projectionToCell(Coordinate coordinate, FunctionType functionType)
-        {
-            double[] c = affineCoefficients(functionType);
-            double rw, cl;
-            if (c[2] == 0 && c[4] == 0)
-            {
-                // If we don't have skew terms, then the x and y terms are independent.
-                // also both of the other tests will have divide by zero errors and fail.
-                cl = (coordinate.X - c[0]) / c[1];
-                rw = (coordinate.Y - c[3]) / c[5];
-            }
-            else if (c[2] != 0)
-            {
-                // this equation will have divide by zero if c[2] is zero, but works if c[4] is zero
-                double div = (c[5] * c[1]) / c[2] - c[4];
-                cl = (c[3] + (c[5] * coordinate.X) / c[2] - (c[5] * c[0]) / c[2] - coordinate.Y) / div;
-                rw = (coordinate.X - c[0] - c[1] * cl) / c[2];
-            }
-            else
-            {
-                double div = (c[4] * c[2] / c[1] - c[5]);
-                rw = (c[3] + c[4] * coordinate.X / c[1] - c[4] * c[0] / c[1] - coordinate.Y) / div;
-                cl = (coordinate.X - c[2] * rw - c[0]) / c[1];
-            }
-            int iRow = (int)Math.Round(rw);
-            int iCol = (int)Math.Round(cl);
-
-            if (iRow < 0 || iCol < 0 || iRow >= countInList(functionType, true) || iCol >= countInList(functionType, false))
-            {
-                throw new Exception("projectionToCell method trying to find uncorrect index");
-            }
-            return new Index(iRow, iCol);
-        }
         private Coordinate cellToProjection(Index index, FunctionType functionType)
         {
             switch (functionType)
@@ -341,13 +285,13 @@ namespace WindStressPRM
                 case FunctionType.FunctionVelocityX:
                 case FunctionType.FunctionVelocityY:
                     {
-                        return Input.PrognosticCells.Cells[index.Row, index.Col].Coordinate;
+                        return Input.PrognosticCells.CellToProjection(index);
                     }
                 case FunctionType.FunctionClimate5:
                 case FunctionType.FunctionClimate10:
                 case FunctionType.FunctionClimate15:
                     {
-                        return Input.ClimateCells.Cells[index.Row, index.Col].Coordinate;
+                        return Input.ClimateCells.CellToProjection(index);
                     }
                 default:
                     break;
@@ -427,6 +371,26 @@ namespace WindStressPRM
             throw new Exception("There is no cell size");
         }
 
+        private Index projectionToCell(Coordinate coords, FunctionType functionType) 
+        {
+            switch(functionType) {
+                case FunctionType.FunctionClimate5:
+                case FunctionType.FunctionClimate10:
+                case FunctionType.FunctionClimate15:
+                    {
+                        return Input.ClimateCells.ProjectionToCell(coords);
+                    }
+                case FunctionType.FunctionVelocityX:
+                case FunctionType.FunctionVelocityY:
+                    {
+                        return Input.PrognosticCells.ProjectionToCell(coords);
+                    }
+                default:
+                    break;
+            }
+            throw new Exception("There is no associated matrix");
+        }
+
         private double interpol(Coordinate coords, FunctionType functionType)
         {
             // select directions for projections
-- 
GitLab