using System;
using System.Collections.Generic;
namespace WindStressPRM
{
/// <summary>
/// Index class for raster lists in list
/// Индекс-класс для растровых массивов
/// </summary>
public class Index
{
/// <summary>
/// Outer index
/// Внешний индекс
/// </summary>
public int Row;
/// <summary>
/// Inner index
/// Внутренний индекс
/// </summary>
public int Col;
/// <summary>
/// designated constructor
/// </summary>
/// <param name="Row">row index</param>
/// <param name="Col">column index</param>
public Index(int Row, int Col)
{
if (Row <= -1 || Col <= -1) {
throw new Exception("Index must be initialized with nonegative integer value");
}
this.Row = Row;
this.Col = Col;
}
}
/// <summary>
/// Coordinate pair
/// </summary>
public class Coordinate
{
/// <summary>
/// easting coordinate
/// Координата по широте
/// </summary>
public double X;
/// <summary>
/// northing coordinate
/// Координата по долготе
/// </summary>
public double Y;
/// <summary>
/// designated constructor
/// Основной конструктор
/// </summary>
/// <param name="X"></param>
/// <param name="Y"></param>
public Coordinate(double X, double Y)
{
if (!(this.CheckValue()))
{
throw new System.ArgumentException("Passed coordinates are not valid!");
}
else
{
this.X = X;
this.Y = Y;
}
}
/// <summary>
/// Проверка на валидность значений
/// </summary>
/// <returns></returns>
public bool CheckValue()
{
bool checker = true;
if (X != null && Y != null)
{
checker = true;
}
else {
checker = false;
}
return checker;
}
}
/// <summary>
/// Cell obj for regular prognostic wind field
/// Объект ячейки для поля прогностического ветра на высоте 10м на регулярной сетке
/// </summary>
public class PrognosticCell
{
/// <summary>
/// U - component of wind velocity, m/s
/// U - компонента скорости ветра, м/с
/// </summary>
public double VelocityX;
/// <summary>
/// V - component of wind velocity, m/s
/// V - компонента скорости ветра, м/с
/// </summary>
public double VelocityY;
/// <summary>
/// Cell center coordinates
/// Координаты центра ячейки
/// </summary>
public Coordinate Coordinate;
/// <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)
{
if (!(this.CheckValue()))
{
throw new System.ArgumentException("Prognostic wind velocities are incorrect");
}
else
{
this.Coordinate = coord;
this.VelocityX = vX;
this.VelocityY = vY;
}
}
/// <summary>
/// Проверка полей на валидность
/// </summary>
/// <returns></returns>
public bool CheckValue()
{
bool checker = false;
///Скорость ветра на высоте 10м от поверхности на Земле не может превышать 70м/c
if (Math.Abs(VelocityX) < 70 && Math.Abs(VelocityY) < 70)
{
checker = true;
}
else {
checker = false;
}
return checker;
}
}
/// <summary>
/// Cell obj for climate wind regular data
/// Объект - ячейка полей ветра определенной повторяемости на 10м на регулярной сетке
/// </summary>
public class ClimateCell
{
/// <summary>
/// once-in-5-years frequency wind, m/s
/// скорость ветра повторяемости один раз в 5 лет, м/с
/// </summary>
public double Wind5;
/// <summary>
/// once-in-10-years frequency wind, m/s
/// скорость ветра повторяемости один раз в 10 лет, м/с
/// </summary>
public double Wind10;
/// <summary>
/// once-in-15-years frequency wind, m/s
/// скорость ветра повторяемости один раз в 15 лет, м/с
/// </summary>
public double Wind15;
/// <summary>
/// Cell center coordinate pair
/// Координаты центра ячейки
/// </summary>
public Coordinate Coordinate;
/// <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)
{
if (!(this.CheckValue()))
{
throw new System.ArgumentException("Climate wind value is not correct");
}
else
{
this.Coordinate = coord;
this.Wind5 = w5;
this.Wind10 = w10;
this.Wind15 = w15;
}
}
/// <summary>
/// Провекра валидности полей класса
/// </summary>
/// <returns></returns>
public bool CheckValue()
{
bool checker = true;
if (Math.Abs(Wind5) < 70 && Math.Abs(Wind10) < 70 && Math.Abs(Wind15) < 70)
{
checker = true;
}
else {
checker = false;
}
return checker;
}
}
/// <summary>
/// Cell Size parameters
/// Параметры ячейки (регулярной сетки)
/// </summary>
public struct CellSize
{
/// <summary>
/// ширина ячейки (расстояние между соседними по широте центрами ячеек)
/// </summary>
public double Width;
/// <summary>
/// высота ячейки (расстояние между соседними по долготе центрами ячеек)
/// </summary>
public double Height;
/// <summary>
/// designated constructor
/// </summary>
/// <param name="wdh">Cell Width</param>
/// <param name="hgh">Cell Height</param>
public CellSize(double wdh, double hgh)
{
this.Width = wdh;
this.Height = hgh;
if (!(this.CheckValue()))
{
throw new System.ArgumentException("Cell width or height values are incorrect!");
}
}
/// <summary>
/// Проверка валидности полей
/// </summary>
/// <returns></returns>
public bool CheckValue()
{
bool checker = true;
if (Width > 0 && Height > 0)
{
checker = true;
}
else { checker = false; }
return checker;
}
}
/// <summary>
/// power line object
/// объект - ЛЭП
/// </summary>
public class Powerline
{
/// <summary>
/// unique id
/// уникальный идентификатор
/// </summary>
public int Identifier { get; set; }
/// <summary>
/// year of construction
/// год постройки ЛЭП
/// </summary>
public int Year;
/// <summary>
/// average height of cable span between two poles, meters
/// средняя высота пролета, м
/// </summary>
public double Height;
/// <summary>
/// power kW for switches
/// Мощность ЛЭП, кВт
/// </summary>
public int Power;
/// <summary>
/// Line vertices coordinate list
/// список координат вершин ЛЭП как линейного объекта
/// </summary>
public List<Coordinate> Coordinates { get; set; }
/// <summary>
/// assigned powerstation/pole
/// идентификатор соответсвующего конца/начала линии (столб, трансформаторная подстанция, понижающая подстанция)
/// </summary>
public int PointFromID;
/// <summary>
/// assigned powerstation/pole
/// идентификатор соответсвующего конца/начала линии (столб, трансформаторная подстанция, понижающая подстанция)
/// </summary>
public int PointToID;
/// <summary>
/// broken/not broken switch
/// сломана (true) или нет (false) линяя
/// </summary>
public bool IsBroken;
/// <summary>
/// power on switch
/// получает (true) или нет (false) линия питание
/// </summary>
public bool IsON;
private int MissingIdValue = -1;
public Powerline()
{
//default constructor body
}
/// <summary>
/// designated constructor
/// </summary>
/// <param name="coord"></param>
/// <param name="id"></param>
/// <param name="yer"></param>
/// <param name="h"></param>
/// <param name="pw"></param>
/// <param name="isbrkn"></param>
/// <param name="ison"></param>
/// <param name="toID"></param>
/// <param name="fromID"></param>
public Powerline(List<Coordinate> coordinates, int id, int year, double height, int power, int toID, int fromID)
{
this.Coordinates = coordinates;
this.Identifier = id;
this.Year = year;
this.Height = height;
this.Power = power;
this.IsBroken = false;
this.IsON = false;
this.PointFromID = fromID;
this.PointToID = toID;
if (!(this.CheckValue()))
{
throw new System.ArgumentException("Powerline object wasn't initialized correctly");
}
}
/// <summary>
/// проверка валидности полей
/// </summary>
/// <returns></returns>
public bool CheckValue()
{
bool checker = true;
if ( Identifier >= 0 && Math.Abs(Year - 1985) < 45 && Math.Abs(Height - 15)<15 && Power > 0 && Power < 1000 && PointFromID >=0 && PointToID >=0)
{
checker = true;
}
else { checker = false; }
return checker;
}
}
/// <summary>
/// powerstation/pole point class
/// класс для трансформаторных подстанций/столбов/понижающих(конечных) подстанций
/// </summary>
public class PowerStation
{
/// <summary>
/// unique id
/// уникальный идентификатор подстанции/столба
/// </summary>
public int Identifier;
/// <summary>
/// Coordinates
/// </summary>
public Coordinate Coordinate;
/// <summary>
/// station name field
/// название подстанции
/// </summary>
public string Name;
/// <summary>
/// power, kW
/// мощность, кВт
/// </summary>
public int Power;
/// <summary>
/// type of point - trans/pole/endpoint
/// тип станции - трансформаторная подстанция/столб/понижающая(конечная)подстанция
/// </summary>
public enum StationType
{
Pole, Trans, Endstat
};
/// <summary>
/// тип подстанции
/// </summary>
public StationType Stationtype;
/// <summary>
/// is point a source?
/// является ли подстаниция источником питания (в случае ТЭЦ или питания от внешних для рассматриваемой цепи ЛЭП)
/// </summary>
public bool IsSource;
/// <summary>
/// power on switch
/// поступает (true) или нет (false) на подстанцию питание
/// </summary>
public bool IsON;
/// <summary>
/// asigned powerlines list
/// список оканчивающихся/начинающихся на подстанции ЛЭП
/// </summary>
public List<Powerline> LineList;
public PowerStation()
{
//default constructor
}
/// <summary>
/// designated constructor
/// </summary>
/// <param name="crds"></param>
/// <param name="id"></param>
/// <param name="stname"></param>
/// <param name="stpower"></param>
/// <param name="sttype"></param>
/// <param name="issource"></param>
/// <param name="ison"></param>
public PowerStation(Coordinate crds, int id, string stname, int stpower, StationType sttype, bool issource)
{
this.Coordinate = crds;
this.Identifier = id;
this.Name = stname;
this.Power = stpower;
this.Stationtype = sttype;
this.IsSource = issource;
this.IsON = false;
}
public bool CheckValue()
{
bool checker = true;
if (Identifier >=0 && Power >0 && Power < 1000 )
{
checker = true;
}
else { checker = false; }
return checker;
}
}
enum FunctionType
{
FunctionVelocityX = 0,
FunctionVelocityY = 1,
FunctionClimate5 = 2,
FunctionClimate10 = 3,
FunctionClimate15 = 4
}
/// <summary>
/// DTO for input
/// </summary>
public class Input
{
/// <summary>
/// prognistic raster info
/// массив прогностического ветра
/// </summary>
public List<List<PrognosticCell>> PrognosticCells { get; set; }
/// <summary>
/// prognostic raster cell info
/// параметры ячеек регулярной сетки прогностического ветра
/// </summary>
public CellSize PrognosticCellSize { get; set; }
/// <summary>
/// affine coefficients from prognostic raster projections
/// коэффициенты аффиного проеобразования из проекции массива прогностического ветра
/// </summary>
public double[] PrognosticAffineCoefficients { get; set; }
/// <summary>
/// climate raster array
/// массив климатических полей скорости ветра заданной повторяемости
/// </summary>
public List<List<ClimateCell>> ClimateCells { get; set; }
/// <summary>
/// climate raster cell info
/// параметры ячеек регулярной сетки климатических полей скорости ветра заданной повторяемости
/// </summary>
public CellSize ClimateCellSize { get; set; }
/// <summary>
/// affine coefficients from climate raster projection
/// коэффициенты аффинного преобразования из проекции массива климатических полей скорости ветра заданной повторяемости
/// </summary>
public double[] ClimateAffineCoefficients { get; set; }
/// <summary>
/// lines list
/// список ЛЭП
/// </summary>
public List<Powerline> PowerLines { get; set; }
/// <summary>
/// stations/poles list
/// список точечных объектов - трансформаторных подстанций/столбов/понижающих(конечных) подстанций
/// </summary>
public List<PowerStation> PowerStations { get; set; }
/// <summary>
/// maximum distance for line segment, meters
/// максимальное расстояние между точками ЛЭП, для которых проверяется сломаются/несломаются под действием ветра
/// </summary>
public double DistThreshold
{
get
{
return 500;
}
}
/// <summary>
/// missing raster value for cells, if you haven't normal value use this one.
/// потерянное значение для растрового объекта; Если Вы не имеете данных - используйте это значение
/// </summary>
public double kMissingRasterValue
{
get
{
return -9999;
}
}
}
/// <summary>
/// Output
/// </summary>
public class Output
{
/// <summary>
/// stations list without power
/// Список подстанций на которые не поступит питание в результате предсказанных поломок ЛЭП в сети
/// </summary>
public List<PowerStation> DisabledStations { get; set; }
/// <summary>
/// broken lines list
/// Список прогнозируемых сломанных ЛЭП в результате ветрового воздействия
/// </summary>
public List<Powerline> DisabledLines { get; set; }
}
/// <summary>
/// main calculations class
/// </summary>
public class StressPowerChecker
{
/// <summary>
/// Input Data
/// </summary>
private Input Input;
/// <summary>
/// Main function for power graph algorithm
/// Основая функция для вычисления распределения питания по графу ЛЭП-подстанции
/// </summary>
public Output CheckPower(Input input)
{
this.Input = input;
//Calculate which lines are broken
List<WindStressPRM.Powerline> prmBrokenLines = brokenPowerLinesAfterCheck();
//get the graph
PreparingPowerItems();
//start from source points
foreach (PowerStation pwstation in input.PowerStations)
{
if (pwstation.IsSource)
{
CheckPowerPointsForStation(pwstation);
}
}
//fill output
Output output = new Output();
output.DisabledStations = new List<PowerStation>();
output.DisabledLines = new List<Powerline>();
foreach (Powerline line in input.PowerLines)
{
if (line.IsBroken)
{
output.DisabledLines.Add(line);
}
}
foreach (PowerStation powerStation in input.PowerStations)
{
//stations of type pole can be disabled if the line is broken
if (!powerStation.IsON && !(powerStation.Stationtype == PowerStation.StationType.Pole))
{
output.DisabledStations.Add(powerStation);
}
}
return output;
}
/// <summary>
/// recursive search function for power graph
/// </summary>
/// <param name="sourcepoint"> powered station to search next powered station from</param>
private void CheckPowerPointsForStation(PowerStation sourcepoint)
{
if (!sourcepoint.IsON)
{
throw new Exception("CheckPowerPointsForStation is called from disabled sourcepoint");
}
// if the point is not a pole - i.e. we know
// it can redistribute power within connected lines
// we turn it ON if any of the connected lines are powered
foreach (Powerline line in sourcepoint.LineList)
{
if (!line.IsBroken && !line.IsON)
{
line.IsON = true;
foreach (PowerStation powerStation in Input.PowerStations)
{
if (powerStation.Identifier != sourcepoint.Identifier && (powerStation.Identifier == line.PointFromID || powerStation.Identifier == line.PointToID))
{
if (!(sourcepoint.Stationtype == PowerStation.StationType.Pole)) {
powerStation.IsON = true;
CheckPowerPointsForStation(powerStation);
}
else {
// if line is a pole we have to check if it's actually able to
// get electricity to other points i.e. no connected lines are broken
bool powerLineCheck = false;
foreach (Powerline powerline in powerStation.LineList)
{
if (powerline.IsBroken) {
powerLineCheck = true;
}
}
if (!powerLineCheck)
{
powerStation.IsON = true;
CheckPowerPointsForStation(powerStation);
}
}
}
}
}
//else we have broken line or already switched ON powerpoint
}
}
/// <summary>
/// preparing powerstations and lists of lines for them to get a graph-like structure
/// </summary>
private void PreparingPowerItems()
{
//First we make sure that all the sources are ON
//and all non sources are OFF
foreach (PowerStation powerStation in Input.PowerStations)
{
if (powerStation.IsSource == true) {
powerStation.IsON = true;
}
else {
powerStation.IsON = false;
}
// for each power station we create a list of powerlines it is attached to
List<Powerline> lines = new List<Powerline>();
foreach (Powerline line in Input.PowerLines)
{
//we also switch OFF all lines
line.IsON = false;
if (line.PointFromID == powerStation.Identifier || line.PointToID == powerStation.Identifier) {
lines.Add(line);
}
}
powerStation.LineList = lines;
}
}
/// <summary>
/// function to chec if powerline is broken by wind
/// </summary>
/// <returns>list of broken powerlines</returns>
private List<Powerline> brokenPowerLinesAfterCheck()
{
List<Powerline> brokenLines = new List<Powerline>();
// actually there are curves in powerLines
foreach (Powerline powerCurve in Input.PowerLines)
{
// get coordinates list
List<Coordinate> points = powerCurve.Coordinates;
List<bool> checkList = new List<bool>();
// cycle throw all points in line
for (int i = 1; i < points.Count; i++)
{
double x1 = points[i - 1].X;
double y1 = points[i - 1].Y;
double x2 = points[i].X;
double y2 = points[i].Y;
bool result = linearLineIsBroken(points[i - 1], points[i], powerCurve.Height, powerCurve.Power);
checkList.Add(result);
}
foreach (bool chkpnt in checkList)
{
if (chkpnt == true)
{
powerCurve.IsBroken = true;
}
}
if (powerCurve.IsBroken)
{
brokenLines.Add(powerCurve);
}
}
return brokenLines;
}
/// <summary>
/// Divides the line between two vertices if needed
/// depending on dist_threshold
/// and checks if line is broken by wind in any of midpoints
/// </summary>
/// <param name="coord1">first vertice coordinate</param>
/// <param name="coord2">second vertice </param>
/// <param name="heightLine">height</param>
/// <param name="power">power for climatology switches</param>
/// <returns>true if line is broken false otherwise</returns>
private bool linearLineIsBroken(Coordinate coord1, Coordinate coord2, double heightLine, int power)
{
double distance = Math.Sqrt((coord2.X - coord1.X) * (coord2.X - coord1.X) + (coord2.Y - coord1.Y) * (coord2.Y - coord1.Y));
double distpropD = distance / Input.DistThreshold;
List<Coordinate> pointlist = new List<Coordinate>();
Coordinate midpoint = new Coordinate(0, 0);
int distpropI = Convert.ToInt32(distpropD);
if (distpropI > 1)
{
double constXdiff = (coord2.X - coord1.X) / distpropI;
double constYdiff = (coord2.Y - coord1.Y) / distpropI;
Coordinate subCoord1 = new Coordinate(coord1.X, coord1.Y);
Coordinate subCoord2 = new Coordinate(coord1.X + constXdiff, coord1.Y + constXdiff);
for (int j = 0; j < distpropI; j++)
{
if (j == 0)
{
midpoint.X = (subCoord1.X + subCoord2.X) / 2;
midpoint.Y = (subCoord1.Y + subCoord2.Y) / 2;
pointlist.Add(midpoint);
}
else
{
// move to next center
subCoord1.X = subCoord2.X;
subCoord1.Y = subCoord2.Y;
subCoord2.X = subCoord1.X + constXdiff;
subCoord2.Y = subCoord1.Y + constYdiff;
midpoint.X = (subCoord1.X + subCoord2.X) / 2;
midpoint.Y = (subCoord1.Y + subCoord2.Y) / 2;
pointlist.Add(midpoint);
}
}
}
else
{
midpoint.X = (coord1.X + coord2.X) / 2;
midpoint.Y = (coord1.Y + coord2.Y) / 2;
pointlist.Add(midpoint);
}
FunctionType climateType;
if (power > 5 && power < 330)
{
climateType = FunctionType.FunctionClimate10;
}
else
{
if (power <= 5) {
climateType = FunctionType.FunctionClimate5;
}
else {
climateType = FunctionType.FunctionClimate15;
}
}
List<bool> checkbool = new List<bool>();
foreach (Coordinate coords in pointlist)
{
bool res = false;
double uwind = interpol(coords, FunctionType.FunctionVelocityX);
double vwind = interpol(coords, FunctionType.FunctionVelocityY);
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;
double sinwind = Math.Sin(anglewind);
double C_height = 1.0;
if (umod < 20)
{ //wind is too low
res = false;
}
else
{ // calculate prognostic and threshold windstress
double p1 = -0.00078501;
double p2 = 0.13431;
double p3 = -2.11112;
double p4 = 19.548;
double qpr = p1 * umod * umod * umod + p2 * umod * umod + p3 * umod + p4;
double qcl = p1 * climwind * climwind * climwind + p2 * climwind * climwind + p3 * climwind + p4;
double Ppr = qpr * C_height * sinwind * sinwind;
double Pcl = qcl * C_height * 1.0;
if (Ppr >= Pcl)
{
// here line is broken
res = true;
}
}
checkbool.Add(res);
}
bool result = false;
foreach (bool lastcheck in checkbool)
{
if (lastcheck == true) { result = true; }
}
return result;
}
private double[] affineCoefficients(FunctionType functionType)
{
switch (functionType)
{
case FunctionType.FunctionVelocityX:
case FunctionType.FunctionVelocityY:
{
return Input.PrognosticAffineCoefficients;
}
case FunctionType.FunctionClimate5:
case FunctionType.FunctionClimate10:
case FunctionType.FunctionClimate15:
{
return Input.ClimateAffineCoefficients;
}
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)
{
case FunctionType.FunctionVelocityX:
case FunctionType.FunctionVelocityY:
{
return Input.PrognosticCells[index.Row][index.Col].Coordinate;
}
case FunctionType.FunctionClimate5:
case FunctionType.FunctionClimate10:
case FunctionType.FunctionClimate15:
{
return Input.ClimateCells[index.Row][index.Col].Coordinate;
}
default:
break;
}
return null;
}
int countInList(FunctionType functionType, bool forRows)
{
switch (functionType)
{
case FunctionType.FunctionVelocityX:
case FunctionType.FunctionVelocityY:
{
return forRows ? Input.PrognosticCells.Count : Input.PrognosticCells[0].Count;
}
case FunctionType.FunctionClimate5:
case FunctionType.FunctionClimate10:
case FunctionType.FunctionClimate15:
{
return forRows ? Input.ClimateCells.Count : Input.ClimateCells[0].Count;
}
default:
break;
}
return 0;
}
private double valueForFunction(FunctionType functionType, Index index)
{
switch (functionType)
{
case FunctionType.FunctionVelocityX:
{
return Input.PrognosticCells[index.Row][index.Col].VelocityX;
}
case FunctionType.FunctionVelocityY:
{
return Input.PrognosticCells[index.Row][index.Col].VelocityY;
}
case FunctionType.FunctionClimate5:
{
return Input.ClimateCells[index.Row][index.Col].Wind5;
}
case FunctionType.FunctionClimate10:
{
return Input.ClimateCells[index.Row][index.Col].Wind10;
}
case FunctionType.FunctionClimate15:
{
return Input.ClimateCells[index.Row][index.Col].Wind15;
}
default:
break;
}
return 0;
}
private CellSize cellSizeForFunction(FunctionType functionType)
{
switch (functionType)
{
case FunctionType.FunctionVelocityX:
case FunctionType.FunctionVelocityY:
{
return Input.PrognosticCellSize;
}
case FunctionType.FunctionClimate5:
case FunctionType.FunctionClimate10:
case FunctionType.FunctionClimate15:
{
return Input.ClimateCellSize;
}
default:
break;
}
throw new Exception("There is no cell size");
}
private double interpol(Coordinate coords, FunctionType functionType)
{
// select directions for projections
const bool normalX = true;// true - East, false West
const bool normalY = false;// true - North, false South
Index rc = projectionToCell(coords, functionType);
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
Index rcBotLeft = new Index(Math.Min(row2, rc.Row), Math.Min(col2, rc.Col));
Index rcBotRight = new Index(Math.Max(row2, rc.Row), Math.Min(col2, rc.Col));
Index rcTopLeft = new Index(Math.Min(row2, rc.Row), Math.Max(col2, rc.Col));
Index rcTopRight = new 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);
Coordinate origin = cellToProjection(rcBotLeft, functionType);
bool testBotLeft = valBotLeft == Input.kMissingRasterValue;
bool testBotRight = valBotRight == Input.kMissingRasterValue;
bool testTopLeft = valTopLeft == Input.kMissingRasterValue;
bool testTopRight = valTopRight == Input.kMissingRasterValue;
if (testBotLeft || testBotRight || testTopLeft || testTopRight)
{
// tests indicates that value at test-cell is missed.
if (testTopRight && testTopLeft && testBotLeft && testBotRight)
{
throw new Exception("Current value in such a bad place, you need to fill these place with raster values");
}
int count = 0;
if (testBotLeft)
{
valBotLeft = 0;
count++;
}
if (testBotRight)
{
valBotRight = 0;
count++;
}
if (testTopLeft)
{
valTopLeft = 0;
count++;
}
if (testTopRight)
{
valTopRight = 0;
count++;
}
//of course there is count not 0;
if (count == 0) {
throw new Exception("Interpolation Logic Error");
}
double average = (valTopLeft + valTopRight + valBotLeft + valBotRight)/count;
if (testBotLeft)
{
valBotLeft = average;
}
if (testBotRight)
{
valBotRight = average;
}
if (testTopLeft)
{
valTopLeft = average;
}
if (testTopRight)
{
valTopRight = average;
}
}
// 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;
}
}
}