package dk.thoerup.datumconversion; public class DatumConverter { static Utm LLtoUTM(int ReferenceEllipsoid, Ll ll) { //converts lat/long to UTM coords. Equations from USGS Bulletin 1532 //East Longitudes are positive, West longitudes are negative. //North latitudes are positive, South latitudes are negative //Lat and Long are in decimal degrees //Written by Chuck Gantz- chuck.gantz@globalstar.com double UTMNorthing; double UTMEasting; String UTMZone; double Lat = ll.lattitude; double Long = ll.longitude; double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius; double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared; double k0 = 0.9996; double LongOrigin; double eccPrimeSquared; double N, T, C, A, M; //Make sure the longitude is between -180.00 .. 179.9 double LongTemp = (Long+180) - ((int)(Long+180)/360)*360-180; // -180.00 .. 179.9; double LatRad = Lat * Constants.deg2rad; double LongRad = LongTemp * Constants.deg2rad; double LongOriginRad; int ZoneNumber; ZoneNumber = ( (int)((LongTemp + 180)/6) ) + 1; if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 ) ZoneNumber = 32; // Special zones for Svalbard if( Lat >= 72.0 && Lat < 84.0 ) { if( LongTemp >= 0.0 && LongTemp < 9.0 ) ZoneNumber = 31; else if( LongTemp >= 9.0 && LongTemp < 21.0 ) ZoneNumber = 33; else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber = 35; else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber = 37; } LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone LongOriginRad = LongOrigin * Constants.deg2rad; //compute the UTM Zone from the latitude and longitude UTMZone = ""+ ZoneNumber + UTMLetterDesignator(Lat); eccPrimeSquared = (eccSquared)/(1-eccSquared); N = a/Math.sqrt(1-eccSquared * Math.sin(LatRad) * Math.sin(LatRad)); T = Math.tan(LatRad) * Math.tan(LatRad); C = eccPrimeSquared * Math.cos(LatRad) * Math.cos(LatRad); A = Math.cos(LatRad)*(LongRad-LongOriginRad); M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(2*LatRad) + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*Math.sin(4*LatRad) - (35*eccSquared*eccSquared*eccSquared/3072)*Math.sin(6*LatRad)); UTMEasting = (double)(k0*N*(A+(1-T+C)*A*A*A/6 + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120) + 500000.0); UTMNorthing = (double)(k0*(M+N*Math.tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720))); if(Lat < 0) UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere return new Utm(UTMNorthing, UTMEasting,UTMZone); } static char UTMLetterDesignator(double Lat) { //This routine determines the correct UTM letter designator for the given latitude //returns 'Z' if latitude is outside the UTM limits of 84N to 80S //Written by Chuck Gantz- chuck.gantz@globalstar.com char LetterDesignator; if((84 >= Lat) && (Lat >= 72)) LetterDesignator = 'X'; else if((72 > Lat) && (Lat >= 64)) LetterDesignator = 'W'; else if((64 > Lat) && (Lat >= 56)) LetterDesignator = 'V'; else if((56 > Lat) && (Lat >= 48)) LetterDesignator = 'U'; else if((48 > Lat) && (Lat >= 40)) LetterDesignator = 'T'; else if((40 > Lat) && (Lat >= 32)) LetterDesignator = 'S'; else if((32 > Lat) && (Lat >= 24)) LetterDesignator = 'R'; else if((24 > Lat) && (Lat >= 16)) LetterDesignator = 'Q'; else if((16 > Lat) && (Lat >= 8)) LetterDesignator = 'P'; else if(( 8 > Lat) && (Lat >= 0)) LetterDesignator = 'N'; else if(( 0 > Lat) && (Lat >= -8)) LetterDesignator = 'M'; else if((-8> Lat) && (Lat >= -16)) LetterDesignator = 'L'; else if((-16 > Lat) && (Lat >= -24)) LetterDesignator = 'K'; else if((-24 > Lat) && (Lat >= -32)) LetterDesignator = 'J'; else if((-32 > Lat) && (Lat >= -40)) LetterDesignator = 'H'; else if((-40 > Lat) && (Lat >= -48)) LetterDesignator = 'G'; else if((-48 > Lat) && (Lat >= -56)) LetterDesignator = 'F'; else if((-56 > Lat) && (Lat >= -64)) LetterDesignator = 'E'; else if((-64 > Lat) && (Lat >= -72)) LetterDesignator = 'D'; else if((-72 > Lat) && (Lat >= -80)) LetterDesignator = 'C'; else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits return LetterDesignator; } static Ll UTMtoLL(int ReferenceEllipsoid, Utm utm ) { //converts UTM coords to lat/long. Equations from USGS Bulletin 1532 //East Longitudes are positive, West longitudes are negative. //North latitudes are positive, South latitudes are negative //Lat and Long are in decimal degrees. //Written by Chuck Gantz- chuck.gantz@globalstar.com double k0 = 0.9996; double a = Constants.ellipsoid[ReferenceEllipsoid].EquatorialRadius; double eccSquared = Constants.ellipsoid[ReferenceEllipsoid].eccentricitySquared; double eccPrimeSquared; double e1 = (1-Math.sqrt(1-eccSquared))/(1+Math.sqrt(1-eccSquared)); double N1, T1, C1, R1, D, M; double LongOrigin; double mu, phi1, phi1Rad; double x, y; int ZoneNumber; char ZoneLetter; int NorthernHemisphere; //1 for northern hemispher, 0 for southern double Lat; double Long; double UTMNorthing = utm.northing; double UTMEasting = utm.easting; String UTMZone = utm.utmZone; x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude y = UTMNorthing; ZoneNumber = Integer.parseInt(UTMZone.substring(0,UTMZone.length()-1)); ZoneLetter = UTMZone.charAt(UTMZone.length()-1); //ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10); if((ZoneLetter - 'N') >= 0) NorthernHemisphere = 1;//point is in northern hemisphere else { NorthernHemisphere = 0;//point is in southern hemisphere y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere } LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone eccPrimeSquared = (eccSquared)/(1-eccSquared); M = y / k0; mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256)); phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*Math.sin(2*mu) + (21*e1*e1/16-55*e1*e1*e1*e1/32)*Math.sin(4*mu) +(151*e1*e1*e1/96)*Math.sin(6*mu); phi1 = phi1Rad*Constants.rad2deg; N1 = a/Math.sqrt(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad)); T1 = Math.tan(phi1Rad)*Math.tan(phi1Rad); C1 = eccPrimeSquared*Math.cos(phi1Rad)*Math.cos(phi1Rad); R1 = a*(1-eccSquared)/Math.pow(1-eccSquared*Math.sin(phi1Rad)*Math.sin(phi1Rad), 1.5); D = x/(N1*k0); 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 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720); Lat = Lat * Constants.rad2deg; Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1) *D*D*D*D*D/120)/Math.cos(phi1Rad); Long = LongOrigin + Long * Constants.rad2deg; return new Ll(Lat,Long); } /////////////////////////////////////////////////////////////////////////////////////////////////////// static Swiss LLtoSwissGrid(double Lat, double Long) { //converts lat/long to Swiss Grid coords. Equations from "Supplementary PROJ.4 Notes- //Swiss Oblique Mercator Projection", August 5, 1995, Release 4.3.3, by Gerald I. Evenden //Lat and Long are in decimal degrees //This transformation is, of course, only valid in Switzerland //Written by Chuck Gantz- chuck.gantz@globalstar.com double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid double eccSquared = Constants.ellipsoid[3].eccentricitySquared; double ecc = Math.sqrt(eccSquared); double LongOrigin = 7.43958333; //E7d26'22.500" double LatOrigin = 46.95240556; //N46d57'8.660" double SwissNorthing; double SwissEasting; double LatRad = Lat*Constants.deg2rad; double LongRad = Long*Constants.deg2rad; double LatOriginRad = LatOrigin*Constants.deg2rad; double LongOriginRad = LongOrigin*Constants.deg2rad; double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c); //eqn. 1 double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2)) -c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2)) - ecc/2 * Math.log((1+ecc*Math.sin(LatOriginRad)) / (1-ecc*Math.sin(LatOriginRad)))); double LongRadPrime = c*(LongRad - LongOriginRad); //eqn 2 double w = c*(Math.log(Math.tan(Constants.FOURTHPI + LatRad/2)) - ecc/2 * Math.log((1+ecc*Math.sin(LatRad)) / (1-ecc*Math.sin(LatRad)))) + K; //eqn 1 double LatRadPrime = 2 * (Math.atan(Math.exp(w)) - Constants.FOURTHPI); //eqn 1 //eqn 3 double sinLatDoublePrime = Math.cos(equivLatOrgRadPrime) * Math.sin(LatRadPrime) - Math.sin(equivLatOrgRadPrime) * Math.cos(LatRadPrime) * Math.cos(LongRadPrime); double LatRadDoublePrime = Math.asin(sinLatDoublePrime); //eqn 4 double sinLongDoublePrime = Math.cos(LatRadPrime)*Math.sin(LongRadPrime) / Math.cos(LatRadDoublePrime); double LongRadDoublePrime = Math.asin(sinLongDoublePrime); double R = a*Math.sqrt(1-eccSquared) / (1-eccSquared*Math.sin(LatOriginRad) * Math.sin(LatOriginRad)); SwissNorthing = R*Math.log(Math.tan(Constants.FOURTHPI + LatRadDoublePrime/2)) + 200000.0; //eqn 5 SwissEasting = R*LongRadDoublePrime + 600000.0; //eqn 6 return new Swiss(SwissNorthing, SwissEasting); } static Ll SwissGridtoLL(double SwissNorthing, double SwissEasting) { double Lat; double Long; double a = Constants.ellipsoid[3].EquatorialRadius; //Bessel ellipsoid double eccSquared = Constants.ellipsoid[3].eccentricitySquared; double ecc = Math.sqrt(eccSquared); double LongOrigin = 7.43958333; //E7d26'22.500" double LatOrigin = 46.95240556; //N46d57'8.660" double LatOriginRad = LatOrigin*Constants.deg2rad; double LongOriginRad = LongOrigin*Constants.deg2rad; double R = a*Math.sqrt(1-eccSquared) / (1-eccSquared*Math.sin(LatOriginRad) * Math.sin(LatOriginRad)); double LatRadDoublePrime = 2*(Math.atan(Math.exp((SwissNorthing - 200000.0)/R)) - Constants.FOURTHPI); //eqn. 7 double LongRadDoublePrime = (SwissEasting - 600000.0)/R; //eqn. 8 with equation corrected double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c); double sinLatRadPrime = Math.cos(equivLatOrgRadPrime)*Math.sin(LatRadDoublePrime) + Math.sin(equivLatOrgRadPrime)*Math.cos(LatRadDoublePrime)*Math.cos(LongRadDoublePrime); double LatRadPrime = Math.asin(sinLatRadPrime); double sinLongRadPrime = Math.cos(LatRadDoublePrime)*Math.sin(LongRadDoublePrime)/Math.cos(LatRadPrime); double LongRadPrime = Math.asin(sinLongRadPrime); Long = (LongRadPrime/c + LongOriginRad) * Constants.rad2deg; Lat = NewtonRaphson(LatRadPrime) * Constants.rad2deg; return new Ll(Lat, Long); } static double NewtonRaphson( double initEstimate) { double Estimate = initEstimate; double tol = 0.00001; double corr; double eccSquared = Constants.ellipsoid[3].eccentricitySquared; double ecc = Math.sqrt(eccSquared); double LatOrigin = 46.95240556; //N46d57'8.660" double LatOriginRad = LatOrigin*Constants.deg2rad; double c = Math.sqrt(1+((eccSquared * Math.pow(Math.cos(LatOriginRad), 4)) / (1-eccSquared))); double equivLatOrgRadPrime = Math.asin(Math.sin(LatOriginRad) / c); //eqn. 1 double K = Math.log(Math.tan(Constants.FOURTHPI + equivLatOrgRadPrime/2)) -c*(Math.log(Math.tan(Constants.FOURTHPI + LatOriginRad/2)) - ecc/2 * Math.log((1+ecc*Math.sin(LatOriginRad)) / (1-ecc*Math.sin(LatOriginRad)))); double C = (K - Math.log(Math.tan(Constants.FOURTHPI + initEstimate/2)))/c; do { corr = CorrRatio(Estimate, C); Estimate = Estimate - corr; } while (Math.abs(corr) > tol); return Estimate; } static double CorrRatio(double LatRad, double C) { double eccSquared = Constants.ellipsoid[3].eccentricitySquared; double ecc = Math.sqrt(eccSquared); double corr = (C + Math.log(Math.tan(Constants.FOURTHPI + LatRad/2)) - 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)); return corr; } }