WSL/SLF GitLab Repository

Coords.cc 75.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
/***********************************************************************************/
/*  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/>.
*/
18
#include <cmath>
19
#include <cstdio>
20
#include <iomanip>
21

22
#include <meteoio/dataClasses/Coords.h>
23
#include <meteoio/MathOptim.h>
24
#include <meteoio/meteoLaws/Meteoconst.h> //for math constants
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

#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]
45
 * - 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]
46
47
48
49
50
51
52
53
 * - 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
 *
54
55
56
 * 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):
57
58
 * @code
 * COORDSYS	= PROJ4
59
 * COORDPARAM	= 21781
60
61
62
63
 * @endcode
 *
 */

64
65
//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.
66
67
68
69
70
71
72
73
74
75
76

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
};

/**
77
78
79
* @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).
80
81
82
83
* @param[in] in Coord object to compare to
* @return true or false
*/
bool Coords::operator==(const Coords& in) const {
84
	//check on lat/lon
85
86
	if(latitude!=IOUtils::nodata && longitude!=IOUtils::nodata) {
		const bool comparison = ( IOUtils::checkEpsilonEquality(getLat(), in.getLat(), IOUtils::lat_epsilon) &&
87
88
		                          IOUtils::checkEpsilonEquality(getLon(), in.getLon(), IOUtils::lon_epsilon) &&
		                          IOUtils::checkEpsilonEquality(getAltitude(), in.getAltitude(), IOUtils::grid_epsilon) );
89
90
		return comparison;
	}
91
	//check on easting/northing
92
	if(easting!=IOUtils::nodata && northing!=IOUtils::nodata) {
93
94
		//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
95
96
97
		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) &&
98
		                          getEPSG()==in.getEPSG());
99
100
101
		return comparison;
	}
	if(grid_i!=IOUtils::nodata && grid_j!=IOUtils::nodata && grid_k!=IOUtils::nodata) {
102
103
104
105
		//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;
	}
106
107
108
	//every field is nodata... the objects can only be equal if both are nodata
	if(in.isNodata()==true) return true;
	else return false;
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
}

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

139
140
141
142
143
144
145
146
147
148
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;
}

149
150
151
152
153
154
///< 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)
155
void Coords::moveByBearing(const double& i_bearing, const double& i_distance) {
156
157
158
159
	double new_lat, new_lon;

	switch(distance_algo) {
		case GEO_COSINE:
160
			cosineInverse(latitude, longitude, i_distance, i_bearing, new_lat, new_lon);
161
162
			break;
		case GEO_VINCENTY:
163
			VincentyInverse(latitude, longitude, i_distance, i_bearing, new_lat, new_lon);
164
165
166
167
168
169
170
171
			break;
		default:
			throw InvalidArgumentException("Unrecognized geodesic distance algorithm selected", AT);
	}

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

172
173
/**
* @brief Simple merge strategy.
174
* If some fields of the first argument are empty, they will be filled by the matching field from the
175
176
177
178
179
180
181
182
* 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);
Fierz's avatar
Fierz committed
183
	return tmp;
184
185
186
187
}

/**
* @brief Simple merge strategy.
188
* If some fields of the current object are empty, they will be filled by the matching field from the
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
* 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
212
	if(coordsystem=="LOCAL" && !coordparam.empty()) {
213
214
		parseLatLon(coordparam, ref_latitude, ref_longitude);
	}
215
216
217
218
219
220
221
222
	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);
	}
}

223
224
225
226
/**
* @brief Print the content of the Coords object (usefull for debugging)
* The Coords is bound by "<Coords>" and "</Coords>" on separate lines
*/
227
const std::string Coords::toString() const {
228
	std::ostringstream os;
229
	os << "<Coords>\n";
230
231
232
	os << "Altitude\t" << altitude << "\n";
	os << "Lat/Long\t" << printLatLon() << "\n";
	os << "Lat/Long\t" << "(" << getLat() << " , " << getLon() << ")" << "\n";
233
	std::streamsize p = os.precision();
234
	os << "X/Y_coords\t" << std::fixed << std::setprecision(0) << "(" << getEasting() << " , " << getNorthing() << ")" << "\n";
235
236
	os << std::resetiosflags(std::ios_base::fixed|std::ios_base::floatfield);
	os.precision(p);
237
238
239
	os << "I/J_indices\t" << "(" << getGridI() << " , " << getGridJ() << ")" << "\n";
	os << "Projection\t" << coordsystem << " " << coordparam << "\n";
	os << "EPSG\t\t" << getEPSG() << "\n";
240
	os << "</Coords>\n";
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
	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));
264
265
266
	return os;
}

267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
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));
288
	return is;
289
290
}

291
292
293
294
295
/**
* @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...
*/
296
297
298
299
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),
300
                   coordsystem("NULL"), coordparam("NULL"), distance_algo(GEO_COSINE) {}
301
302
303

