/* * Copyright (C) 2013 Dr. John Lindsay <jlindsay@uoguelph.ca> * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program 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 General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ package whitebox.georeference; /* * Based on Steven Dutch (2010) http://www.uwgb.edu/dutchs/UsefulData/ConvertUTMNoOZ.HTM; * * The conversion from latitude and longitude to UTM is probably accurate to a * meter or so. The reverse conversion is a bit less accurate, but is * generally accurate within ten meters. * * Check accuracy for your purposes by inputting one set of coordinates, then * converting back. See how close the values are to the original values */ /** * * @author johnlindsay */ public class UTM2LL { /* * Based on Steven Dutch (2010) http://www.uwgb.edu/dutchs/UsefulData/ConvertUTMNoOZ.HTM; */ // global variables private Ellipsoid ellipsoid; private final double drad = Math.PI / 180; private double latitude; private double longitude; private boolean isNorthern = true; private String utmZone; private double utmz; // constructors public UTM2LL() { // no-args constructor } public UTM2LL(Ellipsoid ellipsoid, String utmZone) { this.ellipsoid = ellipsoid; this.utmZone = utmZone; if (utmZone.toLowerCase().contains("n")) { isNorthern = true; utmz = Double.parseDouble(utmZone.toLowerCase().replace("n", "")); } else { isNorthern = false; utmz = Double.parseDouble(utmZone.toLowerCase().replace("s", "")); } } // Properties public Ellipsoid getEllipsoid() { return ellipsoid; } public void setEllipsoid(Ellipsoid ellipsoid) { this.ellipsoid = ellipsoid; } public double getLatitude() { return latitude; } public double getLongitude() { return longitude; } public String getUtmZone() { return utmZone; } public void setUtmZone(String utmZone) { this.utmZone = utmZone; if (utmZone.toLowerCase().contains("n")) { isNorthern = true; utmz = Double.parseDouble(utmZone.toLowerCase().replace("n", "")); } else { isNorthern = false; utmz = Double.parseDouble(utmZone.toLowerCase().replace("s", "")); } } // Methods public void convertUTMCoordinates(double x, double y) { double M, a, b, f, e, e0, esq, e0sq, zcm, e1, M0, mu, phi1, C1, T1, N1, R1, D, phi; double k0 = 0.9996; double k = 1;//local scale a = ellipsoid.majorAxis(); b = ellipsoid.minorAxis(); e = ellipsoid.firstEccentricity(); e0 = e / Math.sqrt(1 - e * e);//Called e prime in reference esq = (1 - (b / a) * (b/a));//e squared for use in expansions e0sq = e * e / (1 - e * e);// e0 squared - always even powers if (x < 160000 || x > 840000){ //Outside permissible range of easting values. Results may be unreliable Use with caution; return; } if (y < 0) { //Negative values not allowed. Results may be unreliable. Use with caution return; } if (y > 10000000) { //Northing may not exceed 10,000,000. Results may be unreliable. Use with caution return; } zcm = 3 + 6*(utmz-1) - 180;//Central meridian of zone e1 = (1 - Math.sqrt(1 - e*e))/(1 + Math.sqrt(1 - e*e));//Called e1 in USGS PP 1395 also M0 = 0;//In case origin other than zero lat - not needed for standard UTM M = M0 + y/k0;//Arc length along standard meridian. if (!isNorthern){ M = M0 + (y - 10000000) / k; } mu = M / (a * (1 - esq * (1 / 4.0 + esq * (3.0 / 64.0 + 5 * esq / 256.0)))); phi1 = mu + e1*(3.0 / 2.0 - 27*e1*e1 / 32.0)*Math.sin(2*mu) + e1*e1*(21/16 -55*e1*e1/32)*Math.sin(4*mu);//Footprint Latitude phi1 = phi1 + e1 * e1 * e1 * (Math.sin(6 * mu) * 151.0 / 96.0 + e1 * Math.sin(8.0 * mu) * 1097.0 / 512.0); C1 = e0sq * Math.pow(Math.cos(phi1), 2); T1 = Math.pow(Math.tan(phi1), 2); N1 = a / Math.sqrt(1 - Math.pow(e * Math.sin(phi1), 2)); R1 = N1 * (1 - e * e) / (1 - Math.pow(e * Math.sin(phi1), 2)); D = (x - 500000) / (N1 * k0); phi = (D * D) * (1.0 / 2.0 - D * D * (5 + 3 * T1 + 10 * C1 - 4.0 * C1 * C1 - 9 * e0sq) / 24.0); phi = phi + Math.pow(D , 6) * (61.0 + 90.0 * T1 + 298.0 * C1 + 45.0 * T1 * T1 - 252.0 * e0sq - 3 * C1 * C1) / 720.0; phi = phi1 - (N1 * Math.tan(phi1) / R1) * phi; latitude = Math.floor(1000000 * phi / drad) / 1000000; double lng = D * (1 + D * D * ((-1 - 2 * T1 - C1) / 6.0 + D * D * (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * C1 * C1 + 8.0 * e0sq + 24.0 * T1 * T1) / 120.0)) / Math.cos(phi1); double lngd = zcm + lng / drad; longitude = Math.floor(1000000 * lngd) / 1000000; } public static void main(String[] args) { UTM2LL utm2ll = new UTM2LL(Ellipsoid.WGS_84, "18N"); double easting = 627103.0885902103; double northing = 4484335.479356929; String utmZone = "18N"; utm2ll.setUtmZone(utmZone); utm2ll.convertUTMCoordinates(easting, northing); System.out.println("Latitude: " + utm2ll.latitude + "\tLongitude: " + utm2ll.longitude); easting = 560560.5471820442; northing = 4822000.781513004; utmZone = "17N"; utm2ll.setUtmZone(utmZone); utm2ll.convertUTMCoordinates(easting, northing); System.out.println("Latitude: " + utm2ll.latitude + "\tLongitude: " + utm2ll.longitude); } }