WSL/SLF GitLab Repository

Coords.cc 54.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
/***********************************************************************************/
/*  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 <meteoio/Coords.h>

#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]
 * - 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
 *
47
48
49
 * 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):
50
51
 * @code
 * COORDSYS	= PROJ4
52
 * COORDPARAM	= 21781
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 * @endcode
 *
 */

map<std::string, convfunc> Coords::to_wgs84;
map<std::string, convfunc> Coords::from_wgs84;
const bool Coords::__init = Coords::initializeMaps();

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

bool Coords::initializeMaps() {
	//Please don't forget to mirror the keywords here in the documentation in Coords.h!!
	to_wgs84["CH1903"]   = &Coords::CH1903_to_WGS84;
	from_wgs84["CH1903"] = &Coords::WGS84_to_CH1903;
74
75
76
77
78
79
80
81
82
	to_wgs84["UTM"]      = &Coords::UTM_to_WGS84;
	from_wgs84["UTM"]    = &Coords::WGS84_to_UTM;
	to_wgs84["PROJ4"]    = &Coords::PROJ4_to_WGS84;
	from_wgs84["PROJ4"]  = &Coords::WGS84_to_PROJ4;
	to_wgs84["LOCAL"]    = &Coords::local_to_WGS84;
	from_wgs84["LOCAL"]  = &Coords::WGS84_to_local;
	to_wgs84["NULL"]     = &Coords::NULL_to_WGS84;
	from_wgs84["NULL"]   = &Coords::WGS84_to_NULL;

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
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
139
140
	return true;
}

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

	if(latitude==IOUtils::nodata || longitude==IOUtils::nodata || easting==IOUtils::nodata || northing==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;
	} else {
		const bool comparison = ( IOUtils::checkEpsilonEquality(getLat(), in.getLat(), IOUtils::lat_epsilon) &&
				IOUtils::checkEpsilonEquality(getLon(), in.getLon(), IOUtils::lon_epsilon) );
		return comparison;
	}
}

/**
* @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;
		setFunctionPointers();
	}
	return *this;
}

/**
* @brief Print the content of the Coords object (usefull for debugging)
* The Coords is bound by "<Coords>" and "</Coords>" on separate lines
*/
std::ostream& operator<<(std::ostream &os, const Coords& coord)
{
	os << "<Coords>\n";
141
	os << "Altitude\t" << coord.altitude << "\n";
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
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
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
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
	os << "Lat/Long\t" << coord.printLatLon() << "\n";
	os << "X/Y_coords\t" << "(" << coord.getEasting() << " , " << coord.getNorthing() << ")" << "\n";
	os << "I/J_indices\t" << "(" << coord.getGridI() << " , " << coord.getGridJ() << ")" << "\n";
	os << "Projection\t" << coord.coordsystem << " " << coord.coordparam << "\n";
	os << "</Coords>\n";
	return os;
}

/**
* @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() {
	setDefaultValues();
	setProj("NULL", "NULL");
}

/**
* @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)
*
* See setProj() for a full description of these strings
*/
Coords::Coords(const std::string& _coordinatesystem, const std::string& _parameters) {
	setDefaultValues();
	setProj(_coordinatesystem, _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] _lat_ref latitude of the reference point
* @param[in] _long_ref longitude of the reference point
*/
Coords::Coords(const double& _lat_ref, const double& _long_ref)
{
	setDefaultValues();
	setLocalRef(_lat_ref, _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 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::stringstream 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::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] _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 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
* @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;
	if(_altitude!=IOUtils::nodata) {
		altitude = _altitude;
	}
	if(coordsystem!="NULL" && _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] _easting easting to set
* @param[in] _northing northing to set
* @param[in] _altitude altitude to set
* @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;
	if(_altitude!=IOUtils::nodata) {
		altitude = _altitude;
	}
	if(coordsystem!="NULL" && _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] _grid_i grid index along the X direction
* @param[in] _grid_j grid index along the Y direction
* @param[in] _grid_k grid index along the Z direction
* @param[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 _grid_i, const int _grid_j, const int _grid_k, const bool _invalidate) {
//HACK TODO make grid_i,j,k friends of Grid2/3DObject -> remove _invalidate and ALWAYS invalidate
	grid_i = _grid_i;
	grid_j = _grid_j;
	grid_k = _grid_k;
	if(_invalidate==true) {
		latitude = IOUtils::nodata;
		longitude = IOUtils::nodata;
		easting = IOUtils::nodata;
		northing = IOUtils::nodata;
		altitude = IOUtils::nodata;
	}
}

/**
* @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)
*
*  \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& _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") && ((latitude==IOUtils::nodata) || (longitude==IOUtils::nodata))) {
		convert_to_WGS84(easting, northing, latitude, longitude);
	}

	if(_coordinatesystem == "") {
		coordsystem = std::string("NULL");
	} else {
		coordsystem = _coordinatesystem;
	}
	coordparam  = _parameters;
	setFunctionPointers();

	//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);
	}
	ref_latitude = _ref_latitude;
	ref_longitude = _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] _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);
	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] _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") {
		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 _update should the necessary coordinates be updated? (default=true)
*/
void Coords::copyProj(const Coords& source, const bool _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;
		}
		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);
			}
		}
	}
}