/**
* @brief Regular constructor: usually, this is the constructor to use
304
305
* @param[in] in_coordinatesystem string identifying the coordinate system to use
* @param[in] in_parameters string giving some additional parameters for the projection (optional)
306
307
308
*
* See setProj() for a full description of these strings
*/
309
310
311
312
313
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),
314
                   coordsystem(in_coordinatesystem), coordparam(in_parameters), distance_algo(GEO_COSINE)
315
{
316
	setProj(in_coordinatesystem, in_parameters);
317
318
319
320
321
322
}

/**
* @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.
323
324
* @param[in] in_lat_ref latitude of the reference point
* @param[in] in_long_ref longitude of the reference point
325
*/
326
327
328
329
330
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),
331
                   coordsystem("LOCAL"), coordparam(), distance_algo(GEO_COSINE)
332
{
333
	setLocalRef(in_lat_ref, in_long_ref);
334
335
336
	setProj("LOCAL", "");
}

337
338
339
340
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),
341
                   coordsystem(c.coordsystem), coordparam(c.coordparam), distance_algo(c.distance_algo) {}
342

343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
/**
* @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);
	}
}

377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
/**
* @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") {
449
		std::ostringstream dms;
450
451
452
453
454
455
456
457
458
459
460
461
		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 {
462
	std::ostringstream dms;
463
464
465
466
467
468
469
470
471
472
	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().
473
474
475
* @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)
476
*/
477
void Coords::setLatLon(const std::string& in_coordinates, const double in_altitude, const bool in_update) {
478
	double lat, lon;
479
480
	parseLatLon(in_coordinates, lat, lon);
	setLatLon(lat, lon, in_altitude, in_update);
481
482
483
484
485
486
487
}

/**
* @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().
488
489
490
491
* @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)
492
*/
493
494
495
496
497
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;
498
	}
499
	if(coordsystem!="NULL" && in_update==true) {
500
501
502
503
504
505
506
507
508
509
		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().
510
511
512
513
* @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)
514
*/
515
516
517
518
519
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;
520
	}
521
	if(coordsystem!="NULL" && in_update==true) {
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
		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>
540
541
542
543
* @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)
544
*/
545
546
547
548
549
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) {
550
551
552
553
554
555
556
557
		latitude = IOUtils::nodata;
		longitude = IOUtils::nodata;
		easting = IOUtils::nodata;
		northing = IOUtils::nodata;
		altitude = IOUtils::nodata;
	}
}

558
/**
559
* @brief Set altitude at a given value.
560
561
562
563
564
565
566
567
568
569
570
571
* 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;
	}
}

572
573
574
/**
* @brief Set projection to use
* This projection will be used for converting between lat/lon and East/North
575
576
* @param[in] in_coordinatesystem string identifying the coordinate system to use
* @param[in] in_parameters string giving some additional parameters for the projection (optional)
577
578
579
580
581
582
583
584
*
*  \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)
*/
585
void Coords::setProj(const std::string& in_coordinatesystem, const std::string& in_parameters) {
586
587
588
589
590
591
	//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);
	}

592
	if(in_coordinatesystem.empty()) {
593
594
		coordsystem = std::string("NULL");
	} else {
595
		coordsystem = in_coordinatesystem;
596
	}
597
	coordparam  = in_parameters;
598
	if(coordsystem=="LOCAL" && !coordparam.empty()) {
599
600
		parseLatLon(coordparam, ref_latitude, ref_longitude);
	}
601
602
603
604
605

	//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);
	}
606
607
608
609
610
611
	//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);
	}
612
613
614
615
616
}

