Coords.cc

/***********************************************************************************/
/*  Copyright 2009 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    MeteoIO is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <cmath>
#include <cstdio>
#include <iomanip>

#include <meteoio/dataClasses/Coords.h>
#include <meteoio/MathOptim.h>
#include <meteoio/meteoLaws/Meteoconst.h> //for math constants

#ifdef PROJ4
	#include <proj_api.h>
#endif

using namespace std;

namespace mio {
 /**
 * @page coords Available coordinate systems
 * Geographic coordinates will be transparently and automatically converted to lat/lon and any other coordinate system that
 * the client program uses. However, in order to do so, the input coordinate system must be specified. In order to output
 * geolocalized data, the desired coordinate system must also be specified for the outputs (in the output section).
 * This is done through the use of the COORDIN and COORDPARAM keys (see the documentation for each plugin).
 *
 * There are two ways of supporting a given coordinate system: through the use of an adhoc implementation
 * (that becomes part of MeteoIO) or through the use of an external library, Proj4 [ref: http://trac.osgeo.org/proj/].
 * The current internal implementations are the following (given by their keyword):
 * - CH1903 for coordinates in the Swiss Grid [ref: http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf]
 * - UTM for UTM coordinates (the zone must be specified in the parameters, for example 31T) [ref: http://www.oc.nps.edu/oc2902w/maps/utmups.pdf]
 * - UPS for Universal Polar Stereographic coordinates (the zone, either N or S, must be specified in the parameters). [ref: J. Hager, J. Behensky, B. Drew, <i>THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)</i>, 1989, Defense Mapping Agency, DMATM 8358.2]
 * - LOCAL for local coordinate system (using the horizontal and vertical distance from a reference point, see Coords::geo_distances for the available choice of distance algorithms)
 *
 * Such an example of use is the following:
 * @code
 * COORDSYS	= UTM
 * COORDPARAM	= 31T
 * @endcode
 *
 * On the other hand, when using the Proj4 library for handling the coordinate conversion, the EPSG codes of
 * the chosen projection must be specified (such codes can be found at http://spatialreference.org/ref/epsg/?page=1)
 * as illustrated below (21781 is the EPSG code for the CH1903 coordinate system. Such a code is 32767 at the maximum):
 * @code
 * COORDSYS	= PROJ4
 * COORDPARAM	= 21781
 * @endcode
 *
 */

//Adding a coordinate system is done by editing convert_to_WGS84 and convert_from_WGS84
//and implementing the XXX_to_WGS84 and WGS84_to_XXX methods.

const struct Coords::ELLIPSOID Coords::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 Equality operator that checks that two coordinates objects represent the same 3D point.
* The comparison checks either lat/lon/alt or easting/northing/alt. If both objects have nodata coordinates,
* then they are equal (even if the internal projections might be set to different systems).
* @param[in] in Coord object to compare to
* @return true or false
*/
bool Coords::operator==(const Coords& in) const {
	//check on lat/lon
	if(latitude!=IOUtils::nodata && longitude!=IOUtils::nodata) {
		const bool comparison = ( IOUtils::checkEpsilonEquality(getLat(), in.getLat(), IOUtils::lat_epsilon) &&
		                          IOUtils::checkEpsilonEquality(getLon(), in.getLon(), IOUtils::lon_epsilon) &&
		                          IOUtils::checkEpsilonEquality(getAltitude(), in.getAltitude(), IOUtils::grid_epsilon) );
		return comparison;
	}
	//check on easting/northing
	if(easting!=IOUtils::nodata && northing!=IOUtils::nodata) {
		//in this case, it means that we don't know anything about the projection parameters
		//otherwise the lat/long would have been calculated. So EPSG should be nodata
		const bool comparison = ( IOUtils::checkEpsilonEquality(getEasting(), in.getEasting(), IOUtils::grid_epsilon) &&
		                          IOUtils::checkEpsilonEquality(getNorthing(), in.getNorthing(), IOUtils::grid_epsilon) &&
		                          IOUtils::checkEpsilonEquality(getAltitude(), in.getAltitude(), IOUtils::grid_epsilon) &&
		                          getEPSG()==in.getEPSG());
		return comparison;
	}
	if(grid_i!=IOUtils::nodata && grid_j!=IOUtils::nodata && grid_k!=IOUtils::nodata) {
		//only available information is grid indices
		const bool comparison = ( grid_i==in.grid_i && grid_j==in.grid_j && grid_k==in.grid_k );
		return comparison;
	}
	//every field is nodata... the objects can only be equal if both are nodata
	if(in.isNodata()==true) return true;
	else return false;
}

