WSL/SLF GitLab Repository

Commit 835792f6 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

Better handling of the i,j,k index in Coords objects, better strategy for...

Better handling of the i,j,k index in Coords objects, better strategy for handling invalid POI (the invalid points can be marked as invalid but kept in the vector)
parent 461b49bc
This diff is collapsed.
......@@ -125,6 +125,7 @@ public:
int getGridI() const;
int getGridJ() const;
int getGridK() const;
bool indexIsValid() const;
void getProj(std::string& proj_type, std::string& proj_args) const;
short int getEPSG() const;
......@@ -136,7 +137,7 @@ public:
void setLatLon(const double in_latitude, const double in_longitude, const double in_altitude, const bool in_update=true);
void setLatLon(const std::string& in_coordinates, const double in_altitude, const bool in_update=true);
void setXY(const double in_easting, const double in_northing, const double in_altitude, const bool in_update=true);
void setGridIndex(const int in_grid_i, const int in_grid_j, const int in_grid_k, const bool in_invalidate=true);
void setGridIndex(const int in_grid_i, const int in_grid_j, const int in_grid_k, const bool setValid=false);
void setAltitude(const double in_altitude, const bool in_update=true);
void setProj(const std::string& in_coordinatesystem, const std::string& in_parameters="");
void setLocalRef(const double in_ref_latitude, const double in_ref_longitude);
......@@ -176,6 +177,7 @@ public:
int grid_i; ///<grid index i (please notice that this index is NOT automatically regenerated NOR checked)
int grid_j; ///<grid index j (please notice that this index is NOT automatically regenerated NOR checked)
int grid_k; ///<grid index k (please notice that this index is NOT automatically regenerated NOR checked)
bool validIndex; ///< are grid index invalid?
std::string coordsystem;
std::string coordparam;
......
......@@ -31,6 +31,16 @@
namespace mio {
const struct CoordsAlgorithms::ELLIPSOID CoordsAlgorithms::ellipsoids[6] = {
{ 6378137., 6356752.3142 }, ///< E_WGS84
{ 6378137., 6356752.3141 }, ///< E_GRS80
{ 6377563.396, 6356256.909 }, ///< E_AIRY
{ 6378388., 6356911.946 }, ///< E_INTL1924
{ 6378249.145, 6356514.86955 }, ///< E_CLARKE1880
{ 6378160., 6356774.719 } ///< E_GRS67
};
/**
* @brief Print a nicely formatted lat/lon in degrees, minutes, seconds
* @return lat/lon
......
......@@ -23,7 +23,7 @@ using namespace std;
namespace mio {
Grid2DObject& Grid2DObject::operator=(const Grid2DObject& source) {
if(this != &source) {
if (this != &source) {
grid2D = source.grid2D;
cellsize = source.cellsize;
llcorner = source.llcorner;
......@@ -177,7 +177,7 @@ Grid2DObject::Grid2DObject(const Grid2DObject& i_grid2Dobj, const size_t& i_nx,
//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;
if( (llcorner.getEasting()!=IOUtils::nodata) && (llcorner.getNorthing()!=IOUtils::nodata) ) {
if ( (llcorner.getEasting()!=IOUtils::nodata) && (llcorner.getNorthing()!=IOUtils::nodata) ) {
llcorner.setXY( llcorner.getEasting()+static_cast<double>(i_nx)*i_grid2Dobj.cellsize,
llcorner.getNorthing()+static_cast<double>(i_ny)*i_grid2Dobj.cellsize, IOUtils::nodata);
}
......@@ -189,7 +189,7 @@ bool Grid2DObject::gridify(std::vector<Coords>& vec_points, const bool& keep_inv
std::vector<Coords>::iterator v_Itr = vec_points.begin();
while ( v_Itr != vec_points.end() ) {
if( gridify(*v_Itr)==false ) {
if ( gridify(*v_Itr)==false ) {
if (!keep_invalid) v_Itr = vec_points.erase(v_Itr);
status=false;
} else {
......@@ -203,12 +203,12 @@ bool Grid2DObject::gridify(std::vector<Coords>& vec_points, const bool& keep_inv
bool Grid2DObject::gridify(Coords& point) const {
std::string proj_type, proj_args;
point.getProj(proj_type, proj_args);
if(proj_type=="NULL") {
if (proj_type=="NULL") {
//if the projection was "NULL", we set it to the grid's
point.copyProj(llcorner);
}
if(point.getGridI()!=IOUtils::inodata && point.getGridJ()!=IOUtils::inodata) {
if (point.getGridI()!=IOUtils::inodata && point.getGridJ()!=IOUtils::inodata) {
//we need to compute (easting,northing) and (lat,lon)
return( grid_to_WGS84(point) );
} else {
......@@ -220,21 +220,22 @@ bool Grid2DObject::gridify(Coords& point) const {
bool Grid2DObject::grid_to_WGS84(Coords& point) const {
int i = point.getGridI(), j = point.getGridJ();
if(i==IOUtils::inodata || j==IOUtils::inodata) {
if (i==IOUtils::inodata || j==IOUtils::inodata) {
//the point is invalid (outside the grid or contains nodata)
point.setGridIndex(IOUtils::inodata, IOUtils::inodata, IOUtils::inodata, false); //mark the point as invalid
return false;
}
const int ncols = (signed)getNx();
const int nrows = (signed)getNy();
if(i>ncols || i<0 || j>nrows || j<0) {
if (i>ncols || i<0 || j>nrows || j<0) {
//the point is outside the grid, we reset the indices to the closest values
//still fitting in the grid and return an error
if(i<0) i=0;
if(j<0) j=0;
if(i>ncols) i=ncols;
if(j>nrows) j=nrows;
point.setGridIndex(i, j, IOUtils::inodata, false);
if (i<0) i=0;
if (j<0) j=0;
if (i>ncols) i=ncols;
if (j>nrows) j=nrows;
point.setGridIndex(i, j, IOUtils::inodata, false); //mark the point as invalid
return false;
}
......@@ -242,7 +243,7 @@ bool Grid2DObject::grid_to_WGS84(Coords& point) const {
const double easting = ((double)i) * cellsize + llcorner.getEasting();
const double northing = ((double)j) * cellsize + llcorner.getNorthing();
if(point.isSameProj(llcorner)==true) {
if (point.isSameProj(llcorner)==true) {
//same projection between the grid and the point -> precise, simple and efficient arithmetics
point.setXY(easting, northing, IOUtils::nodata);
} else {
......@@ -254,20 +255,21 @@ bool Grid2DObject::grid_to_WGS84(Coords& point) const {
point.copyProj(tmp_proj); //back to the original projection -> reproject the coordinates
}
point.setGridIndex(i, j, IOUtils::inodata, false);
point.setGridIndex(i, j, IOUtils::inodata, true); //mark the point as valid
return true;
}
bool Grid2DObject::WGS84_to_grid(Coords& point) const {
if(point.getLat()==IOUtils::nodata || point.getLon()==IOUtils::nodata) {
if (point.getLat()==IOUtils::nodata || point.getLon()==IOUtils::nodata) {
//if the point is invalid, there is nothing we can do
point.setGridIndex(IOUtils::inodata, IOUtils::inodata, IOUtils::inodata, false); //mark the point as invalid
return false;
}
bool error_code=true;
bool status=true;
int i,j;
if(point.isSameProj(llcorner)==true) {
if (point.isSameProj(llcorner)==true) {
//same projection between the grid and the point -> precise, simple and efficient arithmetics
i = (int)floor( (point.getEasting()-llcorner.getEasting()) / cellsize );
j = (int)floor( (point.getNorthing()-llcorner.getNorthing()) / cellsize );
......@@ -281,25 +283,25 @@ bool Grid2DObject::WGS84_to_grid(Coords& point) const {
//checking that the calculated indices fit in the grid2D
//and giving them the closest value within the grid if not.
if(i<0) {
if (i<0) {
i=0;
error_code=false;
status=false;
}
if(i>(signed)getNx()) {
if (i>(signed)getNx()) {
i=(signed)getNx();
error_code=false;
status=false;
}
if(j<0) {
if (j<0) {
j=0;
error_code=false;
status=false;
}
if(j>(signed)getNy()) {
if (j>(signed)getNy()) {
j=(signed)getNy();
error_code=false;
status=false;
}
point.setGridIndex(i, j, IOUtils::inodata, false);
return error_code;
point.setGridIndex(i, j, IOUtils::inodata, status); //mark as valid or invalid according to status
return status;
}
void Grid2DObject::set(const size_t& i_ncols, const size_t& i_nrows,
......@@ -358,7 +360,7 @@ void Grid2DObject::setValues(const double& i_cellsize, const Coords& i_llcorner)
bool Grid2DObject::isSameGeolocalization(const Grid2DObject& target) const
{
if( grid2D.getNx()==target.grid2D.getNx() &&
if ( grid2D.getNx()==target.grid2D.getNx() &&
grid2D.getNy()==target.grid2D.getNy() &&
llcorner==target.llcorner &&
cellsize==target.cellsize) {
......@@ -383,7 +385,7 @@ bool Grid2DObject::clusterization(const std::vector<double>& thresholds, const s
if (val!=IOUtils::nodata){
size_t i = 0;
for ( ;i<nscl; i++)
if(thresholds[i] >= val)
if (thresholds[i] >= val)
break;
grid2D(jj) = ids[i];
}
......
......@@ -187,16 +187,14 @@ Grid3DObject::Grid3DObject(const double& i_cellsize, const Coords& i_llcorner, c
cellsize(i_cellsize), z_is_absolute(true) {}
bool Grid3DObject::gridify(std::vector<Coords>& vec_points, std::vector<Coords>& vec_invalid) const
bool Grid3DObject::gridify(std::vector<Coords>& vec_points, const bool& keep_invalid) const
{
bool status=true;
vec_invalid.clear();
std::vector<Coords>::iterator v_Itr = vec_points.begin();
while ( v_Itr != vec_points.end() ) {
if( gridify(*v_Itr)==false ) {
vec_invalid.push_back( *v_Itr );
v_Itr = vec_points.erase(v_Itr);
if ( gridify(*v_Itr)==false ) {
if (!keep_invalid) v_Itr = vec_points.erase(v_Itr);
status=false;
} else {
++v_Itr;
......@@ -206,12 +204,6 @@ bool Grid3DObject::gridify(std::vector<Coords>& vec_points, std::vector<Coords>&
return status;
}
bool Grid3DObject::gridify(std::vector<Coords>& vec_points) const
{
std::vector<Coords> vec_invalid;
return gridify(vec_points, vec_invalid);
}
bool Grid3DObject::gridify(Coords& point) const
{
std::string proj_type, proj_args;
......@@ -236,6 +228,7 @@ bool Grid3DObject::grid_to_WGS84(Coords& point) const
if(i==IOUtils::inodata || j==IOUtils::inodata || k==IOUtils::inodata) {
//the point is invalid (outside the grid or contains nodata)
point.setGridIndex(IOUtils::inodata, IOUtils::inodata, IOUtils::inodata, false); //mark the point as invalid
return false;
}
......@@ -251,7 +244,7 @@ bool Grid3DObject::grid_to_WGS84(Coords& point) const
if(i>ncols) i=ncols;
if(j>nrows) j=nrows;
if(k>ndepths) k=ndepths;
point.setGridIndex(i, j, k, false);
point.setGridIndex(i, j, k, false); //mark the point as invalid
return false;
}
......@@ -271,6 +264,8 @@ bool Grid3DObject::grid_to_WGS84(Coords& point) const
point.setXY(easting, northing, altitude);
point.copyProj(tmp_proj); //back to the original projection -> reproject the coordinates
}
point.setGridIndex(i, j, k, true); //mark the point as valid
return true;
}
......@@ -278,10 +273,10 @@ bool Grid3DObject::WGS84_to_grid(Coords point) const
{
if(point.getLat()==IOUtils::nodata || point.getLon()==IOUtils::nodata || point.getAltitude()==IOUtils::nodata) {
//if the point is invalid, there is nothing we can do
point.setGridIndex(IOUtils::inodata, IOUtils::inodata, IOUtils::inodata, false); //mark the point as invalid
return false;
}
bool error_code=true;
int i,j,k;
if(point.isSameProj(llcorner)==true) {
......@@ -298,35 +293,36 @@ bool Grid3DObject::WGS84_to_grid(Coords point) const
k = (int)floor( (point.getAltitude()-llcorner.getAltitude()) / cellsize );
}
bool status=true;
//checking that the calculated indices fit in the grid2D
//and giving them the closest value within the grid if not.
if(i<0) {
i=0;
error_code=false;
status=false;
}
if(i>(signed)grid3D.getNx()) {
i=(signed)grid3D.getNx();
error_code=false;
status=false;
}
if(j<0) {
j=0;
error_code=false;
status=false;
}
if(j>(signed)grid3D.getNy()) {
j=(signed)grid3D.getNy();
error_code=false;
status=false;
}
if(k<0) {
k=0;
error_code=false;
status=false;
}
if(k>(signed)grid3D.getNz()) {
k=(signed)grid3D.getNz();
error_code=false;
status=false;
}
point.setGridIndex(i, j, k, false);
return error_code;
point.setGridIndex(i, j, k, status); //mark as valid or invalid according to status
return status;
}
void Grid3DObject::set(const size_t& i_ncols, const size_t& i_nrows, const size_t& i_ndepths,
......
......@@ -134,12 +134,12 @@ class Grid3DObject{
* This means that the Coords::point object that is given either contains geographic coordinates or
* grid indices. This method will calculate the missing ones (so that (i,j) match with (lat,lon)
* and (east,north)). Any point that is either invalid or outside the grid is removed from the vector.
* If the given point had a "NULL" projection, it will be set to the grid's.
* @param vec_points vector containing the coordinates to convert
* @param vec_invalid vector containing the rejected coordinates
* @param keep_invalid keep invalid coordinates? (default: false)
* @return false if invalid or external points had to be removed
*/
bool gridify(std::vector<Coords>& vec_points, std::vector<Coords>& vec_invalid) const;
bool gridify(std::vector<Coords>& vec_points) const;
bool gridify(std::vector<Coords>& vec_points, const bool& keep_invalid=false) const;
/**
* @brief check if the current Grid3DObject has the same geolocalization attributes
......
......@@ -562,7 +562,7 @@ void calculate_dimensions(const mio::Grid2DObject& grid, double*& lat_array, dou
// corner to calculate the lat/lon intervals between cells
Coords urcorner(grid.llcorner);
urcorner.setGridIndex(static_cast<int>(grid.getNx() - 1), static_cast<int>(grid.getNy() - 1), IOUtils::nodata, true);
grid.gridify(urcorner);
grid.gridify(urcorner); //no need to check the return value: we know it fits within the grid
const double lat_interval = (urcorner.getLat() - lat_array[0]) / static_cast<double>(grid.getNy()-1);
const double lon_interval = (urcorner.getLon() - lon_array[0]) / static_cast<double>(grid.getNx()-1);
......
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