WSL/SLF GitLab Repository

Commit fd9b0d8e authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A new constructor has been implemented for the GridObjects as well as better...

A new constructor has been implemented for the GridObjects as well as better members alignement. A potential issue with Winstral has been fixed. A wind direction algorithm has been implemented according to Ryan, 1977.
parent 907a717e
......@@ -37,23 +37,30 @@ Grid2DObject& Grid2DObject::operator=(const Grid2DObject& source) {
* Default constructor.
* grid2D attribute is initialized by Array2D default constructor.
*/
Grid2DObject::Grid2DObject() : grid2D(), ncols(0), nrows(0), cellsize(0), llcorner() {}
Grid2DObject::Grid2DObject() : grid2D(), llcorner(), cellsize(0.), ncols(0), nrows(0) {}
Grid2DObject::Grid2DObject(const size_t& i_ncols, const size_t& i_nrows,
const double& i_cellsize, const Coords& i_llcorner) : grid2D(i_ncols, i_nrows, IOUtils::nodata), ncols(i_ncols), nrows(i_nrows), cellsize(i_cellsize), llcorner(i_llcorner) {}
const double& i_cellsize, const Coords& i_llcorner)
: grid2D(i_ncols, i_nrows, IOUtils::nodata), llcorner(i_llcorner),
cellsize(i_cellsize), ncols(i_ncols), nrows(i_nrows) {}
Grid2DObject::Grid2DObject(const size_t& i_ncols, const size_t& i_nrows,
const double& i_cellsize, const Coords& i_llcorner, const Array2D<double>& i_grid2D) : grid2D(i_grid2D), ncols(i_ncols), nrows(i_nrows), cellsize(i_cellsize), llcorner(i_llcorner) {}
const double& i_cellsize, const Coords& i_llcorner, const Array2D<double>& i_grid2D)
: grid2D(i_grid2D), llcorner(i_llcorner), cellsize(i_cellsize), ncols(i_ncols), nrows(i_nrows) {}
Grid2DObject::Grid2DObject(const size_t& i_ncols, const size_t& i_nrows,
const double& i_cellsize, const Coords& i_llcorner, const double& init) : grid2D(i_ncols, i_nrows, init), ncols(i_ncols), nrows(i_nrows), cellsize(i_cellsize), llcorner(i_llcorner) {}
const double& i_cellsize, const Coords& i_llcorner, const double& init)
: grid2D(i_ncols, i_nrows, init), llcorner(i_llcorner), cellsize(i_cellsize), ncols(i_ncols), nrows(i_nrows) {}
Grid2DObject::Grid2DObject(const Grid2DObject& i_grid, const double& init)
: grid2D(i_grid.ncols, i_grid.nrows, init), llcorner(i_grid.llcorner), cellsize(i_grid.cellsize),
ncols(i_grid.ncols), nrows(i_grid.nrows) {}
Grid2DObject::Grid2DObject(const Grid2DObject& i_grid2Dobj, const size_t& i_nx, const size_t& i_ny,
const size_t& i_ncols, const size_t& i_nrows)
: grid2D(i_grid2Dobj.grid2D, i_nx,i_ny, i_ncols,i_nrows), ncols(i_ncols), nrows(i_nrows), cellsize(i_grid2Dobj.cellsize), llcorner(i_grid2Dobj.llcorner)
: grid2D(i_grid2Dobj.grid2D, i_nx,i_ny, i_ncols,i_nrows), llcorner(i_grid2Dobj.llcorner), cellsize(i_grid2Dobj.cellsize),
ncols(i_ncols), nrows(i_nrows)
{
//setValues(i_ncols, i_nrows, i_grid2Dobj.cellsize);
//we take the previous corner (so we use the same projection parameters)
//and we shift it by the correct X and Y distance
//llcorner = i_grid2Dobj.llcorner;
......
......@@ -80,6 +80,8 @@ class Grid2DObject {
Grid2DObject(const size_t& ncols, const size_t& nrows,
const double& cellsize, const Coords& i_llcorner, const double& init);
Grid2DObject(const Grid2DObject& i_grid, const double& init);
Grid2DObject(const size_t& ncols, const size_t& nrows,
const double& cellsize, const Coords& i_llcorner, const Array2D<double>& grid2D_in);
......@@ -179,10 +181,10 @@ class Grid2DObject {
bool clusterization(const std::vector<double>& thresholds, const std::vector<double>& ids);
Array2D<double> grid2D; ///<the grid itself (simple 2D table containing the values for each point)
Coords llcorner; ///<lower left corner of the grid
double cellsize; ///<dimension in meters of a cell (considered to be square)
size_t ncols; ///<number of columns in the grid //HACK should point to grid2D.nx
size_t nrows; ///<number of rows in the grid
double cellsize; ///<dimension in meters of a cell (considered to be square)
Coords llcorner; ///<lower left corner of the grid
protected:
void setValues(const size_t& ncols, const size_t& nrows,
......
......@@ -37,51 +37,41 @@ Grid3DObject& Grid3DObject::operator=(const Grid3DObject& source) {
}
Grid3DObject::Grid3DObject()
: grid3D(), llcorner(), cellsize(0.), z(), ncols(0), nrows(0), ndepths(0), z_is_absolute(true)
{
}
: grid3D(), z(), llcorner(), cellsize(0.), ncols(0), nrows(0), ndepths(0), z_is_absolute(true) {}
Grid3DObject::Grid3DObject(const Grid3DObject& i_grid3Dobj,
Grid3DObject::Grid3DObject(const Grid3DObject& i_grid3D,
const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
const size_t& i_nwidths, const size_t& i_nheights, const size_t& i_ndepths)
: grid3D(i_grid3Dobj.grid3D, i_nx,i_ny,i_nz, i_nwidths,i_nheights,i_ndepths), llcorner(i_grid3Dobj.llcorner),
cellsize(i_grid3Dobj.cellsize), z(), ncols(i_nwidths), nrows(i_nheights), ndepths(i_ndepths), z_is_absolute(true)
: grid3D(i_grid3D.grid3D, i_nx,i_ny,i_nz, i_nwidths,i_nheights,i_ndepths), z(i_grid3D.z), llcorner(i_grid3D.llcorner),
cellsize(i_grid3D.cellsize), ncols(i_nwidths), nrows(i_nheights), ndepths(i_ndepths), z_is_absolute(true)
{
//setValues(i_nwidths, i_nheights, i_ndepths, i_grid3Dobj.cellsize);
//we take the previous corner (so we use the same projection parameters)
//and we shift it by the correct X and Y distance
if( (llcorner.getEasting()!=IOUtils::nodata) && (llcorner.getNorthing()!=IOUtils::nodata) ) {
llcorner.setXY( llcorner.getEasting()+static_cast<double>(i_nx)*i_grid3Dobj.cellsize,
llcorner.getNorthing()+static_cast<double>(i_ny)*i_grid3Dobj.cellsize,
llcorner.getAltitude()+static_cast<double>(i_nz)*i_grid3Dobj.cellsize );
llcorner.setXY( llcorner.getEasting()+static_cast<double>(i_nx)*i_grid3D.cellsize,
llcorner.getNorthing()+static_cast<double>(i_ny)*i_grid3D.cellsize,
llcorner.getAltitude()+static_cast<double>(i_nz)*i_grid3D.cellsize );
}
}
Grid3DObject::Grid3DObject(const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepths,
const double& i_cellsize, const Coords& i_llcorner)
: grid3D(i_ncols, i_nrows, i_ndepths, IOUtils::nodata), llcorner(i_llcorner),
cellsize(i_cellsize), z(), ncols(i_ncols), nrows(i_nrows), ndepths(i_ndepths), z_is_absolute(true)
{
//setValues(i_ncols, i_nrows, i_ndepths, i_cellsize, i_llcorner);
}
: grid3D(i_ncols, i_nrows, i_ndepths, IOUtils::nodata), z(), llcorner(i_llcorner),
cellsize(i_cellsize), ncols(i_ncols), nrows(i_nrows), ndepths(i_ndepths), z_is_absolute(true) {}
Grid3DObject::Grid3DObject(const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepths,
const double& i_cellsize, const Coords& i_llcorner, const double& init)
: grid3D(i_ncols, i_nrows, i_ndepths, init), llcorner(i_llcorner),
cellsize(i_cellsize), z(), ncols(i_ncols), nrows(i_nrows), ndepths(i_ndepths), z_is_absolute(true)
{
//setValues(i_ncols, i_nrows, i_ndepths, i_cellsize, i_llcorner);
}
: grid3D(i_ncols, i_nrows, i_ndepths, init), z(), llcorner(i_llcorner),
cellsize(i_cellsize), ncols(i_ncols), nrows(i_nrows), ndepths(i_ndepths), z_is_absolute(true) {}
Grid3DObject::Grid3DObject(const Grid3DObject& i_grid, const double& init)
: grid3D(i_grid.ncols, i_grid.nrows, i_grid.ndepths, init), z(i_grid.z), llcorner(i_grid.llcorner),
cellsize(i_grid.cellsize), ncols(i_grid.ncols), nrows(i_grid.nrows), ndepths(i_grid.ndepths), z_is_absolute(i_grid.z_is_absolute) {}
Grid3DObject::Grid3DObject(const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepths,
const double& i_cellsize, const Coords& i_llcorner, const Array3D<double>& i_grid3D)
: grid3D(i_grid3D), llcorner(i_llcorner),
cellsize(i_cellsize), z(), ncols(i_ncols), nrows(i_nrows), ndepths(i_ndepths), z_is_absolute(true)
{
//set(i_ncols, i_nrows, i_ndepths, i_cellsize, i_llcorner, i_grid3D);
}
: grid3D(i_grid3D), z(), llcorner(i_llcorner),
cellsize(i_cellsize), ncols(i_ncols), nrows(i_nrows), ndepths(i_ndepths), z_is_absolute(true) {}
bool Grid3DObject::gridify(std::vector<Coords>& vec_points) const
{
......
......@@ -74,9 +74,9 @@ class Grid3DObject{
* one passed as i_grid3Dobj argument. The resulting Grid3DObject is a by value copy of
* a subspace of the space spanned by the i_grid3Dobj
*/
Grid3DObject(const Grid3DObject& i_grid3Dobj,
Grid3DObject(const Grid3DObject& i_grid3D,
const size_t& i_nx, const size_t& i_ny, const size_t& i_nz,
const size_t& i_nwidths, const size_t& i_nheights, const size_t& i_ndepths);
const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepths);
Grid3DObject(const size_t& ncols, const size_t& nrows, const size_t& ndepths,
const double& cellsize, const Coords& i_llcorner);
......@@ -84,6 +84,8 @@ class Grid3DObject{
Grid3DObject(const size_t& ncols, const size_t& nrows, const size_t& ndepths,
const double& cellsize, const Coords& i_llcorner, const double& init);
Grid3DObject(const Grid3DObject& i_grid, const double& init);
Grid3DObject(const size_t& ncols, const size_t& nrows, const size_t& ndepths,
const double& cellsize, const Coords& i_llcorner, const Array3D<double>& grid3D);
......@@ -167,9 +169,9 @@ class Grid3DObject{
void extractLayer(const size_t& i_z, Grid2DObject& layer);
Array3D<double> grid3D;
std::vector<double> z; ///> Vector of depths
Coords llcorner;
double cellsize;
std::vector<double> z; ///> Vector of depths
size_t ncols, nrows, ndepths;
bool z_is_absolute; ///> Are z coordinates absolute or relative to a DEM?
......
......@@ -561,6 +561,32 @@ void Interpol2D::SteepSlopeRedistribution(const DEMObject& dem, const Grid2DObje
}
}
//Compute the wind direction changes by the terrain, see Ryan, "a mathematical model for diagnosis
//and prediction of surface winds in mountainous terrain", 1977, journal of applied meteorology, 16, 6
/**
* @brief compute the change of wind direction by the local terrain
* This is according to Ryan, <i>"a mathematical model for diagnosis and prediction of surface
* winds in mountainous terrain"</i>, 1977, journal of applied meteorology, <b>16</b>, 6.
* @param dem array of elevations (dem). The slope and azimuth must have been updated as they are required for the DEM analysis.
* @param grid 2D array of wind direction to fill
* @author Mathias Bavay
*/
void Interpol2D::RyanWindDir(const DEMObject& dem, Grid2DObject &grid)
{
for(size_t ii=0; ii<grid.getNx()*grid.getNy(); ii++) {
const double azi = dem.azi(ii);
const double slope = dem.slope(ii);
if (azi==IOUtils::nodata || slope==IOUtils::nodata) {
grid(ii) = IOUtils::nodata;
continue;
}
const double Yd = 100.*tan(slope*Cst::to_rad);
const double Fd = -0.225 * std::min(Yd, 100.) * sin(2.*(grid(ii)-azi)*Cst::to_rad);
grid(ii) += Fd;
}
}
//compute the Winstral sx factor for one single direction and one single point (ii,jj) in dem up to dmax distance
double Interpol2D::WinstralSX_core(const Grid2DObject& dem, const double& dmax, const double& bearing, const size_t& i, const size_t& j)
{
......@@ -655,8 +681,9 @@ void Interpol2D::WinstralSX(const DEMObject& dem, const double& dmax, const doub
const double bearing_inc = 5.;
const double bearing_width = 30.;
const double bearing1 = in_bearing - bearing_width/2.;
const double bearing2 = in_bearing + bearing_width/2.;
double bearing1 = fmod( in_bearing - bearing_width/2., 360. );
double bearing2 = fmod( in_bearing + bearing_width/2., 360. );
if (bearing1>bearing2) std::swap(bearing1, bearing2);
const size_t ncols = dem.ncols, nrows = dem.nrows;
for(size_t jj = 0; jj<nrows; jj++) {
......
......@@ -90,6 +90,7 @@ class Interpol2D {
static void steepestDescentDisplacement(const DEMObject& dem, const Grid2DObject& grid, const size_t& ii, const size_t& jj, short &d_i_dest, short &d_j_dest);
static double depositAroundCell(const DEMObject& dem, const size_t& ii, const size_t& jj, const double& precip, Grid2DObject &grid);
static void RyanWindDir(const DEMObject& dem, Grid2DObject &grid);
static double WinstralSX_core(const Grid2DObject& dem, const double& dmax, const double& bearing, const size_t& ii, const size_t& jj);
static double AvgSX_core(const Grid2DObject& dem, const Grid2DObject& sx, const double& dmax, const double& bearing, const size_t& ii, const size_t& jj);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment