| 1 |
package dk.thoerup.datumconversion; |
| 2 |
|
| 3 |
public class DatumConverter { |
| 4 |
|
| 5 |
static Utm LLtoUTM(int ReferenceEllipsoid, Ll ll) |
| 6 |
{ |
| 7 |
//converts lat/long to UTM coords. Equations from USGS Bulletin 1532 |
| 8 |
//East Longitudes are positive, West longitudes are negative. |
| 9 |
//North latitudes are positive, South latitudes are negative |
| 10 |
//Lat and Long are in decimal degrees |
| 11 |
//Written by Chuck Gantz- chuck.gantz@globalstar.com |
| 12 |
|
| 13 |
|
| 14 |
double UTMNorthing; |
| 15 |
double UTMEasting; |
| 16 |
String UTMZone; |
| 17 |
|
| 18 |
double Lat = ll.lattitude; |
| 19 |
double Long = ll.longitude; |
| 20 |
|
| 21 |
|
| 22 |
|
| 23 |
double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius; |
| 24 |
double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared; |
| 25 |
double k0 = 0.9996; |
| 26 |
|
| 27 |
double LongOrigin; |
| 28 |
double eccPrimeSquared; |
| 29 |
double N, T, C, A, M; |
| 30 |
|
| 31 |
//Make sure the longitude is between -180.00 .. 179.9 |
| 32 |
double LongTemp = (Long+180) - ((int)(Long+180)/360)*360-180; // -180.00 .. 179.9; |
| 33 |
|
| 34 |
double LatRad = Lat * Constants.deg2rad; |
| 35 |
double LongRad = LongTemp * Constants.deg2rad; |
| 36 |
double LongOriginRad; |
| 37 |
int ZoneNumber; |
| 38 |
|
| 39 |
ZoneNumber = ( (int)((LongTemp + 180)/6) ) + 1; |
| 40 |
|
| 41 |
if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 ) |
| 42 |
ZoneNumber = 32; |
| 43 |
|
| 44 |
// Special zones for Svalbard |
| 45 |
if( Lat >= 72.0 && Lat < 84.0 ) |
| 46 |
{ |
| 47 |
if( LongTemp >= 0.0 && LongTemp < 9.0 ) ZoneNumber = 31; |
| 48 |
else if( LongTemp >= 9.0 && LongTemp < 21.0 ) ZoneNumber = 33; |
| 49 |
else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber = 35; |
| 50 |
else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber = 37; |
| 51 |
} |
| 52 |
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone |
| 53 |
LongOriginRad = LongOrigin * Constants.deg2rad; |
| 54 |
|
| 55 |
//compute the UTM Zone from the latitude and longitude |
| 56 |
UTMZone = ""+ ZoneNumber + UTMLetterDesignator(Lat); |
| 57 |
|
| 58 |
eccPrimeSquared = (eccSquared)/(1-eccSquared); |
| 59 |
|
| 60 |
N = a/Math.sqrt(1-eccSquared * Math.sin(LatRad) * Math.sin(LatRad)); |
| 61 |
T = Math.tan(LatRad) * Math.tan(LatRad); |
| 62 |
C = eccPrimeSquared * Math.cos(LatRad) * Math.cos(LatRad); |
| 63 |
A = Math.cos(LatRad)*(LongRad-LongOriginRad); |
| 64 |
|
| 65 |
M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad |
| 66 |
- (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(2*LatRad) |
| 67 |
+ (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(4*LatRad) |
| 68 |
- (35*eccSquared*eccSquared*eccSquared/3072)*Math.sin(6*LatRad)); |
| 69 |
|
| 70 |
UTMEasting = (double)(k0*N*(A+(1-T+C)*A*A*A/6 |
| 71 |
+ (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120) |
| 72 |
+ 500000.0); |
| 73 |
|
| 74 |
UTMNorthing = (double)(k0*(M+N*Math.tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 |
| 75 |
+ (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720))); |
| 76 |
if(Lat < 0) |
| 77 |
UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere |
| 78 |
|
| 79 |
return new Utm(UTMNorthing, UTMEasting,UTMZone); |
| 80 |
} |
| 81 |
|
| 82 |
static char UTMLetterDesignator(double Lat) |
| 83 |
{ |
| 84 |
//This routine determines the correct UTM letter designator for the given latitude |
| 85 |
//returns 'Z' if latitude is outside the UTM limits of 84N to 80S |
| 86 |
//Written by Chuck Gantz- chuck.gantz@globalstar.com |
| 87 |
char LetterDesignator; |
| 88 |
|
| 89 |
if((84 >= Lat) && (Lat >= 72)) LetterDesignator = 'X'; |
| 90 |
else if((72 > Lat) && (Lat >= 64)) LetterDesignator = 'W'; |
| 91 |
else if((64 > Lat) && (Lat >= 56)) LetterDesignator = 'V'; |
| 92 |
else if((56 > Lat) && (Lat >= 48)) LetterDesignator = 'U'; |
| 93 |
else if((48 > Lat) && (Lat >= 40)) LetterDesignator = 'T'; |
| 94 |
else if((40 > Lat) && (Lat >= 32)) LetterDesignator = 'S'; |
| 95 |
else if((32 > Lat) && (Lat >= 24)) LetterDesignator = 'R'; |
| 96 |
else if((24 > Lat) && (Lat >= 16)) LetterDesignator = 'Q'; |
| 97 |
else if((16 > Lat) && (Lat >= 8)) LetterDesignator = 'P'; |
| 98 |
else if(( 8 > Lat) && (Lat >= 0)) LetterDesignator = 'N'; |
| 99 |
else if(( 0 > Lat) && (Lat >= -8)) LetterDesignator = 'M'; |
| 100 |
else if((-8> Lat) && (Lat >= -16)) LetterDesignator = 'L'; |
| 101 |
else if((-16 > Lat) && (Lat >= -24)) LetterDesignator = 'K'; |
| 102 |
else if((-24 > Lat) && (Lat >= -32)) LetterDesignator = 'J'; |
| 103 |
else if((-32 > Lat) && (Lat >= -40)) LetterDesignator = 'H'; |
| 104 |
else if((-40 > Lat) && (Lat >= -48)) LetterDesignator = 'G'; |
| 105 |
else if((-48 > Lat) && (Lat >= -56)) LetterDesignator = 'F'; |
| 106 |
else if((-56 > Lat) && (Lat >= -64)) LetterDesignator = 'E'; |
| 107 |
else if((-64 > Lat) && (Lat >= -72)) LetterDesignator = 'D'; |
| 108 |
else if((-72 > Lat) && (Lat >= -80)) LetterDesignator = 'C'; |
| 109 |
else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits |
| 110 |
|
| 111 |
return LetterDesignator; |
| 112 |
} |
| 113 |
|
| 114 |
|
| 115 |
static Ll UTMtoLL(int ReferenceEllipsoid, Utm utm ) |
| 116 |
{ |
| 117 |
//converts UTM coords to lat/long. Equations from USGS Bulletin 1532 |
| 118 |
//East Longitudes are positive, West longitudes are negative. |
| 119 |
//North latitudes are positive, South latitudes are negative |
| 120 |
//Lat and Long are in decimal degrees. |
| 121 |
//Written by Chuck Gantz- chuck.gantz@globalstar.com |
| 122 |
|
| 123 |
double k0 = 0.9996; |
| 124 |
double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius; |
| 125 |
double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared; |
| 126 |
double eccPrimeSquared; |
| 127 |
double e1 = (1-Math.sqrt(1-eccSquared))/(1+Math.sqrt(1-eccSquared)); |
| 128 |
double N1, T1, C1, R1, D, M; |
| 129 |
double LongOrigin; |
| 130 |
double mu, phi1, phi1Rad; |
| 131 |
double x, y; |
| 132 |
int ZoneNumber; |
| 133 |
char ZoneLetter; |
| 134 |
int NorthernHemisphere; //1 for northern hemispher, 0 for southern |
| 135 |
|
| 136 |
double Lat; |
| 137 |
double Long; |
| 138 |
|
| 139 |
double UTMNorthing = utm.northing; |
| 140 |
double UTMEasting = utm.easting; |
| 141 |
String UTMZone = utm.utmZone; |
| 142 |
|
| 143 |
x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude |
| 144 |
y = UTMNorthing; |
| 145 |
|
| 146 |
ZoneNumber = Integer.parseInt(UTMZone.substring(0,UTMZone.length()-1)); |
| 147 |
ZoneLetter = UTMZone.charAt(UTMZone.length()-1); |
| 148 |
//ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10); |
| 149 |
if((ZoneLetter - 'N') >= 0) |
| 150 |
NorthernHemisphere = 1;//point is in northern hemisphere |
| 151 |
else |
| 152 |
{ |
| 153 |
NorthernHemisphere = 0;//point is in southern hemisphere |
| 154 |
y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere |
| 155 |
} |
| 156 |
|
| 157 |
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone |
| 158 |
|
| 159 |
eccPrimeSquared = (eccSquared)/(1-eccSquared); |
| 160 |
|
| 161 |
M = y / k0; |
| 162 |
mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256)); |
| 163 |
|
| 164 |
phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*Math.sin(2*mu) |
| 165 |
+ (21*e1*e1/16-55*e1*e1*e1*e1/32)*Math.sin(4*mu) |
| 166 |
+(151*e1*e1*e1/96)*Math.sin(6*mu); |
| 167 |
phi1 = phi1Rad*Constants.rad2deg; |
| 168 |
|
| 169 |
N1 = a/Math.sqrt(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad)); |
| 170 |
T1 = Math.tan(phi1Rad)*Math.tan(phi1Rad); |
| 171 |
C1 = eccPrimeSquared*Math.cos(phi1Rad)*Math.cos(phi1Rad); |
| 172 |
R1 = a*(1-eccSquared)/Math.pow(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad), 1.5); |
| 173 |
D = x/(N1*k0); |
| 174 |
|
| 175 |
Lat = phi1Rad - (N1*Math.tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24 |
| 176 |
+(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720); |
| 177 |
Lat = Lat * Constants.rad2deg; |
| 178 |
|
| 179 |
Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1) |
| 180 |
*D*D*D*D*D/120)/Math.cos(phi1Rad); |
| 181 |
Long = LongOrigin + Long * Constants.rad2deg; |
| 182 |
|
| 183 |
return new Ll(Lat,Long); |
| 184 |
} |
| 185 |
|
| 186 |
/////////////////////////////////////////////////////////////////////////////////////////////////////// |
| 187 |
|
| 188 |
static Swiss LLtoSwissGrid(double Lat, double Long) |
| 189 |
{ |
| 190 |
//converts lat/long to Swiss Grid coords. Equations from "Supplementary PROJ.4 Notes- |
| 191 |
//Swiss Oblique Mercator Projection", August 5, 1995, Release 4.3.3, by Gerald I. Evenden |
| 192 |
//Lat and Long are in decimal degrees |
| 193 |
//This transformation is, of course, only valid in Switzerland |
| 194 |
//Written by Chuck Gantz- chuck.gantz@globalstar.com |
| 195 |
double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid |
| 196 |
double eccSquared = Constants.ellipsoid[3].eccentricitySquared; |
| 197 |
double ecc = Math.sqrt(eccSquared); |
| 198 |
|
| 199 |
double LongOrigin = 7.43958333; //E7d26'22.500" |
| 200 |
double LatOrigin = 46.95240556; //N46d57'8.660" |
| 201 |
|
| 202 |
double SwissNorthing; |
| 203 |
double SwissEasting; |
| 204 |
|
| 205 |
double LatRad = Lat*Constants.deg2rad; |
| 206 |
double LongRad = Long*Constants.deg2rad; |
| 207 |
double LatOriginRad = LatOrigin*Constants.deg2rad; |
| 208 |
double LongOriginRad = LongOrigin*Constants.deg2rad; |
| 209 |
|
| 210 |
double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); |
| 211 |
|
| 212 |
double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c); |
| 213 |
|
| 214 |
//eqn. 1 |
| 215 |
double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2)) |
| 216 |
-c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2)) |
| 217 |
- ecc/2 * Math.log((1+ecc*Math.sin(LatOriginRad)) / (1-ecc*Math.sin(LatOriginRad)))); |
| 218 |
|
| 219 |
|
| 220 |
double LongRadPrime = c*(LongRad - LongOriginRad); //eqn 2 |
| 221 |
double w = c*(Math.log(Math.tan(Constants.FOURTHPI + LatRad/2)) |
| 222 |
- ecc/2 * Math.log((1+ecc*Math.sin(LatRad)) / (1-ecc*Math.sin(LatRad)))) + K; //eqn 1 |
| 223 |
double LatRadPrime = 2 * (Math.atan(Math.exp(w)) - Constants.FOURTHPI); //eqn 1 |
| 224 |
|
| 225 |
//eqn 3 |
| 226 |
double sinLatDoublePrime = Math.cos(equivLatOrgRadPrime) * Math.sin(LatRadPrime) |
| 227 |
- Math.sin(equivLatOrgRadPrime) * Math.cos(LatRadPrime) * Math.cos(LongRadPrime); |
| 228 |
double LatRadDoublePrime = Math.asin(sinLatDoublePrime); |
| 229 |
|
| 230 |
//eqn 4 |
| 231 |
double sinLongDoublePrime = Math.cos(LatRadPrime)*Math.sin(LongRadPrime) / Math.cos(LatRadDoublePrime); |
| 232 |
double LongRadDoublePrime = Math.asin(sinLongDoublePrime); |
| 233 |
|
| 234 |
double R = a*Math.sqrt(1-eccSquared) / (1-eccSquared*Math.sin(LatOriginRad) * Math.sin(LatOriginRad)); |
| 235 |
|
| 236 |
SwissNorthing = R*Math.log(Math.tan(Constants.FOURTHPI + LatRadDoublePrime/2)) + 200000.0; //eqn 5 |
| 237 |
SwissEasting = R*LongRadDoublePrime + 600000.0; //eqn 6 |
| 238 |
|
| 239 |
return new Swiss(SwissNorthing, SwissEasting); |
| 240 |
} |
| 241 |
|
| 242 |
|
| 243 |
static Ll SwissGridtoLL(double SwissNorthing, double SwissEasting) |
| 244 |
{ |
| 245 |
double Lat; |
| 246 |
double Long; |
| 247 |
|
| 248 |
double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid |
| 249 |
double eccSquared = Constants.ellipsoid[3].eccentricitySquared; |
| 250 |
double ecc = Math.sqrt(eccSquared); |
| 251 |
|
| 252 |
double LongOrigin = 7.43958333; //E7d26'22.500" |
| 253 |
double LatOrigin = 46.95240556; //N46d57'8.660" |
| 254 |
|
| 255 |
double LatOriginRad = LatOrigin*Constants.deg2rad; |
| 256 |
double LongOriginRad = LongOrigin*Constants.deg2rad; |
| 257 |
|
| 258 |
double R = a*Math.sqrt(1-eccSquared) / (1-eccSquared*Math.sin(LatOriginRad) * Math.sin(LatOriginRad)); |
| 259 |
|
| 260 |
double LatRadDoublePrime = 2*(Math.atan(Math.exp((SwissNorthing - 200000.0)/R)) - Constants.FOURTHPI); //eqn. 7 |
| 261 |
double LongRadDoublePrime = (SwissEasting - 600000.0)/R; //eqn. 8 with equation corrected |
| 262 |
|
| 263 |
|
| 264 |
double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); |
| 265 |
double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c); |
| 266 |
|
| 267 |
double sinLatRadPrime = Math.cos(equivLatOrgRadPrime)*Math.sin(LatRadDoublePrime) |
| 268 |
+ Math.sin(equivLatOrgRadPrime)*Math.cos(LatRadDoublePrime)*Math.cos(LongRadDoublePrime); |
| 269 |
double LatRadPrime = Math.asin(sinLatRadPrime); |
| 270 |
|
| 271 |
double sinLongRadPrime = Math.cos(LatRadDoublePrime)*Math.sin(LongRadDoublePrime)/Math.cos(LatRadPrime); |
| 272 |
double LongRadPrime = Math.asin(sinLongRadPrime); |
| 273 |
|
| 274 |
Long = (LongRadPrime/c + LongOriginRad) * Constants.rad2deg; |
| 275 |
|
| 276 |
Lat = NewtonRaphson(LatRadPrime) * Constants.rad2deg; |
| 277 |
|
| 278 |
return new Ll(Lat, Long); |
| 279 |
} |
| 280 |
|
| 281 |
static double NewtonRaphson( double initEstimate) |
| 282 |
{ |
| 283 |
double Estimate = initEstimate; |
| 284 |
double tol = 0.00001; |
| 285 |
double corr; |
| 286 |
|
| 287 |
double eccSquared = Constants.ellipsoid[3].eccentricitySquared; |
| 288 |
double ecc = Math.sqrt(eccSquared); |
| 289 |
|
| 290 |
double LatOrigin = 46.95240556; //N46d57'8.660" |
| 291 |
double LatOriginRad = LatOrigin*Constants.deg2rad; |
| 292 |
|
| 293 |
double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); |
| 294 |
|
| 295 |
double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c); |
| 296 |
|
| 297 |
//eqn. 1 |
| 298 |
double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2)) |
| 299 |
-c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2)) |
| 300 |
- ecc/2 * Math.log((1+ecc*Math.sin(LatOriginRad)) / (1-ecc*Math.sin(LatOriginRad)))); |
| 301 |
double C = (K - Math.log(Math.tan(Constants.FOURTHPI + initEstimate/2)))/c; |
| 302 |
|
| 303 |
do |
| 304 |
{ |
| 305 |
corr = CorrRatio(Estimate, C); |
| 306 |
Estimate = Estimate - corr; |
| 307 |
} |
| 308 |
while (Math.abs(corr) > tol); |
| 309 |
|
| 310 |
return Estimate; |
| 311 |
} |
| 312 |
|
| 313 |
|
| 314 |
|
| 315 |
static double CorrRatio(double LatRad, double C) |
| 316 |
{ |
| 317 |
double eccSquared = Constants.ellipsoid[3].eccentricitySquared; |
| 318 |
double ecc = Math.sqrt(eccSquared); |
| 319 |
double corr = (C + Math.log(Math.tan(Constants.FOURTHPI + LatRad/2)) |
| 320 |
- ecc/2 * Math.log((1+ecc*Math.sin(LatRad)) / (1-ecc*Math.sin(LatRad)))) * (((1-eccSquared*Math.sin(LatRad)*Math.sin(LatRad)) * Math.cos(LatRad)) / (1-eccSquared)); |
| 321 |
|
| 322 |
return corr; |
| 323 |
} |
| 324 |
|
| 325 |
|
| 326 |
} |
Parent Directory
| 
Revision Log