531
532
533
534
535
536
537
538
539
540
/**
* @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;
541
		parseUTMZone(coordparam, zoneLetter, zoneNumber);
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
		if(zoneLetter >= 'M') {
			//northern hemisphere
			return (32600+zoneNumber);
		} else {
			//southern hemisphere
			return (32700+zoneNumber);
		}
	}
	if(coordsystem=="PROJ4") return atoi(coordparam.c_str());

	//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 short int epsg) {
561
//TODO: get rid of the zone letter. This is not part of the standard and redundant (and messy)
562
	bool found=false;
563
	std::string coord_sys, coord_param;
564
565

	if(!found && (epsg==21781)) {
566
567
		coord_sys="CH1903";
		coord_param="";
568
569
570
571
		found=true;
	}
	if(!found && (epsg>=32601) && (epsg<=32660)) {
		//northern hemisphere
572
		coord_sys="UTM";
573
574
575
		const short int zoneNumber = epsg-32600;
		std::ostringstream osstream;
		osstream << zoneNumber << "P";
576
		coord_param=osstream.str();
577
578
579
580
		found=true;
	}
	if(!found && (epsg>=32701) && (epsg<=32760)) {
		//southern hemisphere
581
		coord_sys="UTM";
582
583
584
		const short int zoneNumber = epsg-32700;
		std::ostringstream osstream;
		osstream << zoneNumber << "N";
585
		coord_param=osstream.str();
586
587
588
589
		found=true;
	}
	if(!found) {
		//anything else has to be processed by proj4
590
		coord_sys="PROJ4";
591
592
		std::ostringstream osstream;
		osstream << epsg;
593
		coord_param=osstream.str();
594
	}
595
	setProj(coord_sys, coord_param);
596
597
}

598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
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
717
718
/////////////////////////////////////////////////////private methods
/**
* @brief Method converting towards WGS84
* @param[in] easting easting of the coordinate to convert
* @param[in] northing northing of the coordinate to convert
* @param[out] latitude converted latitude
* @param[out] longitude converted longitude
*/
void Coords::convert_to_WGS84(double easting, double northing, double& latitude, double& longitude) const
{
	if((easting!=IOUtils::nodata) && (northing!=IOUtils::nodata)) {
		(this->*convToWGS84)(easting, northing, latitude, longitude);
	} else {
		latitude = IOUtils::nodata;
		longitude = IOUtils::nodata;
	}
}

