WSL/SLF GitLab Repository

Commit 687c52df authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The new Coords object is now used by all grids (Grid2DObject, Grid3DObject,...

The new Coords object is now used by all grids (Grid2DObject, Grid3DObject, DEMObject) as well as by all plugins (but checking proper operation was not possible for GSNIO because of missing libraries on my system...). The documentation has been updated as well as the code examples. A few extra methods have been added when usage has shown the need/convenience of having them (specially for standardazing some processing called by the plugins). A few small bugs have been fixed (possibility of not updating or improperly updating the coordinates in some rare cases).

The proper marshalling for Grid3DObject has been implemented (it was more or less a dummy method until now). One bug has been fixed for Proj4 support. The code of most of the plugins has been cleaned up in order to get rid of the few remaining "using namespace" in header files. 
parent b40dedb4
......@@ -28,5 +28,9 @@ int main(int argc, char** argv) {
//we print the (lat,long) and (x,y) coordinates once more to check that they still match what we had a the begining
printf("Lat/long to CH1903\n\t(%g , %g) -> (%g , %g)\n", point1.getLat(), point1.getLon(), point1.getEasting(), point1.getNorthing());
//A nice finishing touch: we print nicely formatted lat/lon
cout << "Pretty printing\n";
cout << "\t(" << point1.getLat() << " , " << point1.getLon() << ") = " << point1.printLatLon() << endl;
return 0;
}
......@@ -3,15 +3,16 @@
//This is a basic example of using as dem: the dem is read, the grid coordinates of a point given by its (lat,long) are retrieved
//and a sub-dem is extracted starting at these coordinates and extending dist_x and dist_y and written out.
int main(void) {
const double lat1=46.1592, lon1=8.12993;
const double dist_x=700, dist_y=1200;
DEMObject dem;
IOHandler *io = NULL;
ConfigReader *cfg = NULL;
unsigned int i,j;
try {
ConfigReader cfg("io.ini");
io = new IOHandler(cfg);
//ConfigReader cfg("io.ini");
cfg = new ConfigReader("io.ini");
io = new IOHandler(*cfg);
} catch (IOException& e){
std::cout << "Problem with IOHandler creation, cause: " << e.what() << std::endl;
}
......@@ -26,7 +27,9 @@ int main(void) {
std::cout << "\tmin slope=" << dem.min_slope << " max slope=" << dem.max_slope << std::endl;
//retrieving grid coordinates of a real world point
dem.WGS84_to_grid(lat1, lon1, i,j);
Coords point(*cfg);
point.setLatLon(46.1592, 8.12993);
dem.WGS84_to_grid(point, i,j);
//computing grid distances from real world distances
const unsigned int ncols = (unsigned int)ceil(dist_x/dem.cellsize);
......@@ -38,11 +41,5 @@ int main(void) {
//writing the sub-dem out
io->write2DGrid(sub_dem,"sub_dem.dem");
//Write a message giving the nicely formatted coordinates of the two opposite corners of the sub-dem
double lat2, lon2;
dem.grid_to_WGS84(i+ncols, j+nrows, lat2, lon2);
std::cout << "Sub-dem (" << MapProj::decimal_to_dms(lat1) << "," << MapProj::decimal_to_dms(lon1) << ") - ";
std::cout << "(" << MapProj::decimal_to_dms(lat2) << "," << MapProj::decimal_to_dms(lon2) << ") written" << std::endl;
return 0;
}
......@@ -149,11 +149,7 @@ void A3DIO::readMeteoData(const Date_IO& dateStart, const Date_IO& dateEnd,
void A3DIO::convertUnits(MeteoData& meteo)
{
for(unsigned int ii=0; ii<MeteoData::nrOfParameters; ii++){
//loop through all meteo params and check whether they're nodata values
if (meteo.param(ii)<=plugin_nodata)
meteo.param(ii) = IOUtils::nodata;
}
meteo.standardizeNodata(plugin_nodata);
//converts C to Kelvin, converts RH to [0,1]
if(meteo.ta!=IOUtils::nodata) {
......@@ -206,20 +202,11 @@ void A3DIO::read1DMeteo(const Date_IO& dateStart, const Date_IO& dateEnd,
IOUtils::getValueForKey(header, "Y_Coord", ycoord);
IOUtils::getValueForKey(header, "Altitude", altitude);
std::string coordsys="", coordparam="";
try {
cfg.getValue("COORDIN", coordsys);
cfg.getValue("COORDPARAM", coordparam, ConfigReader::nothrow);
} catch(std::exception& e) {
//problems while reading values for COORDIN or COORDPARAM
std::cerr << "[E] reading configuration file: " << "\t" << e.what() << std::endl;
throw;
}
//compute/check WGS coordinates (considered as the true reference)
Coords location(coordsys, coordparam);
location.setXY(xcoord, ycoord);
location.setLatLon(latitude, longitude);
//HACK!! FIXME!! we need to transfer the plugin nodata for the checks!!!!
//compute/check WGS coordinates (considered as the true reference) according to the projection as defined in cfg
Coords location(cfg);
location.setXY(xcoord, ycoord, false);
location.setLatLon(latitude, longitude, false);
try {
location.check();
} catch(...) {
......@@ -573,15 +560,6 @@ void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,
unsigned int columns = 0;
std::vector<std::string> vec_altitude, vec_xcoord, vec_ycoord, vec_names;
//Build Coords object to convert easting/northing values to lat/long in WGS84
std::string coordsys="", coordparam="";
try {
cfg.getValue("COORDIN", coordsys);
cfg.getValue("COORDPARAM", coordparam, ConfigReader::nothrow);
} catch(std::exception& e){
//problems while reading values for COORDIN or COORDPARAM
}
fin.clear();
fin.open (filename.c_str(), std::ifstream::in);
if (fin.fail()) {
......@@ -611,6 +589,9 @@ void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,
throw InvalidFormatException("Column count doesn't match from line to line in " + filename, AT);
}
//Build Coords object to convert easting/northing values to lat/long in WGS84
Coords coordinate(cfg);
for (unsigned int ii=4; ii<columns; ii++) {
unsigned int stationnr = hashStations[vec_names.at(ii)];
double altitude, easting, northing, stationName;
......@@ -620,16 +601,9 @@ void A3DIO::read2DMeteoHeader(const std::string& filename, std::map<std::string,
|| (!IOUtils::convertString(stationName, vec_names.at(ii), std::dec))) {
throw ConversionFailedException("Conversion of station description failed in " + filename, AT);
}
coordinate.setXY(easting, northing, altitude);
vecS[stationnr-1].stationName = stationName;
vecS[stationnr-1].position.setProj(coordsys, coordparam);
vecS[stationnr-1].position.setXY(easting, northing, altitude, false);
//calculate/check coordinates if necessary
try {
vecS[stationnr-1].position.check();
} catch (...) {
std::cerr << "[E] Too many nodata values for coordinates conversion in file " << filename << " at " << AT << std::endl;
}
vecS[stationnr-1].position = coordinate;
}
cleanup();
......
......@@ -35,8 +35,13 @@ const struct Coords::ELLIPSOID Coords::ellipsoids[] = {
{ 6378160., 6356774.719 } //E_GRS67
};
bool Coords::operator==(const Coords& in) const
{//check that two Coords objects represent the same location
/**
* @brief Equality operator that checks that lat/lon match
* @param[in] in Coord object to compare to
* @return true or false
*/
bool Coords::operator==(const Coords& in) const {
//check that two Coords objects represent the same location
const double earth_radius = 6371e3; //in meters
const double grid_epsilon = 5.; //in meters
const double long_epsilon = grid_epsilon / earth_radius; //in degrees. small angle, so sin(x)=x
......@@ -47,8 +52,12 @@ bool Coords::operator==(const Coords& in) const
return comparison;
}
bool Coords::operator!=(const Coords& in) const
{
/**
* @brief Inequality operator that checks that lat/lon don't match
* @param[in] in Coord object to compare to
* @return true or false
*/
bool Coords::operator!=(const Coords& in) const {
return !(*this==in);
}
......@@ -75,32 +84,46 @@ Coords& Coords::operator=(const Coords& source) {
* This constructor builds a dummy object that performs no conversions but can be used for comparison
* purpose. This is more or less the equivalent of NULL for a pointer...
*/
Coords::Coords()
{
Coords::Coords() {
initializeMaps();
setDefaultValues();
setProj("NULL", "NULL");
}
/**
* \anchor Coordinate_types
* @brief Regular constructor: usually, this is the constructor to use
* @param[in] _coordinatesystem string identifying the coordinate system to use
* @param[in] _parameters string giving some additional parameters for the projection (optional)
*
* The coordinate system can be any of the following:
* - CH1903 for coordinates in the Swiss Grid
* - UTM for UTM coordinates (the zone must be specified in the parameters, for example 31T)
* - PROJ4 for coordinate conversion relying on the Proj4 library
* - LOCAL for local coordinate system (from reference point)
* See setProj() for a full description of these strings
*/
Coords::Coords(const std::string& _coordinatesystem, const std::string& _parameters)
{
Coords::Coords(const std::string& _coordinatesystem, const std::string& _parameters) {
setDefaultValues();
initializeMaps();
setProj(_coordinatesystem, _parameters);
}
/**
* @brief Use the parameters included in a ConfigReader for setting the projection up
* @param[in] cfg ConfigReader containing the necessary parameters for initializing
* the projection
*/
Coords::Coords(const ConfigReader& cfg) {
std::string coordsys="", coordparam="";
try {
cfg.getValue("COORDIN", coordsys);
cfg.getValue("COORDPARAM", coordparam, ConfigReader::nothrow);
} catch(std::exception& e) {
//problems while reading values for COORDIN or COORDPARAM
std::cerr << "[E] reading configuration file: " << "\t" << e.what() << std::endl;
throw;
}
setDefaultValues();
initializeMaps();
setProj(coordsys, coordparam);
}
/**
* @brief Local projection onstructor: this constructor is only suitable for building a local projection.
* Such a projection defines easting and northing as the distance (in meters) to a reference point
......@@ -116,26 +139,67 @@ Coords::Coords(const double& _lat_ref, const double& _long_ref)
setProj("LOCAL", "");
}
/**
* @brief Returns the East coordinate in the configured projection system
* @return easting
*/
double Coords::getEasting() const {
return easting;
}
/**
* @brief Returns the North coordinate in the configured projection system
* @return northing
*/
double Coords::getNorthing() const {
return northing;
}
/**
* @brief Returns the Latitude in the configured projection system
* @return latitude
*/
double Coords::getLat() const {
return latitude;
}
/**
* @brief Returns the Latitude in the configured projection system
* @return longitude
*/
double Coords::getLon() const {
return longitude;
}
/**
* @brief Returns the Altitude. This is currently independent of the configured projection system
* @return altitude
*/
double Coords::getAltitude() const {
return altitude;
}
/**
* @brief Print a nicely formatted lat/lon in degrees, minutes, seconds
* @return lat/lon
*/
std::string Coords::printLatLon() const {
std::stringstream dms;
dms << "(" << decimal_to_dms(latitude) << " , " << decimal_to_dms(longitude) << ")";
return dms.str();
}
/**
* @brief Set latitude and longitude
* The automatic update of the easting/northing can be turned off so that
* both lat/lon and east/north coordinates can be provided in order to thereafter check the
* coordinates by calling check().
* @param[in] _latitude latitude to set
* @param[in] _longitude longitude to set
* @param[in] _altitude altitude to set (optional)
* @param[in] _update should the easting/northing be updated? (default=true)
*/
void Coords::setLatLon(const double _latitude, const double _longitude, const double _altitude, const bool _update) {
latitude = _latitude;
longitude = _longitude;
......@@ -147,6 +211,31 @@ void Coords::setLatLon(const double _latitude, const double _longitude, const do
}
}
/**
* @brief Set latitude and longitude
* The automatic update of the easting/northing can be turned off so that
* both lat/lon and east/north coordinates can be provided in order to thereafter check the
* coordinates by calling check().
* @param[in] _coordinates string containing the lat/lon to read
* @param[in] _altitude altitude to set (optional)
* @param[in] _update should the easting/northing be updated? (default=true)
*/
void Coords::setLatLon(const std::string& _coordinates, const double _altitude, const bool _update) {
double lat, lon;
parseLatLon(_coordinates, lat, lon);
setLatLon(lat, lon, _altitude, _update);
}
/**
* @brief Set easting and northing
* The automatic update of the latitude/longitude can be turned off so that
* both lat/lon and east/north coordinates can be provided in order to thereafter check the
* coordinates by calling check().
* @param[in] _easting easting to set
* @param[in] _northing northing to set
* @param[in] _altitude altitude to set (optional)
* @param[in] _update should the easting/northing be updated? (default=true)
*/
void Coords::setXY(const double _easting, const double _northing, const double _altitude, const bool _update) {
easting = _easting;
northing = _northing;
......@@ -158,11 +247,23 @@ void Coords::setXY(const double _easting, const double _northing, const double _
}
}
/**
* \anchor Coordinate_types
* @brief Set projection to use
* This projection will be used for converting between lat/lon and East/North
* @param[in] _coordinatesystem string identifying the coordinate system to use
* @param[in] _parameters string giving some additional parameters for the projection (optional)
*
* The coordinate system can be any of the following:
* - CH1903 for coordinates in the Swiss Grid
* - UTM for UTM coordinates (the zone must be specified in the parameters, for example 31T)
* - PROJ4 for coordinate conversion relying on the Proj4 library
* - LOCAL for local coordinate system (from reference point)
*/
void Coords::setProj(const std::string& _coordinatesystem, const std::string& _parameters) {
//the latitude/longitude had not been calculated, so we do it first in order to have our reference
//before further conversions (usage scenario: giving a x,y and then converting to anyother x,y in another system
if(coordsystem!="NULL") {
//the convert method will check for nodata
if(coordsystem!="NULL" && (latitude==IOUtils::nodata) || (longitude==IOUtils::nodata) ) {
convert_to_WGS84(easting, northing, latitude, longitude);
}
......@@ -170,10 +271,18 @@ void Coords::setProj(const std::string& _coordinatesystem, const std::string& _p
coordparam = _parameters;
setFunctionPointers();
//since lat/long is our reference, we refresh x,y
convert_from_WGS84(latitude, longitude, easting, northing);
//since lat/long is our reference, we refresh x,y (only if lat/lon exist)
if(latitude!=IOUtils::nodata && longitude!=IOUtils::nodata) {
convert_from_WGS84(latitude, longitude, easting, northing);
}
}
/**
* @brief Set the local projection reference coordinates
* This projection will be used for converting between lat/lon and East/North
* @param[in] _ref_latitude latitude of the local origin
* @param[in] _ref_longitude longitude of the local origin
*/
void Coords::setLocalRef(const double _ref_latitude, const double _ref_longitude) {
if(_ref_latitude==IOUtils::nodata || _ref_longitude==IOUtils::nodata) {
throw InvalidArgumentException("For LOCAL projection, please provide both reference latitude and longitude!", AT);
......@@ -185,6 +294,11 @@ void Coords::setLocalRef(const double _ref_latitude, const double _ref_longitude
}
}
/**
* @brief Set the local projection reference coordinates
* This projection will be used for converting between lat/lon and East/North
* @param[in] _coordparam string containing the (lat,lon) of the local origin
*/
void Coords::setLocalRef(const std::string _coordparam) {
coordparam = _coordparam;
parseLatLon(coordparam, ref_latitude, ref_longitude);
......@@ -193,6 +307,11 @@ void Coords::setLocalRef(const std::string _coordparam) {
}
}
/**
* @brief Set the algorithm to use to compute distances
* Various algorithm exist that offer various precision/complexity tradeoffs.
* @param[in] _algo enum giving the algorithm to be used (see documentation for geo_distances)
*/
void Coords::setDistances(const geo_distances _algo) {
distance_algo = _algo;
if(coordsystem=="LOCAL") {
......@@ -244,12 +363,51 @@ void Coords::check()
* @param destination destination coordinate
* @return distance in meters
*/
double Coords::distance(Coords& destination) {
double Coords::distance(const Coords& destination) const {
double dist, bearing;
distance(destination, dist, bearing);
return dist;
}
/**
* @brief Check if two Coords object are using the same projection
* @param target coordinate to compare to
* @return true or false
*/
bool Coords::isSameProj(const Coords& target) const {
if(coordsystem=="LOCAL") {
return ( target.coordsystem=="LOCAL" && ref_latitude==target.ref_latitude && ref_longitude==target.ref_longitude );
} else {
return ( coordsystem==target.coordsystem && coordparam==target.coordparam );
}
}
/**
* @brief Copy the projection parameters of another Coords object
* @param source source object to copy the projection from
* @param _update should the necessary coordinates be updated? (default=true)
*/
void Coords::copyProj(const Coords& source, const bool _update) {
if(coordsystem=="LOCAL") {
coordsystem="LOCAL";
coordparam=source.coordparam;
ref_latitude=source.ref_latitude;
ref_longitude=source.ref_longitude;
} else {
coordsystem=source.coordsystem;
coordparam=source.coordparam;
}
setFunctionPointers();
if(_update==true) {
if((latitude!=IOUtils::nodata) && (longitude!=IOUtils::nodata)) {
convert_from_WGS84(latitude, longitude, easting, northing);
} else {
convert_to_WGS84(easting, northing, latitude, longitude);
}
}
}
/////////////////////////////////////////////////////private methods
/**
* @brief Method converting towards WGS84
......@@ -342,7 +500,7 @@ double Coords::dms_to_decimal(const std::string& dms) {
* @param[out] lon parsed longitude
*/
void Coords::parseLatLon(const std::string& coordinates, double&
lat, double& lon) {
lat, double& lon) const {
char lat_str[32]="";
char lon_str[32]="";
......@@ -631,7 +789,7 @@ void Coords::UTM_to_WGS84(double east_in, double north_in, double& lat_out, doub
void Coords::WGS84_to_PROJ4(double lat_in, double long_in, double& east_out, double& north_out) const
{
#ifdef PROJ4
const string src_param="+proj=latlong +datum=WGS84 +ellps=WGS84";
const std::string src_param="+proj=latlong +datum=WGS84 +ellps=WGS84";
projPJ pj_latlong, pj_dest;
double x=long_in*DEG_TO_RAD, y=lat_in*DEG_TO_RAD;
......@@ -707,7 +865,7 @@ void Coords::PROJ4_to_WGS84(double east_in, double north_in, double& lat_out, do
#endif
}
void Coords::distance(Coords& destination, double& distance, double& bearing) {
void Coords::distance(const Coords& destination, double& distance, double& bearing) const {
switch(distance_algo) {
case GEO_COSINE:
distance = cosineDistance(latitude, longitude, destination.getLat(), destination.getLon(), bearing);
......
......@@ -20,6 +20,7 @@
#include "IOExceptions.h"
#include "IOUtils.h"
#include "ConfigReader.h"
#include <string>
#include <map>
......@@ -44,6 +45,7 @@ class Coords {
//Constructors
Coords();
Coords(const std::string& coordinatesystem, const std::string& parameters="");
Coords(const ConfigReader& cfg);
Coords(const double& _lat_ref, const double& _long_ref);
//Operators
......@@ -57,21 +59,24 @@ class Coords {
double getLat() const;
double getLon() const;
double getAltitude() const;
std::string printLatLon() const;
//Setter methods
void setLatLon(const double _latitude, const double _longitude, const double _altitude=IOUtils::nodata, const bool _update=true);
void setLatLon(const std::string& _coordinates, const double _altitude=IOUtils::nodata, const bool _update=true);
void setXY(const double _easting, const double _northing, const double _altitude=IOUtils::nodata, const bool _update=true);
void setProj(const std::string& _coordinatesystem, const std::string& _parameters);
void setProj(const std::string& _coordinatesystem, const std::string& _parameters="");
void setLocalRef(const double _ref_latitude, const double _ref_longitude);
void setLocalRef(const std::string _coordparam);
void setDistances(const geo_distances _algo);
void check();
double distance(Coords& destination);
double distance(const Coords& destination) const;
bool isSameProj(const Coords& target) const;
void copyProj(const Coords& target, const bool _update=true);
//Static helper methods
static double dms_to_decimal(const std::string& dms);
static void parseLatLon(const std::string& coordinates, double& lat, double& lon);
static std::string decimal_to_dms(const double& decimal);
private:
......@@ -91,7 +96,7 @@ class Coords {
void NULL_to_WGS84(double east_in, double north_in, double& lat_out, double& long_out) const;
//Distances calculations
void distance(Coords& destination, double& distance, double& bearing);
void distance(const Coords& destination, double& distance, double& bearing) const;
double cosineDistance(const double& lat1, const double& lon1, const double& lat2, const double& lon2, double& alpha) const;
void cosineInverse(const double& lat_ref, const double& lon_ref, const double& distance, const double& bearing, double& lat, double& lon) const;
double VincentyDistance(const double& lat1, const double& lon1, const double& lat2, const double& lon2, double& alpha) const;
......@@ -102,6 +107,7 @@ class Coords {
void initializeMaps();
void setFunctionPointers();
int getUTMZone(const double latitude, const double longitude, std::string& zone_out) const;
void parseLatLon(const std::string& coordinates, double& lat, double& lon) const;
private:
double altitude;
......
......@@ -43,18 +43,13 @@ DEMObject::DEMObject(const slope_type& _algorithm) : Grid2DObject(), slope(), az
* @brief Constructor that sets variables.
* @param _ncols (unsigned int) number of colums in the grid2D
* @param _nrows (unsigned int) number of rows in the grid2D
* @param _xllcorner (double) x-coordinate of lower left corner
* @param _yllcorner (double) y-coordinate of lower left corner
* @param _latitude (double) decimal latitude, can be IOUtils::nodata
* @param _longitude (double) decimal longitude, can be IOUtils::nodata
* @param _cellsize (double) value for cellsize in grid2D
* @param _llcorner lower lower corner point
* @param _algorithm specify the default algorithm to use for slope computation (default=DFLT)
*/
DEMObject::DEMObject(const unsigned int& _ncols, const unsigned int& _nrows,
const double& _xllcorner, const double& _yllcorner,
const double& _latitude, const double& _longitude,
const double& _cellsize, const slope_type& _algorithm)
: Grid2DObject(_ncols, _nrows, _xllcorner, _yllcorner, _latitude, _longitude, _cellsize),
const double& _cellsize, const Coords& _llcorner, const slope_type& _algorithm)
: Grid2DObject(_ncols, _nrows, _cellsize, _llcorner),
slope(), azi(), curvature(), Nx(), Ny(), Nz()
{
min_altitude = min_slope = min_curvature = std::numeric_limits<double>::max();
......@@ -67,21 +62,16 @@ DEMObject::DEMObject(const unsigned int& _ncols, const unsigned int& _nrows,
* @brief Constructor that sets variables.
* @param _ncols (unsigned int) number of colums in the grid2D
* @param _nrows (unsigned int) number of rows in the grid2D
* @param _xllcorner (double) x-coordinate of lower left corner
* @param _yllcorner (double) y-coordinate of lower left corner
* @param _latitude (double) decimal latitude, can be IOUtils::nodata
* @param _longitude (double) decimal longitude, can be IOUtils::nodata
* @param _cellsize (double) value for cellsize in grid2D
* @param _llcorner lower lower corner point
* @param _altitude (Array2D\<double\>&) grid2D of elevations
* @param _update (bool) also update slope/normals/curvatures and their min/max? (default=true)
* @param _algorithm specify the default algorithm to use for slope computation (default=DFLT)
*/
DEMObject::DEMObject(const unsigned int& _ncols, const unsigned int& _nrows,
const double& _xllcorner, const double& _yllcorner,
const double& _latitude, const double& _longitude,
const double& _cellsize, const Array2D<double>& _altitude,
const double& _cellsize, const Coords& _llcorner, const Array2D<double>& _altitude,
const bool& _update, const slope_type& _algorithm)
: Grid2DObject(_ncols, _nrows, _xllcorner, _yllcorner, _latitude, _longitude, _cellsize, _altitude),
: Grid2DObject(_ncols, _nrows, _cellsize, _llcorner, _altitude),
slope(), azi(), curvature(), Nx(), Ny(), Nz()
{
slope_failures = curvature_failures = 0;
......@@ -100,7 +90,7 @@ DEMObject::DEMObject(const unsigned int& _ncols, const unsigned int& _nrows,
* @param _algorithm specify the default algorithm to use for slope computation (default=DFLT)
*/
DEMObject::DEMObject(const Grid2DObject& _dem, const bool& _update, const slope_type& _algorithm)
: Grid2DObject(_dem.ncols, _dem.nrows, _dem.xllcorner, _dem.yllcorner, _dem.latitude, _dem.longitude, _dem.cellsize, _dem.grid2D),
: Grid2DObject(_dem.ncols, _dem.nrows, _dem.cellsize, _dem.llcorner, _dem.grid2D),
slope(), azi(), curvature(), Nx(), Ny(), Nz()
{
slope_failures = curvature_failures = 0;
......@@ -284,22 +274,37 @@ double DEMObject::horizontalDistance(const double& xcoord1, const double& ycoord
return sqrt( IOUtils::pow2(xcoord2-xcoord1) + IOUtils::pow2(ycoord2-ycoord1) );
}
/**
* @brief Computes the horizontal distance between two points in a metric grid
* @param point1 first point (ie: origin)
* @param point2 second point (ie: destination)
* @return horizontal distance in meters
*
*/
double DEMObject::horizontalDistance(Coords point1, const Coords& point2)
{
if(point1.isSameProj(point2)==false) {
point1.copyProj(point2);
}
return horizontalDistance(point1.getEasting(), point1.getNorthing(),
point2.getEasting(), point2.getNorthing() );
}
/**
* @brief Returns the distance *following the terrain* between two coordinates
* @param xcoord1 east coordinate of the first point
* @param ycoord1 north coordinate of the first point
* @param xcoord2 east coordinate of the second point
* @param ycoord2 north coordinate of the second point
* @param point1 first point (ie: origin)
* @param point2 second point (ie: destination)
* @return distance following the terrain in meters
*
*/
double DEMObject::terrainDistance(const double& xcoord1, const double& ycoord1, const double& xcoord2, const double& ycoord2) {
double DEMObject::terrainDistance(Coords point1, const Coords& point2) {
std::vector<POINT> vec_points;