WSL/SLF GitLab Repository

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

The basic operators have been implemented in Grid2DObject and Grid3DObject as...

The basic operators have been implemented in Grid2DObject and Grid3DObject as well as DEMObject (with the proper update strategy). 

In order to suppress the risk of out of date "ncols, nrows" and since these were redundant with Array2D<>.nx/ny, they have been suppressed and it is therefore now mandatory to rely on the getNx/getNy/getNz getters. This is much safer but impacted quite a lot of code... 
parent 07fadabc
......@@ -22,9 +22,9 @@ int main(void) {
io.write2DGrid(dem, MeteoGrids::DEM, Date(0.));
Grid2DObject slope(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.slope);
Grid2DObject slope(dem.cellsize, dem.llcorner, dem.slope);
io.write2DGrid(slope, MeteoGrids::SLOPE, Date(0.));
Grid2DObject azi(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.azi);
Grid2DObject azi(dem.cellsize, dem.llcorner, dem.azi);
io.write2DGrid(azi,"azi.png");
return 0;
......
......@@ -17,12 +17,12 @@ int main(void) {
std::cout << "Initial air temperatures grid: " << grid2.toString() << "\n";
//simple arithmetic operations (the nodata values are preserved)
grid2.grid2D -= 273.15;
grid2 -= 273.15;
std::cout << "Air temperatures grid in celsius: " << grid2.toString() << "\n";
//operations between grids
Grid2DObject grid3(grid1.ncols, grid1.nrows, grid1.cellsize, grid1.llcorner);
grid3.grid2D = grid1.grid2D * (grid2.grid2D / 100.); //the parenthesis are required because of constness
Grid2DObject grid3(grid1.getNx(), grid1.getNy(), grid1.cellsize, grid1.llcorner);
grid3 = grid1 * (grid2 / 100.); //the parenthesis are required because of constness
//std::cout << "Grids multiplication: " << grid3 << "\n";
//locate the cell that contains a specific point in real world coordinates
......@@ -44,7 +44,7 @@ int main(void) {
//it is not possible to make grid operations if they don't have the same size
std::cout << "Q: What happens when we try to divide two grids that don't have the same size?\n\n";
try {
subgrid.grid2D /= grid1.grid2D;
subgrid /= grid1;
} catch(const std::exception &e) {
std::cout << e.what();
std::cout << "\nA: this is not possible and throws an exception!\n";
......
......@@ -413,7 +413,7 @@ const std::string BufferedIOHandler::toString() const
std::map<std::string, Grid2DObject>::const_iterator it1;
for (it1=mapBufferedGrids.begin(); it1 != mapBufferedGrids.end(); ++it1){
os << setw(10) << "Grid" << " = " << it1->first << ", ";
os << (it1->second).ncols << " x " << (it1->second).nrows << " @ " << (it1->second).cellsize << "m\n";
os << (it1->second).getNx() << " x " << (it1->second).getNy() << " @ " << (it1->second).cellsize << "m\n";
}
os << "</BufferedIOHandler>\n";
......
......@@ -277,7 +277,7 @@ double NoneAlgorithm::getQualityRating(const Date& /*i_date*/, const MeteoData::
}
void NoneAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid) {
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, IOUtils::nodata);
grid.set(dem.getNx(), dem.getNy(), dem.cellsize, dem.llcorner, IOUtils::nodata);
}
......@@ -517,8 +517,8 @@ void RHAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
retrend(dem, trend, grid);
//Recompute Rh from the interpolated td
for (size_t jj=0; jj<grid.nrows; jj++) {
for (size_t ii=0; ii<grid.ncols; ii++) {
for (size_t jj=0; jj<grid.getNy(); jj++) {
for (size_t ii=0; ii<grid.getNx(); ii++) {
double &value = grid(ii,jj);
if(value!=IOUtils::nodata)
value = Atmosphere::DewPointtoRh(value, ta(ii,jj), 1);
......@@ -572,8 +572,8 @@ void ILWRAlgorithm::calculate(const DEMObject& dem, Grid2DObject& grid)
retrend(dem, trend, grid);
//Recompute Rh from the interpolated td
for (size_t jj=0; jj<grid.nrows; jj++) {
for (size_t ii=0; ii<grid.ncols; ii++) {
for (size_t jj=0; jj<grid.getNy(); jj++) {
for (size_t ii=0; ii<grid.getNx(); ii++) {
double &value = grid(ii,jj);
value = Atmosphere::blkBody_Radiation(value, ta(ii,jj));
}
......
......@@ -33,8 +33,8 @@ const Grid2DObject ResamplingAlgorithms2D::NearestNeighbour(const Grid2DObject &
throw InvalidArgumentException(ss.str(), AT);
}
const double cellsize = i_grid.cellsize/factor;
const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.ncols)*factor ));
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor ));
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor ));
Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);
NearestNeighbour(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
......@@ -52,8 +52,8 @@ const Grid2DObject ResamplingAlgorithms2D::BilinearResampling(const Grid2DObject
throw InvalidArgumentException(ss.str(), AT);
}
const double cellsize = i_grid.cellsize/factor;
const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.ncols)*factor ));
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor ));
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor ));
Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);
Bilinear(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
......@@ -68,8 +68,8 @@ const Grid2DObject ResamplingAlgorithms2D::cubicBSplineResampling(const Grid2DOb
throw InvalidArgumentException(ss.str(), AT);
}
const double cellsize = i_grid.cellsize/factor;
const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.ncols)*factor ));
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.nrows)*factor ));
const size_t ncols = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNx())*factor ));
const size_t nrows = static_cast<size_t>(Optim::round( static_cast<double>(i_grid.getNy())*factor ));
Grid2DObject o_grid(ncols, nrows, cellsize, i_grid.llcorner);
cubicBSpline(o_grid.grid2D, i_grid.grid2D); //GridObjects always keep nodata
......
......@@ -378,7 +378,7 @@ template<class T> void Array1D<T>::abs() {
}
template<class T> const Array1D<T> Array1D<T>::getAbs() const {
Array1D<T> result = *this; //make a copy
Array1D<T> result(*this); //make a copy
result.abs(); //already implemented
return result;
......@@ -442,7 +442,7 @@ template<class T> Array1D<T>& Array1D<T>::operator+=(const Array1D<T>& rhs)
template<class T> const Array1D<T> Array1D<T>::operator+(const Array1D<T>& rhs)
{
Array1D<T> result = *this; //make a copy
Array1D<T> result(*this); //make a copy
result += rhs; //already implemented
return result;
......@@ -467,7 +467,7 @@ template<class T> Array1D<T>& Array1D<T>::operator+=(const T& rhs)
template<class T> const Array1D<T> Array1D<T>::operator+(const T& rhs)
{
Array1D<T> result = *this;
Array1D<T> result(*this);
result += rhs; //already implemented
return result;
......@@ -502,7 +502,7 @@ template<class T> Array1D<T>& Array1D<T>::operator-=(const Array1D<T>& rhs)
template<class T> const Array1D<T> Array1D<T>::operator-(const Array1D<T>& rhs)
{
Array1D<T> result = *this; //make a copy
Array1D<T> result(*this); //make a copy
result -= rhs; //already implemented
return result;
......@@ -516,7 +516,7 @@ template<class T> Array1D<T>& Array1D<T>::operator-=(const T& rhs)
template<class T> const Array1D<T> Array1D<T>::operator-(const T& rhs)
{
Array1D<T> result = *this;
Array1D<T> result(*this);
result += -rhs; //already implemented
return result;
......@@ -550,7 +550,7 @@ template<class T> Array1D<T>& Array1D<T>::operator*=(const Array1D<T>& rhs)
template<class T> const Array1D<T> Array1D<T>::operator*(const Array1D<T>& rhs)
{
Array1D<T> result = *this; //make a copy
Array1D<T> result(*this); //make a copy
result *= rhs; //already implemented
return result;
......@@ -575,7 +575,7 @@ template<class T> Array1D<T>& Array1D<T>::operator*=(const T& rhs)
template<class T> const Array1D<T> Array1D<T>::operator*(const T& rhs)
{
Array1D<T> result = *this;
Array1D<T> result(*this);
result *= rhs; //already implemented
return result;
......@@ -609,7 +609,7 @@ template<class T> Array1D<T>& Array1D<T>::operator/=(const Array1D<T>& rhs)
template<class T> const Array1D<T> Array1D<T>::operator/(const Array1D<T>& rhs)
{
Array1D<T> result = *this; //make a copy
Array1D<T> result(*this); //make a copy
result /= rhs; //already implemented
return result;
......@@ -623,14 +623,14 @@ template<class T> Array1D<T>& Array1D<T>::operator/=(const T& rhs)
template<class T> const Array1D<T> Array1D<T>::operator/(const T& rhs)
{
Array1D<T> result = *this;
Array1D<T> result(*this);
result *= 1./rhs; //already implemented
return result;
}
template<class T> bool Array1D<T>::operator==(const Array1D<T>& in) const {
size_t in_nx=in.getNx();
const size_t in_nx = in.getNx();
if(nx!=in_nx)
return false;
......
......@@ -533,7 +533,7 @@ template<class T> void Array2D<T>::abs() {
}
template<class T> const Array2D<T> Array2D<T>::getAbs() const {
Array2D<T> result = *this; //make a copy
Array2D<T> result(*this); //make a copy
result.abs(); //already implemented
return result;
......@@ -599,7 +599,7 @@ template<class T> Array2D<T>& Array2D<T>::operator+=(const Array2D<T>& rhs)
template<class T> const Array2D<T> Array2D<T>::operator+(const Array2D<T>& rhs)
{
Array2D<T> result = *this; //make a copy
Array2D<T> result(*this); //make a copy
result += rhs; //already implemented
return result;
......@@ -625,7 +625,7 @@ template<class T> Array2D<T>& Array2D<T>::operator+=(const T& rhs)
template<class T> const Array2D<T> Array2D<T>::operator+(const T& rhs)
{
Array2D<T> result = *this;
Array2D<T> result(*this);
result += rhs; //already implemented
return result;
......@@ -660,7 +660,7 @@ template<class T> Array2D<T>& Array2D<T>::operator-=(const Array2D<T>& rhs)
template<class T> const Array2D<T> Array2D<T>::operator-(const Array2D<T>& rhs)
{
Array2D<T> result = *this; //make a copy
Array2D<T> result(*this); //make a copy
result -= rhs; //already implemented
return result;
......@@ -674,7 +674,7 @@ template<class T> Array2D<T>& Array2D<T>::operator-=(const T& rhs)
template<class T> const Array2D<T> Array2D<T>::operator-(const T& rhs)
{
Array2D<T> result = *this;
Array2D<T> result(*this);
result += -rhs; //already implemented
return result;
......@@ -709,7 +709,7 @@ template<class T> Array2D<T>& Array2D<T>::operator*=(const Array2D<T>& rhs)
template<class T> const Array2D<T> Array2D<T>::operator*(const Array2D<T>& rhs)
{
Array2D<T> result = *this; //make a copy
Array2D<T> result(*this); //make a copy
result *= rhs; //already implemented
return result;
......@@ -735,7 +735,7 @@ template<class T> Array2D<T>& Array2D<T>::operator*=(const T& rhs)
template<class T> const Array2D<T> Array2D<T>::operator*(const T& rhs)
{
Array2D<T> result = *this;
Array2D<T> result(*this);
result *= rhs; //already implemented
return result;
......@@ -770,7 +770,7 @@ template<class T> Array2D<T>& Array2D<T>::operator/=(const Array2D<T>& rhs)
template<class T> const Array2D<T> Array2D<T>::operator/(const Array2D<T>& rhs)
{
Array2D<T> result = *this; //make a copy
Array2D<T> result(*this); //make a copy
result /= rhs; //already implemented
return result;
......@@ -784,14 +784,14 @@ template<class T> Array2D<T>& Array2D<T>::operator/=(const T& rhs)
template<class T> const Array2D<T> Array2D<T>::operator/(const T& rhs)
{
Array2D<T> result = *this;
Array2D<T> result(*this);
result *= (1./rhs); //already implemented
return result;
}
template<class T> bool Array2D<T>::operator==(const Array2D<T>& in) const {
size_t in_nx=in.getNx(), in_ny=in.getNy();
const size_t in_nx=in.getNx(), in_ny=in.getNy();
if(nx!=in_nx || ny!=in_ny)
return false;
......
......@@ -605,7 +605,7 @@ template<class T> void Array3D<T>::abs() {
}
template<class T> const Array3D<T> Array3D<T>::getAbs() const {
Array3D<T> result = *this; //make a copy
Array3D<T> result(*this); //make a copy
result.abs(); //already implemented
return result;
......@@ -672,7 +672,7 @@ template<class T> Array3D<T>& Array3D<T>::operator+=(const Array3D<T>& rhs)
template<class T> const Array3D<T> Array3D<T>::operator+(const Array3D<T>& rhs)
{
Array3D<T> result = *this; //make a copy
Array3D<T> result(*this); //make a copy
result += rhs; //already implemented
return result;
......@@ -698,7 +698,7 @@ template<class T> Array3D<T>& Array3D<T>::operator+=(const T& rhs)
template<class T> const Array3D<T> Array3D<T>::operator+(const T& rhs)
{
Array3D<T> result = *this;
Array3D<T> result(*this);
result += rhs; //already implemented
return result;
......@@ -733,7 +733,7 @@ template<class T> Array3D<T>& Array3D<T>::operator-=(const Array3D<T>& rhs)
template<class T> const Array3D<T> Array3D<T>::operator-(const Array3D<T>& rhs)
{
Array3D<T> result = *this; //make a copy
Array3D<T> result(*this); //make a copy
result -= rhs; //already implemented
return result;
......@@ -747,7 +747,7 @@ template<class T> Array3D<T>& Array3D<T>::operator-=(const T& rhs)
template<class T> const Array3D<T> Array3D<T>::operator-(const T& rhs)
{
Array3D<T> result = *this;
Array3D<T> result(*this);
result += -rhs; //already implemented
return result;
......@@ -782,7 +782,7 @@ template<class T> Array3D<T>& Array3D<T>::operator*=(const Array3D<T>& rhs)
template<class T> const Array3D<T> Array3D<T>::operator*(const Array3D<T>& rhs)
{
Array3D<T> result = *this; //make a copy
Array3D<T> result(*this); //make a copy
result *= rhs; //already implemented
return result;
......@@ -808,7 +808,7 @@ template<class T> Array3D<T>& Array3D<T>::operator*=(const T& rhs)
template<class T> const Array3D<T> Array3D<T>::operator*(const T& rhs)
{
Array3D<T> result = *this;
Array3D<T> result(*this);
result *= rhs; //already implemented
return result;
......@@ -843,7 +843,7 @@ template<class T> Array3D<T>& Array3D<T>::operator/=(const Array3D<T>& rhs)
template<class T> const Array3D<T> Array3D<T>::operator/(const Array3D<T>& rhs)
{
Array3D<T> result = *this; //make a copy
Array3D<T> result(*this); //make a copy
result /= rhs; //already implemented
return result;
......@@ -857,14 +857,14 @@ template<class T> Array3D<T>& Array3D<T>::operator/=(const T& rhs)
template<class T> const Array3D<T> Array3D<T>::operator/(const T& rhs)
{
Array3D<T> result = *this;
Array3D<T> result(*this);
result *= (1./rhs); //already implemented
return result;
}
template<class T> bool Array3D<T>::operator==(const Array3D<T>& in) const {
size_t in_nx=in.getNx(), in_ny=in.getNy(), in_nz=in.getNz();
const size_t in_nx=in.getNx(), in_ny=in.getNy(), in_nz=in.getNz();
if(nx!=in_nx || ny!=in_ny || nz!=in_nz)
return false;
......
......@@ -544,7 +544,7 @@ template<class T> void Array4D<T>::abs() {
template<class T> const Array4D<T> Array4D<T>::getAbs() const {
Array4D<T> result = *this; //make a copy
Array4D<T> result(*this); //make a copy
result.abs(); //already implemented
return result;
......@@ -613,7 +613,7 @@ template<class T> Array4D<T>& Array4D<T>::operator+=(const Array4D<T>& rhs)
template<class T> const Array4D<T> Array4D<T>::operator+(const Array4D<T>& rhs)
{
Array4D<T> result = *this; //make a copy
Array4D<T> result(*this); //make a copy
result += rhs; //already implemented
return result;
......@@ -639,7 +639,7 @@ template<class T> Array4D<T>& Array4D<T>::operator+=(const T& rhs)
template<class T> const Array4D<T> Array4D<T>::operator+(const T& rhs)
{
Array4D<T> result = *this;
Array4D<T> result(*this);
result += rhs; //already implemented
return result;
......@@ -674,7 +674,7 @@ template<class T> Array4D<T>& Array4D<T>::operator-=(const Array4D<T>& rhs)
template<class T> const Array4D<T> Array4D<T>::operator-(const Array4D<T>& rhs)
{
Array4D<T> result = *this; //make a copy
Array4D<T> result(*this); //make a copy
result -= rhs; //already implemented
return result;
......@@ -688,7 +688,7 @@ template<class T> Array4D<T>& Array4D<T>::operator-=(const T& rhs)
template<class T> const Array4D<T> Array4D<T>::operator-(const T& rhs)
{
Array4D<T> result = *this;
Array4D<T> result(*this);
result += -rhs; //already implemented
return result;
......@@ -723,7 +723,7 @@ template<class T> Array4D<T>& Array4D<T>::operator*=(const Array4D<T>& rhs)
template<class T> const Array4D<T> Array4D<T>::operator*(const Array4D<T>& rhs)
{
Array4D<T> result = *this; //make a copy
Array4D<T> result(*this); //make a copy
result *= rhs; //already implemented
return result;
......@@ -749,7 +749,7 @@ template<class T> Array4D<T>& Array4D<T>::operator*=(const T& rhs)
template<class T> const Array4D<T> Array4D<T>::operator*(const T& rhs)
{
Array4D<T> result = *this;
Array4D<T> result(*this);
result *= rhs; //already implemented
return result;
......@@ -784,7 +784,7 @@ template<class T> Array4D<T>& Array4D<T>::operator/=(const Array4D<T>& rhs)
template<class T> const Array4D<T> Array4D<T>::operator/(const Array4D<T>& rhs)
{
Array4D<T> result = *this; //make a copy
Array4D<T> result(*this); //make a copy
result /= rhs; //already implemented
return result;
......@@ -798,14 +798,14 @@ template<class T> Array4D<T>& Array4D<T>::operator/=(const T& rhs)
template<class T> const Array4D<T> Array4D<T>::operator/(const T& rhs)
{
Array4D<T> result = *this;
Array4D<T> result(*this);
result *= (1./rhs); //already implemented
return result;
}
template<class T> bool Array4D<T>::operator==(const Array4D<T>& in) const {
size_t in_nx=in.getNx(), in_ny=in.getNy(), in_nz=in.getNz(), in_nw=in.getNw();
const size_t in_nx=in.getNx(), in_ny=in.getNy(), in_nz=in.getNz(), in_nw=in.getNw();
if(nx!=in_nx || ny!=in_ny || nz!=in_nz || nw!=in_nw)
return false;
......
......@@ -70,18 +70,15 @@ DEMObject::DEMObject(const size_t& i_ncols, const size_t& i_nrows,
/**
* @brief Constructor that sets variables.
* @param i_ncols number of colums in the grid2D
* @param i_nrows number of rows in the grid2D
* @param i_cellsize value for cellsize in grid2D
* @param i_llcorner lower lower corner point
* @param i_altitude grid2D of elevations
* @param i_update also update slope/normals/curvatures and their min/max? (default=true)
* @param i_algorithm specify the default algorithm to use for slope computation (default=DFLT)
*/
DEMObject::DEMObject(const size_t& i_ncols, const size_t& i_nrows,
const double& i_cellsize, const Coords& i_llcorner, const Array2D<double>& i_altitude,
DEMObject::DEMObject(const double& i_cellsize, const Coords& i_llcorner, const Array2D<double>& i_altitude,
const bool& i_update, const slope_type& i_algorithm)
: Grid2DObject(i_ncols, i_nrows, i_cellsize, i_llcorner, i_altitude),
: Grid2DObject(i_cellsize, i_llcorner, i_altitude),
slope(), azi(), curvature(), Nx(), Ny(), Nz(),
min_altitude(Cst::dbl_max), min_slope(Cst::dbl_max), min_curvature(Cst::dbl_max),
max_altitude(Cst::dbl_min), max_slope(Cst::dbl_min), max_curvature(Cst::dbl_min),
......@@ -104,7 +101,7 @@ DEMObject::DEMObject(const size_t& i_ncols, const size_t& i_nrows,
* @param i_algorithm specify the default algorithm to use for slope computation (default=DFLT)
*/
DEMObject::DEMObject(const Grid2DObject& i_dem, const bool& i_update, const slope_type& i_algorithm)
: Grid2DObject(i_dem.ncols, i_dem.nrows, i_dem.cellsize, i_dem.llcorner, i_dem.grid2D),
: Grid2DObject(i_dem.cellsize, i_dem.llcorner, i_dem.grid2D),
slope(), azi(), curvature(), Nx(), Ny(), Nz(),
min_altitude(Cst::dbl_max), min_slope(Cst::dbl_max), min_curvature(Cst::dbl_max),
max_altitude(Cst::dbl_min), max_slope(Cst::dbl_min), max_curvature(Cst::dbl_min),
......@@ -216,16 +213,16 @@ void DEMObject::update(const slope_type& algorithm) {
// Creating tables
if(update_flag&SLOPE) {
slope.resize(ncols, nrows);
azi.resize(ncols, nrows);
slope.resize(getNx(), getNy());
azi.resize(getNx(), getNy());
}
if(update_flag&CURVATURE) {
curvature.resize(ncols, nrows);
curvature.resize(getNx(), getNy());
}
if(update_flag&NORMAL) {
Nx.resize(ncols, nrows);
Ny.resize(ncols, nrows);
Nz.resize(ncols, nrows);
Nx.resize(getNx(), getNy());
Ny.resize(getNx(), getNy());
Nz.resize(getNx(), getNy());
}
CalculateAziSlopeCurve(algorithm);
......@@ -316,6 +313,8 @@ void DEMObject::updateAllMinMax() {
*/
void DEMObject::printFailures() {
bool header=true;
const size_t ncols = getNx();
const size_t nrows = getNy();
if(update_flag&SLOPE) {
for ( size_t j = 0; j < nrows; j++ ) {
......@@ -362,6 +361,9 @@ void DEMObject::printFailures() {
*/
void DEMObject::sanitize() {
if(slope_failures>0 || curvature_failures>0) {
const size_t ncols = getNx();
const size_t nrows = getNy();
for ( size_t j = 0; j < nrows; j++ ) {
for ( size_t i = 0; i < ncols; i++ ) {
if(update_flag&SLOPE) {
......@@ -389,13 +391,15 @@ void DEMObject::sanitize() {
*/
Grid2DObject DEMObject::getHillshade(const double& elev, const double& azimuth) const
{
Grid2DObject hillshade(ncols, nrows, cellsize, llcorner);
if(slope.isEmpty() || azi.isEmpty())
throw InvalidArgumentException("Hillshade computation requires slope and azimuth!", AT);
const double zenith_rad = (90.-elev)*Cst::to_rad;
const double azimuth_rad = azimuth*Cst::to_rad;
const size_t ncols = getNx();
const size_t nrows = getNy();
Grid2DObject hillshade(ncols, nrows, cellsize, llcorner);
for ( size_t j = 0; j < nrows; j++ ) {
for ( size_t i = 0; i < ncols; i++ ) {
......@@ -547,7 +551,7 @@ void DEMObject::getPointsBetween(Coords point1, Coords point2, std::vector<GRID_
pts.ix = ix;
pts.iy = iy;
//make sure we only return points within the dem
if(ix>0 && ix<(signed)ncols && iy>0 && iy<(signed)nrows) {
if(ix>0 && ix<(signed)getNx() && iy>0 && iy<(signed)getNy()) {
vec_points.push_back(pts);
}
}
......@@ -577,17 +581,17 @@ void DEMObject::getPointsBetween(const Coords& point, const double& bearing, std
//define the boundaries according to the quadrant we are in
double xlim, ylim;
if(bear>=0. && bear<90.) {
xlim = (double)(ncols-1);
ylim = (double)(nrows-1);
xlim = (double)(getNx()-1);
ylim = (double)(getNy()-1);
} else if (bear>=90. && bear<180.) {
xlim = (double)(ncols-1);
xlim = (double)(getNx()-1);
ylim = 0.;
} else if (bear>=180. && bear<270.) {
xlim = 0.;
ylim = 0.;
} else {
xlim = 0.;
ylim = (double)(nrows-1);
ylim = (double)(getNy()-1);
}
//calculate the two possible intersections between the bearing line and the boundaries
......@@ -691,8 +695,8 @@ void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
}
//Now, calculate the parameters using the previously defined function pointer
for ( size_t j = 0; j < nrows; j++ ) {
for ( size_t i = 0; i < ncols; i++ ) {
for ( size_t j = 0; j < getNy(); j++ ) {
for ( size_t i = 0; i < getNx(); i++ ) {
if( grid2D(i,j) == IOUtils::nodata ) {
if(update_flag&SLOPE) {
slope(i,j) = azi(i,j) = IOUtils::nodata;
......@@ -726,8 +730,8 @@ void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
}
if((update_flag&SLOPE) && (algorithm==D8)) { //extra processing required: discretization