/**
* @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);
}

Coords& Coords::operator=(const Coords& source) {
	if(this != &source) {
		altitude = source.altitude;
		latitude = source.latitude;
		longitude = source.longitude;
		easting = source.easting;
		northing = source.northing;
		ref_latitude = source.ref_latitude;
		ref_longitude = source.ref_longitude;
		distance_algo = source.distance_algo;
		coordsystem = source.coordsystem;
		coordparam = source.coordparam;
		grid_i = source.grid_i;
		grid_j = source.grid_j;
		grid_k = source.grid_k;
	}
	return *this;
}

bool Coords::isNodata() const {
	if( latitude==IOUtils::nodata && longitude==IOUtils::nodata &&
	    easting==IOUtils::nodata && northing==IOUtils::nodata &&
	    altitude==IOUtils::nodata &&
	    grid_i==IOUtils::nodata && grid_j==IOUtils::nodata && grid_k==IOUtils::nodata) {
		return true;
	}
	return false;
}

///< move the point by the specified distance (in m) along easting and northing
void Coords::moveByXY(const double& x_displacement, const double& y_displacement) {
	setXY(easting+x_displacement, northing+y_displacement, altitude, true);
}

///< move the point by the specified bearing and distance (in m)
void Coords::moveByBearing(const double& i_bearing, const double& i_distance) {
	double new_lat, new_lon;

	switch(distance_algo) {
		case GEO_COSINE:
			cosineInverse(latitude, longitude, i_distance, i_bearing, new_lat, new_lon);
			break;
		case GEO_VINCENTY:
			VincentyInverse(latitude, longitude, i_distance, i_bearing, new_lat, new_lon);
			break;
		default:
			throw InvalidArgumentException("Unrecognized geodesic distance algorithm selected", AT);
	}

	setLatLon(new_lat, new_lon, altitude, true);
}

/**
* @brief Simple merge strategy.
* If some fields of the first argument are empty, they will be filled by the matching field from the
* second argument.
* @param coord1 first Coords to merge, highest priority
* @param coord2 second Coords to merge, lowest priority
* @return new Coords object
*/
Coords Coords::merge(const Coords& coord1, const Coords& coord2) {
	Coords tmp(coord1);
	tmp.merge(coord2);
	return tmp;
}

/**
* @brief Simple merge strategy.
* If some fields of the current object are empty, they will be filled by the matching field from the
* provided argument.
* @param coord2 extra Coords to merge, lowest priority
*/
void Coords::merge(const Coords& coord2) {
	if(altitude==IOUtils::nodata) altitude=coord2.altitude;
	if(latitude==IOUtils::nodata) latitude=coord2.latitude;
	if(longitude==IOUtils::nodata) longitude=coord2.longitude;
	if(easting==IOUtils::nodata) easting=coord2.easting;
	if(northing==IOUtils::nodata) northing=coord2.northing;

	if(grid_i==IOUtils::nodata) grid_i=coord2.grid_i;
	if(grid_j==IOUtils::nodata) grid_j=coord2.grid_j;
	if(grid_k==IOUtils::nodata) grid_k=coord2.grid_k;

	if(ref_latitude==IOUtils::nodata) ref_latitude=coord2.ref_latitude;
	if(ref_longitude==IOUtils::nodata) ref_longitude=coord2.ref_longitude;

	if(coordsystem=="NULL") coordsystem=coord2.coordsystem;
	if(coordparam=="NULL") coordparam=coord2.coordparam;

	if(distance_algo==IOUtils::nodata) distance_algo=coord2.distance_algo;

	//in LOCAL projection, the check for the existence of the ref point will be done in the projection functions
	if(coordsystem=="LOCAL" && !coordparam.empty()) {
		parseLatLon(coordparam, ref_latitude, ref_longitude);
	}
	if(latitude!=IOUtils::nodata && coordsystem!="NULL") {
		convert_from_WGS84(latitude, longitude, easting, northing);
	}
	if(latitude==IOUtils::nodata && coordsystem!="NULL") {
		convert_to_WGS84(easting, northing, latitude, longitude);
	}
}

/**
* @brief Print the content of the Coords object (usefull for debugging)
* The Coords is bound by "<Coords>" and "</Coords>" on separate lines
*/
const std::string Coords::toString() const {
	std::ostringstream os;
	os << "<Coords>\n";
	os << "Altitude\t" << altitude << "\n";
	os << "Lat/Long\t" << printLatLon() << "\n";
	os << "Lat/Long\t" << "(" << getLat() << " , " << getLon() << ")" << "\n";
	std::streamsize p = os.precision();
	os << "X/Y_coords\t" << std::fixed << std::setprecision(0) << "(" << getEasting() << " , " << getNorthing() << ")" << "\n";
	os << std::resetiosflags(std::ios_base::fixed|std::ios_base::floatfield);
	os.precision(p);
	os << "I/J_indices\t" << "(" << getGridI() << " , " << getGridJ() << ")" << "\n";
	os << "Projection\t" << coordsystem << " " << coordparam << "\n";
	os << "EPSG\t\t" << getEPSG() << "\n";
	os << "</Coords>\n";
	return os.str();
}

