1 /**************************************************************************\
2 * Copyright (c) Kongsberg Oil & Gas Technologies AS
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are
7 * met:
8 *
9 * Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 *
16 * Neither the name of the copyright holder nor the names of its
17 * contributors may be used to endorse or promote products derived from
18 * this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 \**************************************************************************/
32
33 #include "SbUTMProjection.h"
34 #include <cstdio>
35
SbUTMProjection(const int utmzone,const SbGeoEllipsoid & ellipsoid,double FE,double FN)36 SbUTMProjection::SbUTMProjection(const int utmzone,
37 const SbGeoEllipsoid & ellipsoid,
38 double FE, double FN)
39 : inherited(ellipsoid, FE, FN),
40 forcedutmzone(utmzone)
41 {
42 }
43
44 SbBool
isUTMProjection(void) const45 SbUTMProjection::isUTMProjection(void) const
46 {
47 return TRUE;
48 }
49
50 void
project(const SbGeoAngle & LatRad,const SbGeoAngle & LongRad,double * easting,double * northing) const51 SbUTMProjection::project(const SbGeoAngle & LatRad,
52 const SbGeoAngle & LongRad,
53 double * easting,
54 double * northing) const
55 {
56 double a = this->ellipsoid.getA();
57 double eccSquared = this->ellipsoid.getEccentricitySquared();
58 double k0 = 0.9996;
59
60 double LongOrigin;
61 double eccPrimeSquared;
62 double N, T, C, A, M;
63
64 //Make sure the longitude is between -180.00 .. 179.9
65 double LongTemp = (int(LongRad.deg())+180)-int((int(LongRad.deg())+180)/360)*360-180; // -180.00 .. 179.9;
66
67 // double LatRad = Lat*deg2rad;
68 // double LongRad = LongTemp*deg2rad;
69 double LongOriginRad;
70 int ZoneNumber;
71
72 ZoneNumber = int((LongTemp + 180)/6) + 1;
73
74 if (this->forcedutmzone != -1) {
75 ZoneNumber = this->forcedutmzone;
76 }
77
78 const double deg2rad = M_PI / 180;
79 //const double rad2deg = 180.0 / M_PI;
80 //const double FOURTHPI = M_PI / 4;
81
82 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
83 LongOriginRad = LongOrigin * deg2rad;
84
85 eccPrimeSquared = (eccSquared)/(1-eccSquared);
86
87 N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
88 T = tan(LatRad)*tan(LatRad);
89 C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
90 A = cos(LatRad)*(LongRad-LongOriginRad);
91
92 M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
93 - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2.0*LatRad)
94 + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4.0*LatRad)
95 - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6.0*LatRad));
96
97 *easting = (double)(k0*N*(A+(1-T+C)*A*A*A/6
98 + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
99 + 500000.0);
100
101 *northing = (double)(k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
102 + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
103 if (double(LatRad) < 0.0)
104 *northing += 10000000.0; //10000000 meter offset for southern hemisphere
105 }
106
107 void
unproject(const double UTMEasting,const double UTMNorthing,SbGeoAngle * Lat,SbGeoAngle * Long) const108 SbUTMProjection::unproject(const double UTMEasting,
109 const double UTMNorthing,
110 SbGeoAngle * Lat,
111 SbGeoAngle * Long) const
112 {
113 double k0 = 0.9996;
114 double a = this->ellipsoid.getA();
115 double eccSquared = this->ellipsoid.getEccentricitySquared();
116 double eccPrimeSquared;
117 double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
118 double N1, T1, C1, R1, D, M;
119 double LongOrigin;
120 double mu,/* phi1,*/ phi1Rad;
121 double x, y;
122 int ZoneNumber;
123 //char* ZoneLetter;
124 //int NorthernHemisphere; //1 for northern hemispher, 0 for southern
125
126 x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
127 y = UTMNorthing;
128
129 ZoneNumber = 34; // default zone
130 if (this->forcedutmzone != -1) {
131 ZoneNumber = this->forcedutmzone;
132 }
133 if (this->ellipsoid.getHemisphere() == 'S') {
134 y -= 10000000.0; //remove 10,000,000 meter offset used for southern hemisphere
135 }
136
137 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
138 eccPrimeSquared = (eccSquared)/(1-eccSquared);
139
140 M = y / k0;
141 mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
142
143 phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
144 + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
145 +(151*e1*e1*e1/96)*sin(6*mu);
146
147 N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
148 T1 = tan(phi1Rad)*tan(phi1Rad);
149 C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
150 R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
151 D = x/(N1*k0);
152
153 *Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
154 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
155
156 double tmp = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
157 *D*D*D*D*D/120)/cos(phi1Rad);
158
159 tmp += LongOrigin * M_PI / 180.0;
160 *Long = tmp;
161 }
162
163 void
setUTMZone(const int utmzone)164 SbUTMProjection::setUTMZone(const int utmzone)
165 {
166 this->forcedutmzone = utmzone;
167 }
168
169 int
getUTMZone(void) const170 SbUTMProjection::getUTMZone(void) const
171 {
172 return this->forcedutmzone;
173 }
174