/**
* @brief Set the local projection reference coordinates
* This projection will be used for converting between lat/lon and East/North
617
618
* @param[in] in_ref_latitude latitude of the local origin
* @param[in] in_ref_longitude longitude of the local origin
619
*/
620
621
void Coords::setLocalRef(const double in_ref_latitude, const double in_ref_longitude) {
	if(in_ref_latitude==IOUtils::nodata || in_ref_longitude==IOUtils::nodata) {
622
623
		throw InvalidArgumentException("For LOCAL projection, please provide both reference latitude and longitude!", AT);
	}
624
625
	ref_latitude = in_ref_latitude;
	ref_longitude = in_ref_longitude;
626
627
628
629
630
631
632
633
	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
634
* @param[in] in_coordparam string containing the (lat,lon) of the local origin
635
*/
636
637
void Coords::setLocalRef(const std::string in_coordparam) {
	coordparam = in_coordparam;
638
639
640
641
642
643
644
645
646
	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.
647
* @param[in] in_algo enum giving the algorithm to be used (see documentation for geo_distances)
648
*/
649
650
void Coords::setDistances(const geo_distances in_algo) {
	distance_algo = in_algo;
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
	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
717
* @param i_update should the necessary coordinates be updated? (default=true)
718
*/
719
void Coords::copyProj(const Coords& source, const bool i_update) {
720
721
722
723
724
725
726
727
728
729
730
731
	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;
		}

732
		if(i_update==true) {
733
734
735
736
737
738
739
740
741
			if((latitude!=IOUtils::nodata) && (longitude!=IOUtils::nodata)) {
				convert_from_WGS84(latitude, longitude, easting, northing);
			} else {
				convert_to_WGS84(easting, northing, latitude, longitude);
			}
		}
	}
}

742
743
744
745
746
747
748
749
750
751
/**
* @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;
752
		parseUTMZone(coordparam, zoneLetter, zoneNumber);
753
		if(zoneLetter >= 'N') {
754
755
			//northern hemisphere. We KNOW it will fit in a short int
			return (static_cast<short int>(32600+zoneNumber));
756
		} else {
757
758
			//southern hemisphere. We KNOW it will fit in a short int
			return (static_cast<short int>(32700+zoneNumber));
759
760
		}
	}
761
762
763
764
765
766
767
768
769
770
	if(coordsystem=="UPS") {
		//UPS Zone
		if(coordparam == "S") {
			//southern hemisphere
			return (32761);
		} else {
			//northern hemisphere
			return (32661);
		}
	}
771
772
773
	if(coordsystem=="PROJ4") {
		const int tmp = atoi(coordparam.c_str());
		if(tmp<0 || tmp>32767) {
774
			std::ostringstream ss;
775
776
777
778
779
			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);
	}
780
781
782
783
784
785
786
787
788

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

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

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

800
	if(!found && (epsg==21781)) {
801
		coord_sys.assign("CH1903");
802
		coord_param.clear();
803
804
805
806
		found=true;
	}
	if(!found && (epsg>=32601) && (epsg<=32660)) {
		//northern hemisphere
807
		coord_sys.assign("UTM");
808
		const int zoneNumber = epsg-32600;
809
		std::ostringstream osstream;
810
		osstream << zoneNumber << "N";
811
		coord_param=osstream.str();
812
813
814
815
		found=true;
	}
	if(!found && (epsg>=32701) && (epsg<=32760)) {
		//southern hemisphere
816
		coord_sys.assign("UTM");
817
		const int zoneNumber = epsg-32700;
818
		std::ostringstream osstream;
819
		osstream << zoneNumber << "M";
820
		coord_param=osstream.str();
821
822
		found=true;
	}
823
	if(!found && (epsg==32761 || epsg==32661)) {
824
		coord_sys.assign("UPS");
825
826
827
		coord_param = (epsg==32761)? "S" : "N";
		found=true;
	}
828
829
	if(!found) {
		//anything else has to be processed by proj4
830
		coord_sys.assign("PROJ4");
831
832
		std::ostringstream osstream;
		osstream << epsg;
833
		coord_param=osstream.str();
834
	}
835
	setProj(coord_sys, coord_param);
836
837
}

838
839
840
/////////////////////////////////////////////////////private methods
/**
* @brief Method converting towards WGS84
841
842
843
844
* @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
845
*/
846
void Coords::convert_to_WGS84(double i_easting, double i_northing, double& o_latitude, double& o_longitude) const
847
{
848
	if((i_easting!=IOUtils::nodata) && (i_northing!=IOUtils::nodata)) {
849
850
851
852
853
854
855
		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);
856
	} else {
857
858
		o_latitude = IOUtils::nodata;
		o_longitude = IOUtils::nodata;
859
860
861
862
863
	}
}

/**
* @brief Method converting towards WGS84
864
865
866
867
* @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
868
*/
869
void Coords::convert_from_WGS84(double i_latitude, double i_longitude, double& o_easting, double& o_northing) const
870
{
871
	if((i_latitude!=IOUtils::nodata) && (i_longitude!=IOUtils::nodata)) {
872
873
874
875
876
877
878
		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);
879
	} else {
880
881
		o_easting = IOUtils::nodata;
		o_northing = IOUtils::nodata;
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
	}
}

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

919
	decimal = fabs(d);
920
921
922
923
924
925
926
	if(m!=IOUtils::nodata) {
		decimal += m/60.;
	}
	if(s!=IOUtils::nodata) {
		decimal += s/3600.;
	}

927
928
	if(d<0.) return (-decimal);
	else return decimal;
929
930
931
932
933
934
935
936
937
938
939
940
941
942
}

/**
* @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&
943
944
945
946
947
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]="";
948

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

953
	if     ((sscanf(coordinates.c_str(), "%[0-9.,°d'\"-] %[0-9.,°d'\"-]", lat_str, lon_str) < 2) &&
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
		(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) {
971
	std::ostringstream dms;
972
973
974
975
	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;
976

977
978
	if(decimal<0.)
		dms << "-";
979
	dms << d << "°" << m << "'" << fixed << setprecision(6) << s << "\"";
980
981
982
	return dms.str();
}

983
984
985
986
987
988
989
990
991
992
993
994
/**
* @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

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

/**
* @brief Lenght of one degree of longitude
For faster browsing, not all history is shown. View entire blame