std::iostream& operator<<(std::iostream& os, const Coords& coord) {
	os.write(reinterpret_cast<const char*>(&coord.ref_latitude), sizeof(coord.ref_latitude));
	os.write(reinterpret_cast<const char*>(&coord.ref_longitude), sizeof(coord.ref_longitude));
	os.write(reinterpret_cast<const char*>(&coord.altitude), sizeof(coord.altitude));
	os.write(reinterpret_cast<const char*>(&coord.latitude), sizeof(coord.latitude));
	os.write(reinterpret_cast<const char*>(&coord.longitude), sizeof(coord.longitude));
	os.write(reinterpret_cast<const char*>(&coord.easting), sizeof(coord.easting));
	os.write(reinterpret_cast<const char*>(&coord.northing), sizeof(coord.northing));
	os.write(reinterpret_cast<const char*>(&coord.grid_i), sizeof(coord.grid_i));
	os.write(reinterpret_cast<const char*>(&coord.grid_j), sizeof(coord.grid_j));
	os.write(reinterpret_cast<const char*>(&coord.grid_k), sizeof(coord.grid_k));

	const size_t s_coordsystem = coord.coordsystem.size();
	os.write(reinterpret_cast<const char*>(&s_coordsystem), sizeof(size_t));
	os.write(reinterpret_cast<const char*>(&coord.coordsystem[0]), s_coordsystem*sizeof(coord.coordsystem[0]));
	const size_t s_coordparam = coord.coordparam.size();
	os.write(reinterpret_cast<const char*>(&s_coordparam), sizeof(size_t));
	os.write(reinterpret_cast<const char*>(&coord.coordparam[0]), s_coordparam*sizeof(coord.coordparam[0]));

	os.write(reinterpret_cast<const char*>(&coord.distance_algo), sizeof(coord.distance_algo));
	return os;
}

std::iostream& operator>>(std::iostream& is, Coords& coord) {
	is.read(reinterpret_cast<char*>(&coord.ref_latitude), sizeof(coord.ref_latitude));
	is.read(reinterpret_cast<char*>(&coord.ref_longitude), sizeof(coord.ref_longitude));
	is.read(reinterpret_cast<char*>(&coord.altitude), sizeof(coord.altitude));
	is.read(reinterpret_cast<char*>(&coord.latitude), sizeof(coord.latitude));
	is.read(reinterpret_cast<char*>(&coord.longitude), sizeof(coord.longitude));
	is.read(reinterpret_cast<char*>(&coord.easting), sizeof(coord.easting));
	is.read(reinterpret_cast<char*>(&coord.northing), sizeof(coord.northing));
	is.read(reinterpret_cast<char*>(&coord.grid_i), sizeof(coord.grid_i));
	is.read(reinterpret_cast<char*>(&coord.grid_j), sizeof(coord.grid_j));
	is.read(reinterpret_cast<char*>(&coord.grid_k), sizeof(coord.grid_k));

	size_t s_coordsystem, s_coordparam;
	is.read(reinterpret_cast<char*>(&s_coordsystem), sizeof(size_t));
	coord.coordsystem.resize(s_coordsystem);
	is.read(reinterpret_cast<char*>(&coord.coordsystem[0]), s_coordsystem*sizeof(coord.coordsystem[0]));
	is.read(reinterpret_cast<char*>(&s_coordparam), sizeof(size_t));
	coord.coordparam.resize(s_coordparam);
	is.read(reinterpret_cast<char*>(&coord.coordparam[0]), s_coordparam*sizeof(coord.coordparam[0]));

	is.read(reinterpret_cast<char*>(&coord.distance_algo), sizeof(coord.distance_algo));
	return is;
}

/**
* @brief Default constructor
* 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() : ref_latitude(IOUtils::nodata), ref_longitude(IOUtils::nodata),
                   altitude(IOUtils::nodata), latitude(IOUtils::nodata), longitude(IOUtils::nodata),
                   easting(IOUtils::nodata), northing(IOUtils::nodata),
                   grid_i(IOUtils::inodata), grid_j(IOUtils::inodata), grid_k(IOUtils::inodata),
                   coordsystem("NULL"), coordparam("NULL"), distance_algo(GEO_COSINE) {}

/**
* @brief Regular constructor: usually, this is the constructor to use
* @param[in] in_coordinatesystem string identifying the coordinate system to use
* @param[in] in_parameters string giving some additional parameters for the projection (optional)
*
* See setProj() for a full description of these strings
*/
Coords::Coords(const std::string& in_coordinatesystem, const std::string& in_parameters) :
                   ref_latitude(IOUtils::nodata), ref_longitude(IOUtils::nodata),
                   altitude(IOUtils::nodata), latitude(IOUtils::nodata), longitude(IOUtils::nodata),
                   easting(IOUtils::nodata), northing(IOUtils::nodata),
                   grid_i(IOUtils::inodata), grid_j(IOUtils::inodata), grid_k(IOUtils::inodata),
                   coordsystem(in_coordinatesystem), coordparam(in_parameters), distance_algo(GEO_COSINE)
{
	setProj(in_coordinatesystem, in_parameters);
}