/**
* @brief Method converting towards WGS84
* @param[in] latitude latitude of the coordinate to convert
* @param[in] longitude longitude of the coordinate to convert
* @param[out] easting converted easting
* @param[out] northing converted northing
*/
void Coords::convert_from_WGS84(double latitude, double longitude, double& easting, double& northing) const
{
	if((latitude!=IOUtils::nodata) && (longitude!=IOUtils::nodata)) {
		(this->*convFromWGS84)(latitude, longitude, easting, northing);
	} else {
		easting = IOUtils::nodata;
		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 = d;
	if(m!=IOUtils::nodata) {
		decimal += m/60.;
	}
	if(s!=IOUtils::nodata) {
		decimal += s/3600.;
	}

	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) {
	char lat_str[32]="";
	char lon_str[32]="";

	if(coordinates.size()>(32+32)) {
		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::stringstream dms;
	const int d = (int)floor(decimal);
	const int m = (int)floor( (decimal - (double)d)*60. );
719
	const double s = 3600.*(decimal - (double)d) - 60.*(double)m;
720

721
	dms << d << "°" << m << "'" << fixed << setprecision(6) << s << "\"";
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
	return dms.str();
}

/**
* @brief Coordinate conversion: from WGS84 Lat/Long to Swiss grid
* See http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf for more.
* @param[in] lat_in Decimal Latitude
* @param[in] long_in Decimal Longitude
* @param[out] east_out easting coordinate (Swiss system)
* @param[out] north_out northing coordinate (Swiss system)
*/
void Coords::WGS84_to_CH1903(double lat_in, double long_in, double& east_out, double& north_out) const
{
	//converts WGS84 coordinates (lat,long) to the Swiss coordinates. See http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf
	//The elevation is supposed to be above sea level, so it does not require any conversion
	//lat and long must be decimal (and they will be converted to seconds)
	const double phi_p = (lat_in*3600. - 169028.66) / 10000.;
	const double lambda_p = (long_in*3600. - 26782.5) / 10000.;

	east_out = 600072.37
		+ 211455.93	* lambda_p
		- 10938.51	* lambda_p * phi_p
		- 0.36		* lambda_p * (phi_p*phi_p)
		- 44.54		* (lambda_p*lambda_p*lambda_p);

	north_out = 200147.07
		+ 308807.95	* phi_p
		+ 3745.25	* (lambda_p*lambda_p)
		+ 76.63		* (phi_p*phi_p)
		- 194.56	* (lambda_p*lambda_p) * phi_p
		+ 119.79	* (phi_p*phi_p*phi_p);

	/*// if necessary for the elevation, uncomment this block
	h_out = h_in - 49.55
		+ 2.73		* lambda_p
		+ 6.94		* phi_p;
	*/
}

/**
* @brief Coordinate conversion: from Swiss grid to WGS84 Lat/Long
* See http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf for more.
* @param[in] east_in easting coordinate (Swiss system)
* @param[in] north_in northing coordinate (Swiss system)
* @param[out] lat_out Decimal Latitude
* @param[out] long_out Decimal Longitude
*/
void Coords::CH1903_to_WGS84(double east_in, double north_in, double& lat_out, double& long_out) const
{
	//converts Swiss coordinates to WGS84 coordinates (lat,long). See http://geomatics.ladetto.ch/ch1903_wgs84_de.pdf
	//The elevation is supposed to be above sea level, so it does not require any conversion
	//lat and long are decimal
	const double y_p = (east_in - 600000.) / 1000000.;
	const double x_p = (north_in - 200000.) / 1000000.;

	const double lambda_p = 2.6779094
		+ 4.728982	* y_p
		+ 0.791484	* y_p * x_p
		+ 0.1306	* y_p * (x_p*x_p)
		- 0.0436	* (y_p*y_p*y_p);

	const double phi_p = 16.9023892
		+ 3.238272	* x_p
		- 0.270978	* (y_p*y_p)
		- 0.002528	* (x_p*x_p)
		- 0.0447	* (y_p*y_p) * x_p
		- 0.0140	* (x_p*x_p*x_p);

	lat_out = phi_p * 100./36.;
	long_out = lambda_p * 100./36.;

	/*// if necessary for the elevation, uncomment this block
	h_out = h_in + 49.55
		- 12.60		* y_p
		- 22.64		* x_p;
	*/
}

int Coords::getUTMZone(const double _latitude, const double _longitude, std::string& zone_out) const
{//This routine determines the correct UTM letter designator for the given latitude
//UTM limits its coverage to [80S , 84N], outside of this, returns Y/Z/A/B for the zone

	//computing zone number, assuming longitude in [-180. ; 180[
	int ZoneNumber = int((_longitude + 180.)/6.) + 1;

	// Special zones for Scandinavia
	if( _latitude >= 72.0 && _latitude < 84.0 ) {
		if(      _longitude >= 0.0  && _longitude <  9.0 ) ZoneNumber = 31;
		else if( _longitude >= 9.0  && _longitude < 21.0 ) ZoneNumber = 33;
		else if( _longitude >= 21.0 && _longitude < 33.0 ) ZoneNumber = 35;
		else if( _longitude >= 33.0 && _longitude < 42.0 ) ZoneNumber = 37;
	 }
	if( latitude >= 56.0 && _latitude < 64.0 && _longitude >= 3.0 && _longitude < 12.0 ) {
		ZoneNumber = 32;
	}

	//getting zone letter
	char zoneLetter='Z';
	if     ((0 >= _longitude) && (_latitude >  84)) zoneLetter = 'Y';
	else if((0 <  _longitude) && (_latitude >  84)) zoneLetter = 'Z';
	else if((84 >= _latitude) && (_latitude >= 72)) zoneLetter = 'X';
	else if((72 > _latitude) && (_latitude >= 64)) zoneLetter = 'W';
	else if((64 > _latitude) && (_latitude >= 56)) zoneLetter = 'V';
	else if((56 > _latitude) && (_latitude >= 48)) zoneLetter = 'U';
	else if((48 > _latitude) && (_latitude >= 40)) zoneLetter = 'T';
	else if((40 > _latitude) && (_latitude >= 32)) zoneLetter = 'S';
	else if((32 > _latitude) && (_latitude >= 24)) zoneLetter = 'R';
	else if((24 > _latitude) && (_latitude >= 16)) zoneLetter = 'Q';
	else if((16 > _latitude) && (_latitude >= 8)) zoneLetter = 'P';
	else if(( 8 > _latitude) && (_latitude >= 0)) zoneLetter = 'N';
	else if(( 0 > _latitude) && (_latitude >= -8)) zoneLetter = 'M';
	else if((-8 > _latitude) && (_latitude >= -16)) zoneLetter = 'L';
	else if((-16 > _latitude) && (_latitude >= -24)) zoneLetter = 'K';
	else if((-24 > _latitude) && (_latitude >= -32)) zoneLetter = 'J';
	else if((-32 > _latitude) && (_latitude >= -40)) zoneLetter = 'H';
	else if((-40 > _latitude) && (_latitude >= -48)) zoneLetter = 'G';
	else if((-48 > _latitude) && (_latitude >= -56)) zoneLetter = 'F';
	else if((-56 > _latitude) && (_latitude >= -64)) zoneLetter = 'E';
	else if((-64 > _latitude) && (_latitude >= -72)) zoneLetter = 'D';
	else if((-72 > _latitude) && (_latitude >= -80)) zoneLetter = 'C';
	else if((0 >= _longitude) && (latitude <= -80)) zoneLetter = 'A';
	else if((0 <  _longitude) && (latitude <= -80)) zoneLetter = 'B';

	std::stringstream zone;
	zone << ZoneNumber << zoneLetter;
	zone_out = zone.str();

	return ZoneNumber;
}

/**
* @brief Coordinate conversion: from WGS84 Lat/Long to UTM grid
* See http://www.oc.nps.edu/oc2902w/maps/utmups.pdf for more.
* @param[in] lat_in Decimal Latitude
* @param[in] long_in Decimal Longitude
* @param[out] east_out easting coordinate (Swiss system)
* @param[out] north_out northing coordinate (Swiss system)
*/
void Coords::WGS84_to_UTM(double lat_in, double long_in, double& east_out, double& north_out) const
{//converts WGS84 coordinates (lat,long) to UTM coordinates.
//See USGS Bulletin 1532 or http://earth-info.nga.mil/GandG/publications/tm8358.2/TM8358_2.pdf
//also http://www.uwgb.edu/dutchs/usefuldata/UTMFormulas.HTM
//also http://www.oc.nps.edu/oc2902w/maps/utmups.pdf or Chuck Gantz (http://www.gpsy.com/gpsinfo/geotoutm/)
	//Geometric constants
	const double to_rad = M_PI / 180.0;
	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 eP2 = e2 / (1.-e2);	//second ellispoid eccentricity, squared (=(a²-b²)/b²)
	const double k0 = 0.9996;	//scale factor for the projection

	//getting posistion parameters
	std::string zone;
	long_in = fmod(long_in+360.+180., 360.) - 180.; //normalized to [-180. ; 180.[
	const double Long = long_in * to_rad;
	const double Lat = lat_in * to_rad;
878
879
880
881
882
883
884
885
886
	int zoneNumber = getUTMZone(lat_in, long_in, zone);
	short int in_zoneNumber;
	char in_zoneLetter;
	parseUTMZone(coordparam, in_zoneLetter, in_zoneNumber);
	if(in_zoneNumber!=zoneNumber) {
		std::cerr << "[W] requested UTM zone is not appropriate for the given coordinates. Normally, It should be zone ";
		std::cerr << zoneNumber << "\n";
		zoneNumber = in_zoneNumber;
	}
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
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
	const double long0 = (double)((zoneNumber - 1)*6 - 180 + 3) * to_rad; //+3 puts origin in middle of zone

	//Geometrical parameters
	const double nu = a / sqrt(1.-e2*IOUtils::pow2(sin(Lat))); //radius of curvature of the earth perpendicular to the meridian plane
	const double p = (Long-long0);

	//calculating first the coefficients of the series, then the Meridional Arc M itself
	const double n = (a-b)/(a+b);
	const double n2=n*n, n3=n*n*n, n4=n*n*n*n, n5=n*n*n*n*n, n6=n*n*n*n*n*n;
	const double A = a           * (1. - n + 5./4.*(n2 - n3) + 81./64.*(n4 - n5));
	const double B = (3./2.*a)   * (n - n2 + 7./8.*(n3 - n4) + 55./64.*(n5 - n6));
	const double C = (15./16.*a) * (n2 - n3 + 3./4.*(n4 - n5));
	const double D = (35./48.*a) * (n3 - n4 + 11./16.*(n5 - n6));
	const double E = (315./512.*a) * (n4 - n5); //correctrion of ~0.03mm
	const double M = A*Lat - B*sin(2.*Lat) + C*sin(4.*Lat) - D*sin(6.*Lat) + E*sin(8.*Lat);

	//calculating the coefficients for the series
	const double K1 = M*k0;
	const double K2 = 1./4.*k0*nu*sin(2.*Lat);
	const double K3 = (k0*nu*sin(Lat)*IOUtils::pow3(cos(Lat))*1./24.) * (5. - IOUtils::pow2(tan(Lat)) + 9.*eP2*IOUtils::pow2(cos(Lat)) + 4.*eP2*eP2*IOUtils::pow4(cos(Lat)));
	const double K4 = k0*nu*cos(Lat);
	const double K5 = (k0*nu*IOUtils::pow3(cos(Lat))*1./6.) * (1. - IOUtils::pow2(tan(Lat)) + eP2*IOUtils::pow2(cos(Lat)));

	north_out = K1 + K2*p*p + K3*p*p*p*p;
	east_out = K4*p + K5*p*p*p + 500000.0;

	if(Lat < 0)
		north_out += 10000000.0; //offset for southern hemisphere
}

/**
* @brief Coordinate conversion: from UTM grid to WGS84 Lat/Long
* See http://www.oc.nps.edu/oc2902w/maps/utmups.pdf for more.
* @param[in] east_in easting coordinate (UTM)
* @param[in] north_in northing coordinate (UTM)
* @param[out] lat_out Decimal Latitude
* @param[out] long_out Decimal Longitude
*/
void Coords::UTM_to_WGS84(double east_in, double north_in, double& lat_out, double& long_out) const
{//converts UTM coordinates to WGS84 coordinates (lat,long).
//See USGS Bulletin 1532 or http://earth-info.nga.mil/GandG/publications/tm8358.2/TM8358_2.pdf
//also http://www.uwgb.edu/dutchs/usefuldata/UTMFormulas.HTM
//also http://www.oc.nps.edu/oc2902w/maps/utmups.pdf or Chuck Gantz (http://www.gpsy.com/gpsinfo/geotoutm/)
	//Geometric constants
	const double to_deg = 180.0 / M_PI;
	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 eP2 = e2 / (1.-e2);	//second ellispoid eccentricity, squared (=(a²-b²)/b²)
	const double k0 = 0.9996;		//scale factor for the projection

	//UTM Zone information
939
	short int zoneNumber;
940
	char zoneLetter;
941
	parseUTMZone(coordparam, zoneLetter, zoneNumber);
942
943
944

	//set reference parameters: central meridian of the zone, true northing and easting
	//please note that the special zones still use the reference meridian as given by their zone number (ie: even if it might not be central anymore)
945
	const int long0 = ((int)zoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in "middle" of zone as required for the projection meridian (might not be the middle for special zones)
946
947
948
949
950
	if(zoneLetter<='N') {
		north_in -= 10000000.0; //offset used for southern hemisphere
	}
	east_in -= 500000.0; //longitude offset: x coordinate is relative to central meridian

951
	//calculating footprint latitude fp (it should be done using a few iterations)
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
	const double arc = north_in/k0; //Meridional arc
	const double mu = arc / (a*(1.-e2/4.-3.*e2*e2/64.-5.*e2*e2*e2/256.));
	const double e1 = (1.-b/a) / (1.+b/a); //simplification of [1 - (1 - e2)1/2]/[1 + (1 - e2)1/2]
	const double J1 = (3./2.*e1 - 27./32.*e1*e1*e1);
	const double J2 = (21./16.*e1*e1 - 55./32.*e1*e1*e1*e1);
	const double J3 = (151./96.*e1*e1*e1);
	const double J4 = (1097./512.*e1*e1*e1*e1);
	const double fp = mu + J1*sin(2.*mu) + J2*sin(4.*mu) + J3*sin(6.*mu) + J4*sin(8.*mu);

	//calculating the parameters
	const double C1 = eP2 * IOUtils::pow2(cos(fp));
	const double T1 = IOUtils::pow2( tan(fp) );
	const double R1 = a*(1.-e2) / pow((1.-e2*IOUtils::pow2(sin(fp))), 1.5);
	const double N1 = a / sqrt(1.-e2*IOUtils::pow2(sin(fp)));
	const double D = east_in / (N1*k0);

	//calculating the coefficients of the series for latitude and longitude
	const double Q1 = N1*tan(fp)/R1;
	const double Q2 = 0.5*D*D;
	const double Q3 = (5. + 3.*T1 + 10.*C1 - 4.*C1*C1 - 9.*eP2) * 1./24.*D*D*D*D;
	const double Q4 = (61. + 90.*T1 + 298.*C1 + 45.*T1*T1 - 3.*C1*C1 - 252.*eP2) * 1./720.*D*D*D*D*D*D;
973
974
	//const double Q4extra = (1385. + 3633.*T1 + 4095.*T1*T1 + 1575.*T1*T1*T1) * 1./40320.*D*D*D*D*D*D*D*D;

975
976
977
	const double Q5 = D;
	const double Q6 = (1. + 2.*T1 + C1) * 1./6.*D*D*D;
	const double Q7 = (5. - 2.*C1 + 28.*T1 - 3.*C1*C1 + 8.*eP2 + 24.*T1*T1) * 1./120.*D*D*D*D*D;
978
979
980
981
982
	//const double Q7extra = (61. + 662.*T1 + 1320.*T1*T1 +720.*T1*T1*T1) * 1./5040.*D*D*D*D*D*D*D;

	lat_out = (fp - Q1 * (Q2 - Q3 + Q4 /*+Q4extra*/))*to_deg;
	long_out = (double)long0 + ((Q5 - Q6 + Q7 /*-Q7extra*/)/cos(fp))*to_deg;
}
983

984
985
986
987
988
989
990
991
992
993
994
void Coords::parseUTMZone(const std::string& zone_info, char& zoneLetter, short int& zoneNumber) const
{ //helper method: parse a UTM zone specification string into letter and number
	if ((sscanf(zone_info.c_str(), "%hd%c", &zoneNumber, &zoneLetter) < 2) &&
		(sscanf(zone_info.c_str(), "%hd %c)", &zoneNumber, &zoneLetter) < 2)) {
			throw InvalidFormatException("Can not parse given UTM zone: "+zone_info,AT);
	}
	zoneLetter = toupper(zoneLetter); //just in case... (sorry for the pun!)
	if(zoneLetter=='Y' || zoneLetter=='Z' || zoneLetter=='A' || zoneLetter=='B') {
			//Special zones for the poles: we should NOT use UTM in these regions!
			throw InvalidFormatException("Invalid UTM zone: "+zone_info+" (trying to use UTM in polar regions)",AT);
	}
995
996
997
998
999
1000
}

/**
* @brief Coordinate conversion: from WGS84 Lat/Long to proj4 parameters
* @param lat_in Decimal Latitude
* @param long_in Decimal Longitude
For faster browsing, not all history is shown. View entire blame