WSL/SLF GitLab Repository

Commit 97ebe091 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

A typo had slipped through the last commit for Array4D. The Corripio slope...

A typo had slipped through the last commit for Array4D. The Corripio slope computation was till wrong, this has finally been fixed and checked. The slope default algorithm selection was not kept through BufferedIOHandler, this has been fixed. Moreover, unecessary copies of grids were being made in BufferedIOHandler and have been removed. The tests that were not passing anymore because of the slope changes have been updated, including improved error reporting in dem_reading.
parent dbc74377
......@@ -528,12 +528,12 @@ template<class T> void Array4D<T>::abs() {
if(std::numeric_limits<T>::is_signed) {
const unsigned int nwyz = nwnxny*nz;
if(keep_nodata==false) {
for (unsigned int ii=0; ii<nwyz; ii++) {
for (unsigned int jj=0; jj<nwyz; jj++) {
T& val = vecData[jj];
if(val<0) val=-val;
}
} else {
for (unsigned int ii=0; ii<nwyz; ii++) {
for (unsigned int jj=0; jj<nwyz; jj++) {
T& val = vecData[jj];
if(val<0 && val!=IOUtils::nodata) val=-val;
}
......
......@@ -70,31 +70,27 @@ void BufferedIOHandler::read2DGrid(Grid2DObject& in_grid2Dobj, const std::string
return;
}
Grid2DObject tmpgrid2D; //HACK: why don't we read directly into in_grid2Dobj?
iohandler.read2DGrid(tmpgrid2D, in_filename);
bufferGrid(tmpgrid2D, in_filename);
in_grid2Dobj = tmpgrid2D;
iohandler.read2DGrid(in_grid2Dobj, in_filename);
bufferGrid(in_grid2Dobj, in_filename); //the STL containers make a copy
} else {
iohandler.read2DGrid(in_grid2Dobj, in_filename);
}
}
void BufferedIOHandler::read2DGrid(Grid2DObject& grid_out, const MeteoGrids::Parameters& parameter, const Date& date)
void BufferedIOHandler::read2DGrid(Grid2DObject& in_grid2Dobj, const MeteoGrids::Parameters& parameter, const Date& date)
{
if(max_grids>0) {
const string buffer_name = date.toString(Date::ISO)+"::"+MeteoGrids::getParameterName(parameter);
const std::map<std::string, Grid2DObject>::const_iterator it = mapBufferedGrids.find(buffer_name);
if (it != mapBufferedGrids.end()) { //already in map
grid_out = (*it).second;
in_grid2Dobj = (*it).second;
return;
}
Grid2DObject tmpgrid2D;
iohandler.read2DGrid(tmpgrid2D, parameter, date);
bufferGrid(tmpgrid2D, buffer_name);
grid_out = tmpgrid2D;
iohandler.read2DGrid(in_grid2Dobj, parameter, date);
bufferGrid(in_grid2Dobj, buffer_name); //the STL containers make a copy
} else {
iohandler.read2DGrid(grid_out, parameter, date);
iohandler.read2DGrid(in_grid2Dobj, parameter, date);
}
}
......@@ -106,21 +102,20 @@ void BufferedIOHandler::readDEM(DEMObject& in_grid2Dobj)
//already in map. If the update properties have changed,
//we copy the ones given in input and force the update of the object
const DEMObject::update_type in_ppt = (DEMObject::update_type)in_grid2Dobj.getUpdatePpt();
const DEMObject::slope_type in_slope_alg = (DEMObject::slope_type)in_grid2Dobj.getDefaultAlgorithm();
in_grid2Dobj = (*it).second;
const DEMObject::update_type buff_ppt = (DEMObject::update_type)in_grid2Dobj.getUpdatePpt();
if(in_ppt!=buff_ppt) {
const DEMObject::slope_type buff_slope_alg = (DEMObject::slope_type)in_grid2Dobj.getDefaultAlgorithm();
if(in_ppt!=buff_ppt || in_slope_alg!=buff_slope_alg) {
in_grid2Dobj.setDefaultAlgorithm(in_slope_alg);
in_grid2Dobj.setUpdatePpt(in_ppt);
in_grid2Dobj.update();
}
return;
}
DEMObject tmpgrid2D;
//copy the updating policy of the destination
tmpgrid2D.setUpdatePpt((DEMObject::update_type)in_grid2Dobj.getUpdatePpt());
iohandler.readDEM(tmpgrid2D);
mapBufferedGrids["/:DEM"] = tmpgrid2D;
in_grid2Dobj = tmpgrid2D;
iohandler.readDEM(in_grid2Dobj);
mapBufferedGrids["/:DEM"] = in_grid2Dobj; //the STL containers make a copy
} else {
iohandler.readDEM(in_grid2Dobj);
}
......@@ -135,10 +130,8 @@ void BufferedIOHandler::readLanduse(Grid2DObject& in_grid2Dobj)
return;
}
Grid2DObject tmpgrid2D;
iohandler.readLanduse(tmpgrid2D);
mapBufferedGrids["/:LANDUSE"] = tmpgrid2D;
in_grid2Dobj = tmpgrid2D;
iohandler.readLanduse(in_grid2Dobj);
mapBufferedGrids["/:LANDUSE"] = in_grid2Dobj; //the STL containers make a copy
} else {
iohandler.readLanduse(in_grid2Dobj);
}
......@@ -154,10 +147,8 @@ void BufferedIOHandler::readAssimilationData(const Date& in_date, Grid2DObject&
return;
}
Grid2DObject tmpgrid2D;
iohandler.readAssimilationData(in_date, tmpgrid2D);
mapBufferedGrids["/:ASSIMILATIONDATA" + in_date.toString(Date::FULL)] = tmpgrid2D;
in_grid2Dobj = tmpgrid2D;
iohandler.readAssimilationData(in_date, in_grid2Dobj);
mapBufferedGrids["/:ASSIMILATIONDATA" + in_date.toString(Date::FULL)] = in_grid2Dobj; //the STL containers make a copy
} else {
iohandler.readAssimilationData(in_date, in_grid2Dobj);
}
......
......@@ -285,6 +285,13 @@ void DEMObject::setDefaultAlgorithm(const slope_type& i_algorithm) {
}
}
/**
* @brief Get the default slope calculation algorithm
* @return default algorithm to use for slope computation
*/
int DEMObject::getDefaultAlgorithm() const {
return dflt_algorithm;
}
/**
* @brief Recomputes the min/max of altitude, slope and curvature
* It return +/- std::numeric_limits\<double\>\:\:max() for a given parameter if its grid was empty/undefined
......@@ -401,7 +408,7 @@ double DEMObject::horizontalDistance(Coords point1, const Coords& point2)
point1.copyProj(point2);
}
return horizontalDistance(point1.getEasting(), point1.getNorthing(),
point2.getEasting(), point2.getNorthing() );
point2.getEasting(), point2.getNorthing() );
}
......@@ -630,7 +637,7 @@ void DEMObject::getHorizon(const Coords& point, const double& increment, std::ve
void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
//This computes the slope and the aspect at a given cell as well as the x and y components of the normal vector
double A[4][4]; //table to store neigbouring heights: 3x3 matrix but we want to start at [1][1]
//we use matrix notation: A[y][x]
if(algorithm==DFLT) {
algorithm = dflt_algorithm;
}
......@@ -805,11 +812,11 @@ void DEMObject::CalculateHorn(double A[4][4], double& o_slope, double& o_Nx, dou
}
void DEMObject::CalculateCorripio(double A[4][4], double& o_slope, double& o_Nx, double& o_Ny, double& o_Nz) {
//This calculates the surface normal vector using the two triangle method given in Corripio (2003) but using a 3x3 grid instead of 2x2
if ( A[1][1]!=IOUtils::nodata && A[1][3]!=IOUtils::nodata && A[3][3]!=IOUtils::nodata && A[3][1]!=IOUtils::nodata) {
// See Corripio (2003), knowing that here we normalize the result (divided by Nz=cellsize*cellsize)
o_Nx = (0.5 * (A[1][3] - A[1][1] + A[3][1] - A[3][3]) ) / cellsize;
o_Ny = (0.5 * (A[1][3] - A[1][1] + A[1][3] - A[3][3]) ) / cellsize;
//This calculates the surface normal vector using the two triangle method given in Corripio (2003) but cell centered instead of node centered (ie using a 3x3 grid instead of 2x2)
if ( A[1][1]!=IOUtils::nodata && A[1][3]!=IOUtils::nodata && A[3][1]!=IOUtils::nodata && A[3][3]!=IOUtils::nodata) {
// See Corripio (2003), knowing that here we normalize the result (divided by Nz=cellsize*cellsize) and that we are cell centered instead of node centered
o_Nx = (A[3][1] + A[1][1] - A[3][3] - A[1][3]) / (2.*2.*cellsize);
o_Ny = (A[3][1] - A[1][1] + A[3][3] - A[1][3]) / (2.*2.*cellsize);
o_Nz = 1.;
//There is no difference between slope = acos(n_z/|n|) and slope = atan(sqrt(sx*sx+sy*sy))
//slope = acos( (Nz / sqrt( Nx*Nx + Ny*Ny + Nz*Nz )) );
......
......@@ -91,6 +91,7 @@ class DEMObject : public Grid2DObject {
const bool& i_update=true, const slope_type& i_algorithm=DFLT);
void setDefaultAlgorithm(const slope_type& i_algorithm);
int getDefaultAlgorithm() const;
void setUpdatePpt(const update_type& in_update_flag);
int getUpdatePpt() const;
......
......@@ -14,58 +14,59 @@ string prefix_ref("ref_");
// --- Results ---
// Header: min max mean s. min s. max [s. == slope]
double r_DEM[] ={ 193, 4204.81, 1302.38, 0., 51.6848};
double r_SUB_DEM[] ={ 403.4, 4027.3, 1291.28, 0., 51.6848};
double r_DEM[] ={ 193., 4204.81, 1302.38, 0., 43.486};
double r_SUB_DEM[] ={ 403.4, 4027.3, 1291.28, 0., 40.0628};
// controll basic values on dem
bool simpleDEMcontroll(DEMObject& dem, double results[]){
bool status = true;
if(!IOUtils::checkEpsilonEquality(dem.grid2D.getMin(), results[0], epsilon)){
cerr << "error on Min : " << dem.grid2D.getMin() << " != " << results[0] << endl;
exit(1);
status=false;
}
if(!IOUtils::checkEpsilonEquality(dem.grid2D.getMax(), results[1], epsilon)){
cerr << "error on Max : " << dem.grid2D.getMax() << " != " << results[1]<< endl;
exit(1);
status=false;
}
if(!IOUtils::checkEpsilonEquality(dem.grid2D.getMean(), results[2], 1.0e-2)){ // HACK adapt epsilon for test...
cerr << "error on Mean : " << dem.grid2D.getMean() << " != " << results[2]<< endl;
exit(1);
status=false;
}
if(!IOUtils::checkEpsilonEquality(dem.min_slope, results[3], epsilon)){
cerr << "error on Slope Min : " << dem.min_slope << " != " << results[3]<< endl;
exit(1);
status=false;
}
if(!IOUtils::checkEpsilonEquality(dem.max_slope, results[4], 1.0e-4)){ // HACK adapt epsilon for test ....
cerr << "error on Slope Max : " << dem.max_slope << " != " << results[4]<< endl;
exit(1);
status=false;
}
return true;
return status;
}
// Make output files
bool makeDEMfiles(){
cout << " ----- Read DEM, make subfiles and some basic controll \n";
cout << " ---- Init Values \n";
DEMObject dem;
Config cfg("io.ini");
IOManager io(cfg);
cout << " ---- If output files exist, empty them \n"; // HACK need to make this ???
for(int i = 0; i < n_files; i++){
ofstream ofs(files[i].c_str(), ofstream::out | ofstream::trunc);
ofs << " " ;
ofs.close();
}
cout << " ---- Read DEM \n";
dem.setUpdatePpt(DEMObject::SLOPE);
io.readDEM(dem);
......@@ -73,11 +74,11 @@ bool makeDEMfiles(){
cout << " ---- Generate file with slope values \n";
Grid2DObject slope(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.slope);
io.write2DGrid(slope, MeteoGrids::SLOPE, Date(0.));
cout << " ---- Generate file with azimut values \n";
Grid2DObject azi(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner, dem.azi);
io.write2DGrid(azi, MeteoGrids::AZI, Date(0.));
cout << " ---- Gridify \n"; // HACK how to controll, wath are points possible to controll
//retrieving grid coordinates of a real world point
Coords point;
......@@ -85,7 +86,7 @@ bool makeDEMfiles(){
point.setLatLon(46.232103, 7.362185, IOUtils::nodata);
dem.gridify(point);
cout << " ---- Get out a Sub dem and writh it in file \n";
cout << " ---- Get out a Sub dem and writh it in file \n";
//computing grid distances from real world distances
const double dist_x=70000., dist_y=120000.;
const unsigned int ncols = (unsigned int)ceil(dist_x/dem.cellsize);
......@@ -93,69 +94,71 @@ bool makeDEMfiles(){
DEMObject sub_dem(dem, point.getGridI(), point.getGridJ(), ncols, nrows);
sub_dem.grid2D.setKeepNodata(true); //HACK this removes error, but should it not take it in constructor ???
io.write2DGrid(sub_dem, MeteoGrids::DEM, Date(0.));
cout << " --- Controll some basic data from DEM\n";
bool status=true;
cout << " --- Controll some basic data from DEM\n";
if(!simpleDEMcontroll(dem, r_DEM)){
cerr << " Error on controlling basic values of DEM object !!! \n";
exit(1);
status=false;
}
cout << " --- Controll some basic data from SUB_DEM\n";
if(!simpleDEMcontroll(sub_dem, r_SUB_DEM)){
cerr << " Error on controlling basic values of DEM object !!! \n";
exit(1);
status=false;
}
return true;
return status;
}
bool compareFiles(){
bool status=true;
for(int i = 0; i < n_files; i++){
// ------ Compare reference file with generated results ---------
cout << " --- Start comparing reference file with output file generated before for : " << files[i] << endl;
ifstream ifref((prefix_ref+files[i]).c_str());
ifstream ifout(files[i].c_str());
string l_ref, l_out;
while (!ifref.eof())
{
if(ifout.eof()){
cerr << "Not enough lines generated as result!!!" << endl;
exit(1);
return false;
}
getline(ifref,l_ref);
getline(ifout,l_out);
if (l_ref!=l_out) {
cerr << " ERROR, Sun generatet error at following point an error " << endl;
cerr << "ref : \n " << l_ref << endl;
cerr << "out : \n " << l_out << endl;
exit(1);
status=false;
}
}
if(!ifout.eof()){
cerr << "To much lines generated as result!!!" << endl;
exit(1);
status=false;
}
ifout.close();
ifout.close();
ifref.close();
}
return true;
return status;
}
int main(void) {
if(!makeDEMfiles()){
exit(1);
}
if(!compareFiles()){
exit(1);
}
return 0;
}
This diff is collapsed.
This diff is collapsed.
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