/**
* @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
* which coordinates have to be provided here.
* @param[in] in_lat_ref latitude of the reference point
* @param[in] in_long_ref longitude of the reference point
*/
Coords::Coords(const double& in_lat_ref, const double& in_long_ref) :
                   ref_latitude(in_lat_ref), ref_longitude(in_long_ref),
                   altitude(IOUtils::nodata), latitude(IOUtils::nodata), longitude(IOUtils::nodata),
                   easting(IOUtils::nodata), northing(IOUtils::nodata),
                   grid_i(IOUtils::inodata), grid_j(IOUtils::inodata), grid_k(IOUtils::inodata),
                   coordsystem("LOCAL"), coordparam(), distance_algo(GEO_COSINE)
{
	setLocalRef(in_lat_ref, in_long_ref);
	setProj("LOCAL", "");
}

Coords::Coords(const Coords& c) : ref_latitude(c.ref_latitude), ref_longitude(c.ref_longitude),
                   altitude(c.altitude), latitude(c.latitude), longitude(c.longitude),
                   easting(c.easting), northing(c.northing),
                   grid_i(c.grid_i), grid_j(c.grid_j), grid_k(c.grid_k),
                   coordsystem(c.coordsystem), coordparam(c.coordparam), distance_algo(c.distance_algo) {}

/**
* @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
* which coordinates have to be provided here.
* @param[in] in_coordinatesystem string identifying the coordinate system to use
* @param[in] in_parameters string giving some additional parameters for the projection (empty string if not applicable)
* @param[in] coord_spec coordinate specification
*
* The coordinate specification is given as either: "easting northing epsg" or "lat lon".
*/
Coords::Coords(const std::string& in_coordinatesystem, const std::string& in_parameters, const std::string& coord_spec)
       : ref_latitude(IOUtils::nodata), ref_longitude(IOUtils::nodata),
         altitude(IOUtils::nodata), latitude(IOUtils::nodata), longitude(IOUtils::nodata),
         easting(IOUtils::nodata), northing(IOUtils::nodata),
         grid_i(IOUtils::inodata), grid_j(IOUtils::inodata), grid_k(IOUtils::inodata),
         coordsystem(in_coordinatesystem), coordparam(in_parameters), distance_algo(GEO_COSINE)
{
	std::istringstream iss(coord_spec);
	double coord1=IOUtils::nodata, coord2=IOUtils::nodata;
	int epsg=IOUtils::inodata;

	iss >> std::skipws >> coord1;
	iss >> std::skipws >> coord2;
	iss >> std::skipws >> epsg;

	if(coord1!=IOUtils::nodata && coord2!=IOUtils::nodata && epsg!=IOUtils::inodata) {
		setEPSG(epsg);
		setXY(coord1, coord2, IOUtils::nodata);
		setProj(in_coordinatesystem, in_parameters);
	} else if(coord1!=IOUtils::nodata && coord2!=IOUtils::nodata) {
		setLatLon(coord1, coord2, IOUtils::nodata);
	}
}

/**
* @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 Returns the grid index along the X axis
* @return grid index i
*/
int Coords::getGridI() const {
	return grid_i;
}

/**
* @brief Returns the grid index along the Y axis
* @return grid index j
*/
int Coords::getGridJ() const {
	return grid_j;
}

/**
* @brief Returns the grid index along the Z axis
* @return grid index k
*/
int Coords::getGridK() const {
	return grid_k;
}

/**
* @brief Returns the projection parameters
* @param[out] proj_type projection type
* @param[out] proj_args optional arguments
*/
void Coords::getProj(std::string& proj_type, std::string& proj_args) const {
	proj_type = coordsystem;
	if(coordsystem=="LOCAL") {
		std::ostringstream dms;
		dms << "(" << decimal_to_dms(ref_latitude) << " , " << decimal_to_dms(ref_longitude) << ")";
		proj_args=dms.str();
	} else {
		proj_args = coordparam;
	}
}

/**
* @brief Print a nicely formatted lat/lon in degrees, minutes, seconds
* @return lat/lon
*/
std::string Coords::printLatLon() const {
	std::ostringstream 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] in_coordinates string containing the lat/lon to read
* @param[in] in_altitude altitude to set (optional)
* @param[in] in_update should the easting/northing be updated? (default=true)
*/
void Coords::setLatLon(const std::string& in_coordinates, const double in_altitude, const bool in_update) {
	double lat, lon;
	parseLatLon(in_coordinates, lat, lon);
	setLatLon(lat, lon, in_altitude, in_update);
}

/**
* @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] in_latitude latitude to set
* @param[in] in_longitude longitude to set
* @param[in] in_altitude altitude to set
* @param[in] in_update should the easting/northing be updated? (default=true)
*/
void Coords::setLatLon(const double in_latitude, const double in_longitude, const double in_altitude, const bool in_update) {
	latitude = in_latitude;
	longitude = in_longitude;
	if(in_altitude!=IOUtils::nodata) {
		altitude = in_altitude;
	}
	if(coordsystem!="NULL" && in_update==true) {
		convert_from_WGS84(latitude, longitude, easting, northing);
	}
	grid_i = grid_j = grid_k = IOUtils::inodata;
}

/**
* @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] in_easting easting to set
* @param[in] in_northing northing to set
* @param[in] in_altitude altitude to set
* @param[in] in_update should the easting/northing be updated? (default=true)
*/
void Coords::setXY(const double in_easting, const double in_northing, const double in_altitude, const bool in_update) {
	easting = in_easting;
	northing = in_northing;
	if(in_altitude!=IOUtils::nodata) {
		altitude = in_altitude;
	}
	if(coordsystem!="NULL" && in_update==true) {
		convert_to_WGS84(easting, northing, latitude, longitude);
	}
	grid_i = grid_j = grid_k = IOUtils::inodata;
}

/**
* @brief Set grid indices
* This index represent the position in a cartesian grid. It can not be automatically matched with
* a set of geographic coordinates because it needs the information about the said grid.
* Therefore, the coordinate object needs to be given to a grid object that will either set (i,j) or
* (lat,lon)/(easting,northing) as well as the grid's matching altitude if it is not already set.
* Any subsequent change of either (lat,lon) or (easting,northing) will reset these indexes to IOUtils::inodata.
* By default, setting (i,j) will invalidate (ie: delete) ALL geographic coordinates for the object (since we can not
* convert from grid indices to/from geographic coordinates in the current object by lack of information).
* Finally, the given indices are <b>NOT checked for validity</b>: such check must be done
* by calling Grid2DObject::gridify or Grid3DObject::gridify .
*
* @note To make it short: <b>use this method with caution!!</b>
* @param[in] in_grid_i grid index along the X direction
* @param[in] in_grid_j grid index along the Y direction
* @param[in] in_grid_k grid index along the Z direction
* @param[in] in_invalidate should the geographic coordinates be invalidated? (default=true, this flag should be false ONLY when calling from Grid2/3DObject)
*/
void Coords::setGridIndex(const int in_grid_i, const int in_grid_j, const int in_grid_k, const bool in_invalidate) {
	grid_i = in_grid_i;
	grid_j = in_grid_j;
	grid_k = in_grid_k;
	if(in_invalidate==true) {
		latitude = IOUtils::nodata;
		longitude = IOUtils::nodata;
		easting = IOUtils::nodata;
		northing = IOUtils::nodata;
		altitude = IOUtils::nodata;
	}
}

/**
* @brief Set altitude at a given value.
* If the i,j,k indices were set, reset them to inodata,
* except if specified otherwise with in_update=false.
* @param[in] in_altitude altitude above sea level, in meters
* @param[in] in_update should the indices be (if necessary) recalculated? (default=true)
*/
void Coords::setAltitude(const double in_altitude, const bool in_update) {
	altitude = in_altitude;
	if(in_update==true) {
		grid_i = grid_j = grid_k = IOUtils::inodata;
	}
}

/**
* @brief Set projection to use
* This projection will be used for converting between lat/lon and East/North
* @param[in] in_coordinatesystem string identifying the coordinate system to use
* @param[in] in_parameters string giving some additional parameters for the projection (optional)
*
*  \anchor Coordinate_types
* The coordinate system can be any of the following:
* - CH1903 for coordinates in the Swiss Grid [ref: http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf]
* - UTM for UTM coordinates (the zone must be specified in the parameters, for example 31T) [ref: http://www.oc.nps.edu/oc2902w/maps/utmups.pdf]
* - PROJ4 for coordinate conversion relying on the Proj4 library [ref: http://trac.osgeo.org/proj/]
* - LOCAL for local coordinate system (using the horizontal and vertical distance from a reference point, see Coords::geo_distances for the available choice of distance algorithms)
*/
void Coords::setProj(const std::string& in_coordinatesystem, const std::string& in_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") && ((latitude==IOUtils::nodata) || (longitude==IOUtils::nodata))) {
		convert_to_WGS84(easting, northing, latitude, longitude);
	}

	if(in_coordinatesystem.empty()) {
		coordsystem = std::string("NULL");
	} else {
		coordsystem = in_coordinatesystem;
	}
	coordparam  = in_parameters;
	if(coordsystem=="LOCAL" && !coordparam.empty()) {
		parseLatLon(coordparam, ref_latitude, ref_longitude);
	}

	//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);
	}
	//if we only had x/y but not even a coord system, we could not compute lat/long. We now do it
	if( (latitude==IOUtils::nodata || longitude==IOUtils::nodata) &&
	    (easting!=IOUtils::nodata && northing!=IOUtils::nodata) &&
	    (coordsystem != "NULL") ) {
		convert_to_WGS84(easting, northing, latitude, longitude);
	}
}

/**
* @brief Set the local projection reference coordinates
* This projection will be used for converting between lat/lon and East/North
* @param[in] in_ref_latitude latitude of the local origin
* @param[in] in_ref_longitude longitude of the local origin
*/
void Coords::setLocalRef(const double in_ref_latitude, const double in_ref_longitude) {
	if(in_ref_latitude==IOUtils::nodata || in_ref_longitude==IOUtils::nodata) {
		throw InvalidArgumentException("For LOCAL projection, please provide both reference latitude and longitude!", AT);
	}
	ref_latitude = in_ref_latitude;
	ref_longitude = in_ref_longitude;
	if(coordsystem=="LOCAL") {
		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] in_coordparam string containing the (lat,lon) of the local origin
*/
void Coords::setLocalRef(const std::string in_coordparam) {
	coordparam = in_coordparam;
	parseLatLon(coordparam, ref_latitude, ref_longitude);
	if(coordsystem=="LOCAL") {
		convert_from_WGS84(latitude, longitude, easting, northing);
	}
}

/**
* @brief Set the algorithm to use to compute distances
* Various algorithm exist that offer various precision/complexity tradeoffs.
* @param[in] in_algo enum giving the algorithm to be used (see documentation for geo_distances)
*/
void Coords::setDistances(const geo_distances in_algo) {
	distance_algo = in_algo;
	if(coordsystem=="LOCAL") {
		convert_from_WGS84(latitude, longitude, easting, northing);
	}
}

/**
* @brief Check consistency of coordinates
* When both latitude/longitude and easting/northing are given, there is
* a risk of inconsistency between these two sets of coordinates.
* This method checks that enough information is available (ie: at least one set
* of coordinates is present) and if more than one is present, that it is consistent (within 5 meters)
* It throws and exception if something is not right.
*/
void Coords::check()
{
	//calculate/check coordinates if necessary
	if(coordsystem=="LOCAL" && (ref_latitude==IOUtils::nodata || ref_longitude==IOUtils::nodata)) {
		throw InvalidArgumentException("please define a reference point for LOCAL coordinate system", AT);
	}

	if(latitude==IOUtils::nodata || longitude==IOUtils::nodata) {
		if(easting==IOUtils::nodata || northing==IOUtils::nodata) {
			throw InvalidArgumentException("missing positional parameters (easting,northing) or (lat,long) for coordinate", AT);
		}
		convert_to_WGS84(easting, northing, latitude, longitude);
	} else {
		if(easting==IOUtils::nodata || northing==IOUtils::nodata) {
			convert_from_WGS84(latitude, longitude, easting, northing);
		} else {
			double tmp_lat, tmp_lon;
			convert_to_WGS84(easting, northing, tmp_lat, tmp_lon);

			if(!IOUtils::checkEpsilonEquality(latitude, tmp_lat, IOUtils::lat_epsilon) || !IOUtils::checkEpsilonEquality(longitude, tmp_lon, IOUtils::lon_epsilon)) {
				throw InvalidArgumentException("Latitude/longitude and xllcorner/yllcorner don't match for coordinate", AT);
			}
		}
	}
}

/**
* @brief Calculate the distance between two points
* @param destination destination coordinate
* @return distance in meters
*/
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 i_update should the necessary coordinates be updated? (default=true)
*/
void Coords::copyProj(const Coords& source, const bool i_update) {
	if(!isSameProj(source)) {
		//we only do a copy if we are not already using the same projection
		if(source.coordsystem=="LOCAL") {
			coordsystem="LOCAL";
			coordparam=source.coordparam;
			ref_latitude=source.ref_latitude;
			ref_longitude=source.ref_longitude;
		} else {
			coordsystem=source.coordsystem;
			coordparam=source.coordparam;
		}

		if(i_update==true) {
			if((latitude!=IOUtils::nodata) && (longitude!=IOUtils::nodata)) {
				convert_from_WGS84(latitude, longitude, easting, northing);
			} else {
				convert_to_WGS84(easting, northing, latitude, longitude);
			}
		}
	}
}

/**
* @brief returns the epsg code of the current projection
* @return epsg code
*/
short int Coords::getEPSG() const {
	if(coordsystem=="CH1903") return 21781;
	if(coordsystem=="UTM") {
		//UTM Zone information
		short int zoneNumber;
		char zoneLetter;
		parseUTMZone(coordparam, zoneLetter, zoneNumber);
		if(zoneLetter >= 'N') {
			//northern hemisphere. We KNOW it will fit in a short int
			return (static_cast<short int>(32600+zoneNumber));
		} else {
			//southern hemisphere. We KNOW it will fit in a short int
			return (static_cast<short int>(32700+zoneNumber));
		}
	}
	if(coordsystem=="UPS") {
		//UPS Zone
		if(coordparam == "S") {
			//southern hemisphere
			return (32761);
		} else {
			//northern hemisphere
			return (32661);
		}
	}
	if(coordsystem=="PROJ4") {
		const int tmp = atoi(coordparam.c_str());
		if(tmp<0 || tmp>32767) {
			std::ostringstream ss;
			ss << "Invalid EPSG code argument: " << tmp << ". It should be between 0 and 32767! (please check EPSG registry)";
			throw InvalidArgumentException(ss.str(), AT);
		}
		return static_cast<short>(tmp);
	}

	//all others have no associated EPSG code
	return -1;
}

/**
* @brief set the current projection to a given EPSG-defined projection
* @param epsg epsg code
*/
void Coords::setEPSG(const int epsg) {
//TODO: get rid of the zone letter. This is not part of the standard and redundant (and messy)
	bool found=false;
	std::string coord_sys, coord_param;

	if(epsg<0 || epsg>32767) {
		std::ostringstream ss;
		ss << "Invalid epsg code " << epsg << " (it should be between 0 and 32767)!";
		throw InvalidArgumentException(ss.str(), AT);
	}

	if(!found && (epsg==21781)) {
		coord_sys.assign("CH1903");
		coord_param.clear();
		found=true;
	}
	if(!found && (epsg>=32601) && (epsg<=32660)) {
		//northern hemisphere
		coord_sys.assign("UTM");
		const int zoneNumber = epsg-32600;
		std::ostringstream osstream;
		osstream << zoneNumber << "N";
		coord_param=osstream.str();
		found=true;
	}
	if(!found && (epsg>=32701) && (epsg<=32760)) {
		//southern hemisphere
		coord_sys.assign("UTM");
		const int zoneNumber = epsg-32700;
		std::ostringstream osstream;
		osstream << zoneNumber << "M";
		coord_param=osstream.str();
		found=true;
	}
	if(!found && (epsg==32761 || epsg==32661)) {
		coord_sys.assign("UPS");
		coord_param = (epsg==32761)? "S" : "N";
		found=true;
	}
	if(!found) {
		//anything else has to be processed by proj4
		coord_sys.assign("PROJ4");
		std::ostringstream osstream;
		osstream << epsg;
		coord_param=osstream.str();
	}
	setProj(coord_sys, coord_param);
}

/////////////////////////////////////////////////////private methods
/**
* @brief Method converting towards WGS84
* @param[in] i_easting easting of the coordinate to convert
* @param[in] i_northing northing of the coordinate to convert
* @param[out] o_latitude converted latitude
* @param[out] o_longitude converted longitude
*/
void Coords::convert_to_WGS84(double i_easting, double i_northing, double& o_latitude, double& o_longitude) const
{
	if((i_easting!=IOUtils::nodata) && (i_northing!=IOUtils::nodata)) {
		if(coordsystem=="UTM") UTM_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
		else if(coordsystem=="UPS") UPS_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
		else if(coordsystem=="CH1903") CH1903_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
		else if(coordsystem=="LOCAL") local_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
		else if(coordsystem=="PROJ4") PROJ4_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
		else if(coordsystem=="NULL") NULL_to_WGS84(i_easting, i_northing, o_latitude, o_longitude);
		else throw UnknownValueException("Unknown coordinate system \""+coordsystem+"\"", AT);
	} else {
		o_latitude = IOUtils::nodata;
		o_longitude = IOUtils::nodata;
	}
}

/**
* @brief Method converting towards WGS84
* @param[in] i_latitude latitude of the coordinate to convert
* @param[in] i_longitude longitude of the coordinate to convert
* @param[out] o_easting converted easting
* @param[out] o_northing converted northing
*/
void Coords::convert_from_WGS84(double i_latitude, double i_longitude, double& o_easting, double& o_northing) const
{
	if((i_latitude!=IOUtils::nodata) && (i_longitude!=IOUtils::nodata)) {
		if(coordsystem=="UTM") WGS84_to_UTM(i_latitude, i_longitude, o_easting, o_northing);
		else if(coordsystem=="UPS") WGS84_to_UPS(i_latitude, i_longitude, o_easting, o_northing);
		else if(coordsystem=="CH1903") WGS84_to_CH1903(i_latitude, i_longitude, o_easting, o_northing);
		else if(coordsystem=="LOCAL") WGS84_to_local(i_latitude, i_longitude, o_easting, o_northing);
		else if(coordsystem=="PROJ4") WGS84_to_PROJ4(i_latitude, i_longitude, o_easting, o_northing);
		else if(coordsystem=="NULL") WGS84_to_NULL(i_latitude, i_longitude, o_easting, o_northing);
		else throw UnknownValueException("Unknown coordinate system \""+coordsystem+"\"", AT);
	} else {
		o_easting = IOUtils::nodata;
		o_northing = IOUtils::nodata;
	}
}

/**
* @brief Parse a latitude or longitude
* It can be formatted as any of the following examples:
* - 46° 48' 03" (with or without spaces, decimal or integer numbers)
* - 46d 48' 03" (with or without spaces, decimal or integer numbers)
* - 46 48' 03" (with spaces, decimal or integer numbers)
* - 46° 48.02'(with or without spaces, decimal or integer numbers)
* - 46d 48.02'(with or without spaces, decimal or integer numbers)
* - 46 48.02'(with spaces, decimal or integer numbers)
* - 46.802°
* - 46.802d
* - 46.802
* @param[in] dms string containing the coordinate
* @return coordinate in decimal
*/
double Coords::dms_to_decimal(const std::string& dms) {
	double d=IOUtils::nodata, m=IOUtils::nodata, s=IOUtils::nodata, decimal=IOUtils::nodata;

	if 	((sscanf(dms.c_str(), "%lf°%lf'%lf\"", &d, &m ,&s) < 3) &&
		(sscanf(dms.c_str(), "%lf° %lf' %lf\"", &d, &m ,&s) < 3) &&
		(sscanf(dms.c_str(), "%lfd%lf'%lf\"", &d, &m ,&s) < 3) &&
		(sscanf(dms.c_str(), "%lfd %lf' %lf\"", &d, &m ,&s) < 3) &&
		(sscanf(dms.c_str(), "%lf %lf' %lf\"", &d, &m ,&s) < 3) &&
		(sscanf(dms.c_str(), "%lf° %lf'", &d, &m) < 2) &&
		(sscanf(dms.c_str(), "%lf°%lf'", &d, &m) < 2) &&
		(sscanf(dms.c_str(), "%lfd %lf'", &d, &m) < 2) &&
		(sscanf(dms.c_str(), "%lfd%lf'", &d, &m) < 2) &&
		(sscanf(dms.c_str(), "%lf %lf'", &d, &m) < 2) &&
		(sscanf(dms.c_str(), "%lf°", &d) < 1) &&
		(sscanf(dms.c_str(), "%lfd", &d) < 1) &&
		(sscanf(dms.c_str(), "%lf", &d) < 1)) {
			throw InvalidFormatException("Can not parse given latitude or longitude: "+dms,AT);
	}

	decimal = fabs(d);
	if(m!=IOUtils::nodata) {
		decimal += m/60.;
	}
	if(s!=IOUtils::nodata) {
		decimal += s/3600.;
	}

	if(d<0.) return (-decimal);
	else return decimal;
}

/**
* @brief Parse a latitude-longitude pair
* It can be formatted as any of the following examples:
* - lat lon (without any spaces in the latitude or longitude string)
* - lat/lon
* - (lat;lon)
* - (lat,lon)
* @param[in] coordinates string containing the coordinates
* @param[out] lat parsed latitude
* @param[out] lon parsed longitude
*/
void Coords::parseLatLon(const std::string& coordinates, double&
lat, double& lon)
{
	const size_t len=64;
	char lat_str[len]=""; //each string must be able to accomodate the whole length to avoid buffer overflow
	char lon_str[len]="";

	if(coordinates.size()>len) {
		throw InvalidFormatException("Given lat/lon string is too long! ",AT);
	}

	if     ((sscanf(coordinates.c_str(), "%[0-9.,°d'\"-] %[0-9.,°d'\"-]", lat_str, lon_str) < 2) &&
		(sscanf(coordinates.c_str(), "%[0-9.,°d'\"- ]/%[0-9.,°d'\"- ]", lat_str, lon_str) < 2) &&
		(sscanf(coordinates.c_str(), "(%[0-9.,°d'\"- ];%[0-9.,°d'\"- ])", lat_str, lon_str) < 2) &&
		(sscanf(coordinates.c_str(), "(%[0-9.°d'\"- ],%[0-9.°d'\"- ])", lat_str, lon_str) < 2)) {
			throw InvalidFormatException("Can not parse given lat/lon: "+coordinates,AT);
	}

	lat = dms_to_decimal(std::string(lat_str));
	lon = dms_to_decimal(std::string(lon_str));
}

/**
* @brief Converts a decimal latitude or longitude to degrees, minutes, seconds
* It formats its arguments as in the following example: 46°48'03"
* @param[in] decimal decimal coordinate to convert
* @return string containing the formatted coordinate
*/
std::string Coords::decimal_to_dms(const double& decimal) {
	std::ostringstream dms;
	const double abs_dec = fabs(decimal);
	const int d = (int)floor(abs_dec);
	const int m = (int)floor( (abs_dec - (double)d)*60. );
	const double s = 3600.*(abs_dec - (double)d) - 60.*(double)m;

	if(decimal<0.)
		dms << "-";
	dms << d << "°" << m << "'" << fixed << setprecision(6) << s << "\"";
	return dms.str();
}

/**
* @brief Lenght of one degree of latitude
* This returns the lenght in meters of one degree of latitude around the given latitude
* (ie: latitude-.5, latitude+.5). See https://en.wikipedia.org/wiki/Latitude#The_length_of_a_degree_of_latitude
* @param[in] latitude latitude where to perform the computation
* @return lenght of one degree of latitude
*/
double Coords::lat_degree_lenght(const double& latitude) {
	const double a = ellipsoids[E_WGS84].a; //major ellipsoid semi-axis
	const double b = ellipsoids[E_WGS84].b;	//minor ellipsoid semi-axis
	const double e2 = (a*a-b*b) / (a*a);	//ellispoid eccentricity, squared

	const double degree_length = (Cst::PI*a*(1.-e2)) / ( 180.*pow(1.-e2*Optim::pow2(sin(latitude*Cst::to_rad)), 1.5) );
	return fabs( degree_length );
}

/**
* @brief Lenght of one degree of longitude