GeographicLib  1.30
UTMUPS.cpp
Go to the documentation of this file.
1 /**
2  * \file UTMUPS.cpp
3  * \brief Implementation for GeographicLib::UTMUPS class
4  *
5  * Copyright (c) Charles Karney (2008-2012) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
10 #include <GeographicLib/UTMUPS.hpp>
11 #include <GeographicLib/MGRS.hpp>
12 #include <GeographicLib/PolarStereographic.hpp>
13 #include <GeographicLib/TransverseMercator.hpp>
14 #include <GeographicLib/Utility.hpp>
15 
16 namespace GeographicLib {
17 
18  using namespace std;
19 
20  const Math::real UTMUPS::falseeasting_[4] =
21  { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_,
22  MGRS::utmeasting_ * MGRS::tile_, MGRS::utmeasting_ * MGRS::tile_ };
23  const Math::real UTMUPS::falsenorthing_[4] =
24  { MGRS::upseasting_ * MGRS::tile_, MGRS::upseasting_ * MGRS::tile_,
25  MGRS::maxutmSrow_ * MGRS::tile_, MGRS::minutmNrow_ * MGRS::tile_ };
26  const Math::real UTMUPS::mineasting_[4] =
27  { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_,
28  MGRS::minutmcol_ * MGRS::tile_, MGRS::minutmcol_ * MGRS::tile_ };
29  const Math::real UTMUPS::maxeasting_[4] =
30  { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_,
31  MGRS::maxutmcol_ * MGRS::tile_, MGRS::maxutmcol_ * MGRS::tile_ };
32  const Math::real UTMUPS::minnorthing_[4] =
33  { MGRS::minupsSind_ * MGRS::tile_, MGRS::minupsNind_ * MGRS::tile_,
34  MGRS::minutmSrow_ * MGRS::tile_,
35  (MGRS::minutmNrow_ + MGRS::minutmSrow_ - MGRS::maxutmSrow_)
36  * MGRS::tile_ };
37  const Math::real UTMUPS::maxnorthing_[4] =
38  { MGRS::maxupsSind_ * MGRS::tile_, MGRS::maxupsNind_ * MGRS::tile_,
39  (MGRS::maxutmSrow_ + MGRS::maxutmNrow_ - MGRS::minutmNrow_) * MGRS::tile_,
40  MGRS::maxutmNrow_ * MGRS::tile_ };
41 
42  int UTMUPS::StandardZone(real lat, real lon, int setzone) {
43  if (!(setzone >= MINPSEUDOZONE && setzone <= MAXZONE))
44  throw GeographicErr("Illegal zone requested " + Utility::str(setzone));
45  if (setzone >= MINZONE || setzone == INVALID)
46  return setzone;
47  if (Math::isnan(lat) || Math::isnan(lon)) // Check if lat or lon is a NaN
48  return INVALID;
49  if (setzone == UTM || (lat >= -80 && lat < 84)) {
50  // Assume lon is in [-540, 540).
51  int ilon = int(floor(lon));
52  if (ilon >= 180)
53  ilon -= 360;
54  else if (ilon < -180)
55  ilon += 360;
56  int zone = (ilon + 186)/6;
57  int band = MGRS::LatitudeBand(lat);
58  if (band == 7 && zone == 31 && ilon >= 3)
59  zone = 32;
60  else if (band == 9 && ilon >= 0 && ilon < 42)
61  zone = 2 * ((ilon + 183)/12) + 1;
62  return zone;
63  } else
64  return UPS;
65  }
66 
67  void UTMUPS::Forward(real lat, real lon,
68  int& zone, bool& northp, real& x, real& y,
69  real& gamma, real& k,
70  int setzone, bool mgrslimits) {
71  CheckLatLon(lat, lon);
72  bool northp1 = lat >= 0;
73  int zone1 = StandardZone(lat, lon, setzone);
74  if (zone1 == INVALID) {
75  zone = zone1;
76  northp = northp1;
77  x = y = gamma = k = Math::NaN<real>();
78  return;
79  }
80  real x1, y1, gamma1, k1;
81  bool utmp = zone1 != UPS;
82  if (utmp) {
83  real
84  lon0 = CentralMeridian(zone1),
85  dlon = lon - lon0;
86  dlon = abs(dlon - 360 * floor((dlon + 180)/360));
87  if (dlon > 60)
88  // Check isn't really necessary because CheckCoords catches this case.
89  // But this allows a more meaningful error message to be given.
90  throw GeographicErr("Longitude " + Utility::str(lon)
91  + "d more than 60d from center of UTM zone "
92  + Utility::str(zone1));
93  TransverseMercator::UTM.Forward(lon0, lat, lon, x1, y1, gamma1, k1);
94  } else {
95  if (abs(lat) < 70)
96  // Check isn't really necessary ... (see above).
97  throw GeographicErr("Latitude " + Utility::str(lat)
98  + "d more than 20d from "
99  + (northp1 ? "N" : "S") + " pole");
100  PolarStereographic::UPS.Forward(northp1, lat, lon, x1, y1, gamma1, k1);
101  }
102  int ind = (utmp ? 2 : 0) + (northp1 ? 1 : 0);
103  x1 += falseeasting_[ind];
104  y1 += falsenorthing_[ind];
105  if (! CheckCoords(zone1 != UPS, northp1, x1, y1, mgrslimits, false) )
106  throw GeographicErr("Latitude " + Utility::str(lat)
107  + ", longitude " + Utility::str(lon)
108  + " out of legal range for "
109  + (utmp ? "UTM zone " + Utility::str(zone1) : "UPS"));
110  zone = zone1;
111  northp = northp1;
112  x = x1;
113  y = y1;
114  gamma = gamma1;
115  k = k1;
116  }
117 
118  void UTMUPS::Reverse(int zone, bool northp, real x, real y,
119  real& lat, real& lon, real& gamma, real& k,
120  bool mgrslimits) {
121  if (zone == INVALID || Math::isnan(x) || Math::isnan(y)) {
122  lat = lon = gamma = k = Math::NaN<real>();
123  return;
124  }
125  if (!(zone >= MINZONE && zone <= MAXZONE))
126  throw GeographicErr("Zone " + Utility::str(zone)
127  + " not in range [0, 60]");
128  bool utmp = zone != UPS;
129  CheckCoords(utmp, northp, x, y, mgrslimits);
130  int ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
131  x -= falseeasting_[ind];
132  y -= falsenorthing_[ind];
133  if (utmp)
134  TransverseMercator::UTM.Reverse(CentralMeridian(zone),
135  x, y, lat, lon, gamma, k);
136  else
137  PolarStereographic::UPS.Reverse(northp, x, y, lat, lon, gamma, k);
138  }
139 
140  void UTMUPS::CheckLatLon(real lat, real lon) {
141  if (abs(lat) > 90)
142  throw GeographicErr("Latitude " + Utility::str(lat)
143  + "d not in [-90d, 90d]");
144  if (lon < -540 || lon >= 540)
145  throw GeographicErr("Longitude " + Utility::str(lon)
146  + "d not in [-540d, 540d)");
147  }
148 
149  bool UTMUPS::CheckCoords(bool utmp, bool northp, real x, real y,
150  bool mgrslimits, bool throwp) {
151  // Limits are all multiples of 100km and are all closed on the both ends.
152  // Failure tests are such that NaNs succeed.
153  real slop = mgrslimits ? 0 : MGRS::tile_;
154  int ind = (utmp ? 2 : 0) + (northp ? 1 : 0);
155  if (x < mineasting_[ind] - slop || x > maxeasting_[ind] + slop) {
156  if (!throwp) return false;
157  throw GeographicErr("Easting " + Utility::str(x/1000) + "km not in "
158  + (mgrslimits ? "MGRS/" : "")
159  + (utmp ? "UTM" : "UPS") + " range for "
160  + (northp ? "N" : "S" ) + " hemisphere ["
161  + Utility::str((mineasting_[ind] - slop)/1000)
162  + "km, "
163  + Utility::str((maxeasting_[ind] + slop)/1000)
164  + "km]");
165  }
166  if (y < minnorthing_[ind] - slop || y > maxnorthing_[ind] + slop) {
167  if (!throwp) return false;
168  throw GeographicErr("Northing " + Utility::str(y/1000) + "km not in "
169  + (mgrslimits ? "MGRS/" : "")
170  + (utmp ? "UTM" : "UPS") + " range for "
171  + (northp ? "N" : "S" ) + " hemisphere ["
172  + Utility::str((minnorthing_[ind] - slop)/1000)
173  + "km, "
174  + Utility::str((maxnorthing_[ind] + slop)/1000)
175  + "km]");
176  }
177  return true;
178  }
179 
180  void UTMUPS::Transfer(int zonein, bool northpin, real xin, real yin,
181  int zoneout, bool northpout, real& xout, real& yout,
182  int& zone) {
183  bool northp = northpin;
184  if (zonein != zoneout) {
185  // Determine lat, lon
186  real lat, lon;
187  GeographicLib::UTMUPS::Reverse(zonein, northpin, xin, yin, lat, lon);
188  // Try converting to zoneout
189  real x, y;
190  int zone1;
191  GeographicLib::UTMUPS::Forward(lat, lon, zone1, northp, x, y,
192  zoneout == UTMUPS::MATCH
193  ? zonein : zoneout);
194  if (zone1 == 0 && northp != northpout)
195  throw GeographicErr
196  ("Attempt to transfer UPS coordinates between hemispheres");
197  zone = zone1;
198  xout = x;
199  yout = y;
200  } else {
201  if (zoneout == 0 && northp != northpout)
202  throw GeographicErr
203  ("Attempt to transfer UPS coordinates between hemispheres");
204  zone = zoneout;
205  xout = xin;
206  yout = yin;
207  }
208  if (northp != northpout)
209  // Can't get here if UPS
210  yout += (northpout ? -1 : 1) * MGRS::utmNshift_;
211  return;
212  }
213 
214  void UTMUPS::DecodeZone(const std::string& zonestr, int& zone, bool& northp) {
215  unsigned zlen = unsigned(zonestr.size());
216  if (zlen == 0)
217  throw GeographicErr("Empty zone specification");
218  if (zlen > 3)
219  throw GeographicErr("More than 3 characters in zone specification "
220  + zonestr);
221  if (zlen == 3 &&
222  toupper(zonestr[0]) == 'I' &&
223  toupper(zonestr[1]) == 'N' &&
224  toupper(zonestr[2]) == 'V') {
225  zone = INVALID;
226  northp = false;
227  return;
228  }
229  char hemi = char(toupper(zonestr[zlen - 1]));
230  bool northp1 = hemi == 'N';
231  if (! (northp1 || hemi == 'S'))
232  throw GeographicErr(string("Illegal hemisphere letter ") + hemi + " in "
233  + zonestr + ", specify N or S");
234  if (zlen == 1)
235  zone = UPS;
236  else {
237  const char* c = zonestr.c_str();
238  char* q;
239  int zone1 = strtol(c, &q, 10);
240  if (q == c)
241  throw GeographicErr("No zone number found in " + zonestr);
242  if (q - c != int(zlen) - 1)
243  throw GeographicErr("Extra text " +
244  zonestr.substr(q - c, int(zlen) - 1 - (q - c)) +
245  " in UTM/UPS zone " + zonestr);
246  if (zone1 == UPS)
247  // Don't allow 0N as an alternative to N for UPS coordinates
248  throw GeographicErr("Illegal zone 0 in " + zonestr +
249  ", use just " + hemi + " for UPS");
250  if (!(zone1 >= MINUTMZONE && zone1 <= MAXUTMZONE))
251  throw GeographicErr("Zone " + Utility::str(zone1)
252  + " not in range [1, 60]");
253  zone = zone1;
254  }
255  northp = northp1;
256  }
257 
258  std::string UTMUPS::EncodeZone(int zone, bool northp) {
259  if (zone == INVALID)
260  return string("INV");
261  if (!(zone >= MINZONE && zone <= MAXZONE))
262  throw GeographicErr("Zone " + Utility::str(zone)
263  + " not in range [0, 60]");
264  ostringstream os;
265  if (zone != UPS)
266  os << setfill('0') << setw(2) << zone;
267  os << (northp ? 'N' : 'S');
268  return os.str();
269  }
270 
271  void UTMUPS::DecodeEPSG(int& epsg, int& zone, bool& northp) throw() {
272  if (epsg >= epsg01N && epsg <= epsg60N) {
273  zone = epsg - epsg01N + 1;
274  northp = true;
275  } else if (epsg == epsgN) {
276  zone = UPS;
277  northp = true;
278  } else if (epsg >= epsg01S && epsg <= epsg60S) {
279  zone = epsg - epsg01S + 1;
280  northp = false;
281  } else if (epsg == epsgS) {
282  zone = UPS;
283  northp = false;
284  } else {
285  zone = INVALID;
286  northp = false;
287  }
288  }
289 
290  int UTMUPS::EncodeEPSG(int zone, bool northp) throw() {
291  int epsg = -1;
292  if (zone == UPS)
293  epsg = epsgS;
294  else if (zone >= MINUTMZONE && zone <= MAXUTMZONE)
295  epsg = epsg + (zone - MINUTMZONE) + epsg01S;
296  if (epsg >= 0 && northp)
297  epsg += epsgN - epsgS;
298  return epsg;
299  }
300 
301  Math::real UTMUPS::UTMShift() throw() { return real(MGRS::utmNshift_); }
302 
303 } // namespace GeographicLib
Generated by   doxygen 1.8.2