- refactored interpolation in PRMModule with 1 issue for DotSpatial AffineCoefficients

MES_Wind/PRM_wind.cs

@@ -5,6 +5,17 @@ using System.Text;
 
 namespace MES_Wind
 {
+    public class PRM_index
+    {
+        public int Row;
+        public int Col;
+        public PRM_index(int Row, int Col)
+        {
+            this.Row = Row;
+            this.Col = Col;
+        }
+        public PRM_index() : this(0, 0) { }
+    }
     public class PRM_coordinate
     {
         public double X;
@@ -77,22 +88,32 @@ namespace MES_Wind
 
         
     }
+
+    enum PRMFunctionType
+    { 
+        PRMFunctionVelocityX = 0,
+        PRMFunctionVelocityY = 1,
+        PRMFunctionClimate5 = 2,
+        PRMFunctionClimate10 = 3,
+        PRMFunctionClimate15 = 4
+    }
     class PRM_wind
     {
         //prognistic raster info
         List<List<PRM_raster_cell_prognostic>> prognostic_cells;
         PRM_cell_size prognostic_cellsize;
+        double[] prognostic_AffineCoefficients;
         
         //climate raster info
         List<List<PRM_raster_cell_climate>> climate_cells;
         PRM_cell_size climate_cellsize;
+        double[] climate_AffineCoefficients;
 
         //lines collection
         List<PRM_Line> powerlines;
 
-        #region "control parameters"
+        // must be implemented. AffineCoefficients like in DotSpatial using for projections
         public double dist_threshold;
-        #endregion
 
 
         List<PRM_Line> brokenPowerLinesAfterCheck()
@@ -160,9 +181,12 @@ namespace MES_Wind
             }
             else
             {
-                double uwind = 0; // interpol(coords, Uwind_raster);
-                double vwind = 0; // interpol(coords, Vwind_raster);
-                double climwind = 0; // interpol(coords, clim_layer);
+                // TODO Implement climate15
+                PRMFunctionType climateType = useClimate10 ? PRMFunctionType.PRMFunctionClimate10 : PRMFunctionType.PRMFunctionClimate5;
+                PRM_coordinate coords = new PRM_coordinate((coord1.X + coord2.X) / 2, (coord1.Y + coord2.Y) / 2);
+                double uwind = interpol(coords, PRMFunctionType.PRMFunctionVelocityX);
+                double vwind = interpol(coords, PRMFunctionType.PRMFunctionVelocityY);
+                double climwind = interpol(coords, climateType);
                 double umod = Math.Sqrt(uwind * uwind + vwind * vwind); ;
                 double angleline = Math.Atan2((coord2.Y - coord1.Y), (coord2.X - coord1.X));
                 double anglewind = Math.Atan2(vwind, uwind) - angleline;
@@ -192,6 +216,207 @@ namespace MES_Wind
             return false;
         }
 
+        double[] affineCoefficients(PRMFunctionType functionType)
+        {
+            switch (functionType)
+            {
+                case PRMFunctionType.PRMFunctionVelocityX:
+                case PRMFunctionType.PRMFunctionVelocityY:
+                    {
+                        return prognostic_AffineCoefficients;
+                    }
+                case PRMFunctionType.PRMFunctionClimate5:
+                case PRMFunctionType.PRMFunctionClimate10:
+                case PRMFunctionType.PRMFunctionClimate15:
+                    {
+                        return climate_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.
+        PRM_index projectionToCell(PRM_coordinate coordinate, PRMFunctionType 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))
+            {
+                return new PRM_index();
+            }
 
+            return new PRM_index(iRow, iCol);
+        }
+        PRM_coordinate cellToProjection(PRM_index index, PRMFunctionType functionType)
+        {
+            switch (functionType) {
+                case PRMFunctionType.PRMFunctionVelocityX:
+                case PRMFunctionType.PRMFunctionVelocityY:
+                    {
+                        return prognostic_cells[index.Row][index.Col].coords;
+                    }
+                case PRMFunctionType.PRMFunctionClimate5:
+                case PRMFunctionType.PRMFunctionClimate10:
+                case PRMFunctionType.PRMFunctionClimate15:
+                    {
+                        return climate_cells[index.Row][index.Col].coords;
+                    }
+                default:
+                    break;
+            }
+            return null;
+        }
+
+        int countInList(PRMFunctionType functionType, bool forRows) 
+        {
+            switch (functionType)
+            {
+                case PRMFunctionType.PRMFunctionVelocityX:
+                case PRMFunctionType.PRMFunctionVelocityY:
+                    {
+                        return forRows ? prognostic_cells.Count : prognostic_cells[0].Count;
+                    }
+                case PRMFunctionType.PRMFunctionClimate5:
+                case PRMFunctionType.PRMFunctionClimate10:
+                case PRMFunctionType.PRMFunctionClimate15:
+                    {
+                        return forRows ? climate_cells.Count : climate_cells[0].Count;
+                    }
+                default:
+                    break;
+            }
+            return 0;
+        }
+
+        double valueForFunction(PRMFunctionType functionType, PRM_index index)
+        {
+            switch (functionType)
+            {
+                case PRMFunctionType.PRMFunctionVelocityX:
+                    {
+                        return prognostic_cells[index.Row][index.Col].velocityX;
+                    }
+                case PRMFunctionType.PRMFunctionVelocityY:
+                    {
+                        return prognostic_cells[index.Row][index.Col].velocityY;
+                    }
+                case PRMFunctionType.PRMFunctionClimate5:
+                    {
+                        return climate_cells[index.Row][index.Col].wind5;
+                    }
+                case PRMFunctionType.PRMFunctionClimate10:
+                    {
+                        return climate_cells[index.Row][index.Col].wind10;
+                    }
+                case PRMFunctionType.PRMFunctionClimate15:
+                    {
+                        return climate_cells[index.Row][index.Col].wind15;
+                    }
+                default:
+                    break;
+            }
+            return 0;
+        }
+
+        PRM_cell_size cellSizeForFunction(PRMFunctionType functionType) 
+        {
+            switch (functionType)
+            {
+                case PRMFunctionType.PRMFunctionVelocityX:
+                case PRMFunctionType.PRMFunctionVelocityY:
+                    {
+                        return prognostic_cellsize;
+                    }
+                case PRMFunctionType.PRMFunctionClimate5:
+                case PRMFunctionType.PRMFunctionClimate10:
+                case PRMFunctionType.PRMFunctionClimate15:
+                    {
+                        return climate_cellsize;
+                    }
+                default:
+                    break;
+            }
+            throw new Exception("There is no cell size");
+        }
+
+        double interpol(PRM_coordinate coords, PRMFunctionType functionType)
+        {
+            // select directions for projections
+            const bool normalX = true;// true - East, false West
+            const bool normalY = false;// true - North, false South
+            PRM_index rc = projectionToCell(coords, functionType);
+            PRM_coordinate center = cellToProjection(rc, functionType);
+            double xDiff = coords.X - center.X;
+            double yDiff = coords.Y - center.Y;
+            //calculate second index
+            int row2, col2;
+            if ((xDiff >= 0 && normalX) || (!normalX && xDiff < 0))
+            {
+                row2 = rc.Row >= countInList(functionType, true) - 1 ? rc.Row - 1 : rc.Row + 1;
+            }
+            else
+            {
+                row2 = rc.Row > 0 ? rc.Row - 1 : rc.Row + 1;
+            }
+            if ((yDiff >= 0 && normalY) || (!normalY && yDiff < 0))
+            {
+                col2 = rc.Col >= countInList(functionType, false) - 1 ? rc.Col - 1 : rc.Col + 1;
+            }
+            else
+            {
+                col2 = rc.Col > 0 ? rc.Col - 1 : rc.Col + 1;
+            }
+            // indexes and values at bounds
+            PRM_index rcBotLeft = new PRM_index(Math.Min(row2, rc.Row), Math.Min(col2, rc.Col));
+            PRM_index rcBotRight = new PRM_index(Math.Max(row2, rc.Row), Math.Min(col2, rc.Col));
+            PRM_index rcTopLeft = new PRM_index(Math.Min(row2, rc.Row), Math.Max(col2, rc.Col));
+            PRM_index rcTopRight = new PRM_index(Math.Max(row2, rc.Row), Math.Max(col2, rc.Col));
+            double valBotLeft = valueForFunction(functionType, rcBotLeft);
+            double valBotRight = valueForFunction(functionType, rcBotRight);
+            double valTopLeft = valueForFunction(functionType, rcTopLeft);
+            double valTopRight = valueForFunction(functionType, rcTopRight);
+            PRM_coordinate origin = cellToProjection(rcBotLeft, functionType);
+            //PRM_coordinate last = cellToProjection(rcTopRight, functionType);//test only
+            // sizes for cell
+            double hx = cellSizeForFunction(functionType).width;
+            double hy = cellSizeForFunction(functionType).height;
+            // coefficients
+            double px = (coords.X - origin.X) / hx;
+            double py = (coords.Y - origin.Y) / hy;
+            // inverse directions
+            px *= normalX ? 1 : -1;
+            py *= normalY ? 1 : -1;
+            // interpolation
+            double top = (1 - px) * valTopLeft + px * valTopRight;
+            double bot = (1 - px) * valBotLeft + px * valBotRight;
+            double rval = (1 - py) * bot + py * top;
+            return rval;
+        }
     }
 }