Skip to content
Snippets Groups Projects
Commit aafb55c2 authored by Антон Кудряшов's avatar Антон Кудряшов
Browse files

- removed coordinates from cells

parent a09fbd53
Branches
No related tags found
No related merge requests found
...@@ -91,7 +91,7 @@ namespace MES_Wind ...@@ -91,7 +91,7 @@ namespace MES_Wind
private void btnCalcStress_Click(object sender, EventArgs e) private void btnCalcStress_Click(object sender, EventArgs e)
{ {
const double eps = 0.0001; const double eps = 2;
const double RasterMissingValue = -9999; const double RasterMissingValue = -9999;
//extract prognostic u layer //extract prognostic u layer
IMapRasterLayer uRasterLayer = default(IMapRasterLayer); IMapRasterLayer uRasterLayer = default(IMapRasterLayer);
...@@ -130,16 +130,19 @@ namespace MES_Wind ...@@ -130,16 +130,19 @@ namespace MES_Wind
int rcountPrognostic = uRasterLayer.DataSet.NumRows; int rcountPrognostic = uRasterLayer.DataSet.NumRows;
int ccountPrognostic = uRasterLayer.DataSet.NumColumns; int ccountPrognostic = uRasterLayer.DataSet.NumColumns;
WindStressPRM.Matrix<WindStressPRM.PrognosticCell> prognosticMatrix = new WindStressPRM.Matrix<WindStressPRM.PrognosticCell>(); 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.Size = new WindStressPRM.CellSize(uRasterLayer.DataSet.CellWidth, uRasterLayer.DataSet.CellHeight);
prognosticMatrix.AffineCoefficients = uRasterLayer.Bounds.AffineCoefficients; prognosticMatrix.AffineCoefficients = uRasterLayer.Bounds.AffineCoefficients;
// fill cells of prognostic matrix // fill cells of prognostic matrix
for (int i = 0; i < rcountPrognostic; i++) for (int i = 0; i < rcountPrognostic; i++)
{ {
for (int j = 0; j < ccountPrognostic; j++) for (int j = 0; j < ccountPrognostic; j++)
{
if (i == 0 && j == 0)
{ {
Coordinate dummyRCoords = uRasterLayer.Bounds.CellCenter_ToProj(i, j); Coordinate dummyRCoords = uRasterLayer.Bounds.CellCenter_ToProj(i, j);
WindStressPRM.Coordinate cellCoords =new WindStressPRM.Coordinate(dummyRCoords.X,dummyRCoords.Y); prognosticMatrix.Origin = new WindStressPRM.Coordinate(dummyRCoords.X, dummyRCoords.Y);
}
double uValue = uRasterLayer.DataSet.Value[i, j]; double uValue = uRasterLayer.DataSet.Value[i, j];
double vValue = vRasterLayer.DataSet.Value[i, j]; double vValue = vRasterLayer.DataSet.Value[i, j];
if (Math.Abs(uValue - RasterMissingValue) < eps) { if (Math.Abs(uValue - RasterMissingValue) < eps) {
...@@ -148,15 +151,15 @@ namespace MES_Wind ...@@ -148,15 +151,15 @@ namespace MES_Wind
if (Math.Abs(vValue - RasterMissingValue) < eps) { if (Math.Abs(vValue - RasterMissingValue) < eps) {
vValue = Double.NaN; vValue = Double.NaN;
} }
WindStressPRM.PrognosticCell dummyPrognosticCell = new WindStressPRM.PrognosticCell(cellCoords, uValue, vValue); WindStressPRM.PrognosticCell dummyPrognosticCell = new WindStressPRM.PrognosticCell(uValue, vValue);
prognosticMatrix.Cells[i, j] = dummyPrognosticCell; prognosticMatrix.Cells[j, i] = dummyPrognosticCell;
} }
} }
//Now we create climate raster class //Now we create climate raster class
int rowCountClim = clim5RasterLayer.DataSet.NumRows; int rowCountClim = clim5RasterLayer.DataSet.NumRows;
int columnCountClim = clim5RasterLayer.DataSet.NumColumns; int columnCountClim = clim5RasterLayer.DataSet.NumColumns;
WindStressPRM.Matrix<WindStressPRM.ClimateCell> climateMatrix = new WindStressPRM.Matrix<WindStressPRM.ClimateCell>(); 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.Size = new WindStressPRM.CellSize(clim5RasterLayer.DataSet.CellWidth, clim5RasterLayer.DataSet.CellHeight);
climateMatrix.AffineCoefficients = clim5RasterLayer.Bounds.AffineCoefficients; climateMatrix.AffineCoefficients = clim5RasterLayer.Bounds.AffineCoefficients;
// fill cells of climate matrix // fill cells of climate matrix
...@@ -164,8 +167,10 @@ namespace MES_Wind ...@@ -164,8 +167,10 @@ namespace MES_Wind
{ {
for (int j = 0; j < columnCountClim; j++) for (int j = 0; j < columnCountClim; j++)
{ {
if (i == 0 && j == 0) {
Coordinate dummyCellCoords = clim15RasterLayer.Bounds.CellCenter_ToProj(i, j); Coordinate dummyCellCoords = clim15RasterLayer.Bounds.CellCenter_ToProj(i, j);
WindStressPRM.Coordinate dummyClimCoords = new WindStressPRM.Coordinate(dummyCellCoords.X,dummyCellCoords.Y); climateMatrix.Origin = new WindStressPRM.Coordinate(dummyCellCoords.X, dummyCellCoords.Y);;
}
//add or substruct in range 0 - 27 to change what lines will be broken //add or substruct in range 0 - 27 to change what lines will be broken
double clim5 = clim5RasterLayer.DataSet.Value[i, j]; double clim5 = clim5RasterLayer.DataSet.Value[i, j];
double clim10 = clim10RasterLayer.DataSet.Value[i, j] +1; double clim10 = clim10RasterLayer.DataSet.Value[i, j] +1;
...@@ -182,8 +187,8 @@ namespace MES_Wind ...@@ -182,8 +187,8 @@ namespace MES_Wind
{ {
clim15 = Double.NaN; clim15 = Double.NaN;
} }
WindStressPRM.ClimateCell dummyClim = new WindStressPRM.ClimateCell(dummyClimCoords, clim5, clim10, clim15); WindStressPRM.ClimateCell dummyClim = new WindStressPRM.ClimateCell(clim5, clim10, clim15);
climateMatrix.Cells[i, j] = dummyClim; climateMatrix.Cells[j, i] = dummyClim;
} }
} }
// create PRM_line list to pass to PRM_wind from loaded line layer // create PRM_line list to pass to PRM_wind from loaded line layer
...@@ -265,9 +270,8 @@ namespace MES_Wind ...@@ -265,9 +270,8 @@ namespace MES_Wind
{ {
Coordinate coords = new Coordinate(disabledStation.Coordinate.X, disabledStation.Coordinate.Y); Coordinate coords = new Coordinate(disabledStation.Coordinate.X, disabledStation.Coordinate.Y);
DotSpatial.Topology.Point disabledPoint = new DotSpatial.Topology.Point(coords); DotSpatial.Topology.Point disabledPoint = new DotSpatial.Topology.Point(coords);
Feature disabledPointFeature = new Feature(disabledPoint); IFeature disabledPointFeature = disabledPointSet.AddFeature(disabledPoint);
disabledPointFeature.DataRow[stypecolumn] = disabledStation.Stationtype; disabledPointFeature.DataRow.SetField(stypecolumn, disabledStation.Stationtype);
disabledPointSet.AddFeature(disabledPointFeature);
} }
//add result layer to the map //add result layer to the map
IMapPointLayer resultPointLayer = (MapPointLayer)map1.Layers.Add(disabledPointSet); IMapPointLayer resultPointLayer = (MapPointLayer)map1.Layers.Add(disabledPointSet);
......
...@@ -27,20 +27,14 @@ namespace WindStressPRM ...@@ -27,20 +27,14 @@ namespace WindStressPRM
/// </summary> /// </summary>
public double Wind15 { get; private set; } public double Wind15 { get; private set; }
/// <summary> /// <summary>
/// Cell center coordinate pair
/// Координаты центра ячейки
/// </summary>
public Coordinate Coordinate { get; private set; }
/// <summary>
/// designated constructor /// designated constructor
/// </summary> /// </summary>
/// <param name="coord"></param> /// <param name="coord"></param>
/// <param name="w5"></param> /// <param name="w5"></param>
/// <param name="w10"></param> /// <param name="w10"></param>
/// <param name="w15"></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.Wind5 = w5;
this.Wind10 = w10; this.Wind10 = w10;
this.Wind15 = w15; this.Wind15 = w15;
......
...@@ -22,19 +22,13 @@ namespace WindStressPRM ...@@ -22,19 +22,13 @@ namespace WindStressPRM
/// </summary> /// </summary>
public double VelocityY { get; private set; } public double VelocityY { get; private set; }
/// <summary> /// <summary>
/// Cell center coordinates
/// Координаты центра ячейки
/// </summary>
public Coordinate Coordinate { get; private set; }
/// <summary>
/// designated constructor /// designated constructor
/// </summary> /// </summary>
/// <param name="coord"> Coordinate pair</param> /// <param name="coord"> Coordinate pair</param>
/// <param name="vX"> U component</param> /// <param name="vX"> U component</param>
/// <param name="vY"> V 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.VelocityX = vX;
this.VelocityY = vY; this.VelocityY = vY;
if (!(this.CheckValue())) if (!(this.CheckValue()))
......
...@@ -18,10 +18,28 @@ namespace WindStressPRM ...@@ -18,10 +18,28 @@ namespace WindStressPRM
/// </summary> /// </summary>
/// <param name="index"></param> /// <param name="index"></param>
/// <returns></returns> /// <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); return new Coordinate(Origin.X + index.Row*Size.Width, Origin.Y - index.Col*Size.Height);
} }
/// <summary> /// <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 /// AffineCoefficients for matrix from raster projection
/// </summary> /// </summary>
public double[] AffineCoefficients { get; set; } public double[] AffineCoefficients { get; set; }
......
...@@ -278,62 +278,6 @@ namespace WindStressPRM ...@@ -278,62 +278,6 @@ namespace WindStressPRM
return result; 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) private Coordinate cellToProjection(Index index, FunctionType functionType)
{ {
switch (functionType) switch (functionType)
...@@ -341,13 +285,13 @@ namespace WindStressPRM ...@@ -341,13 +285,13 @@ namespace WindStressPRM
case FunctionType.FunctionVelocityX: case FunctionType.FunctionVelocityX:
case FunctionType.FunctionVelocityY: case FunctionType.FunctionVelocityY:
{ {
return Input.PrognosticCells.Cells[index.Row, index.Col].Coordinate; return Input.PrognosticCells.CellToProjection(index);
} }
case FunctionType.FunctionClimate5: case FunctionType.FunctionClimate5:
case FunctionType.FunctionClimate10: case FunctionType.FunctionClimate10:
case FunctionType.FunctionClimate15: case FunctionType.FunctionClimate15:
{ {
return Input.ClimateCells.Cells[index.Row, index.Col].Coordinate; return Input.ClimateCells.CellToProjection(index);
} }
default: default:
break; break;
...@@ -427,6 +371,26 @@ namespace WindStressPRM ...@@ -427,6 +371,26 @@ namespace WindStressPRM
throw new Exception("There is no cell size"); 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) private double interpol(Coordinate coords, FunctionType functionType)
{ {
// select directions for projections // select directions for projections
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment