/*
* Copyright (C) 2011-2014 GeoForge Project
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package org.geoforge.lang.util.geography;
import org.geoforge.lang.util.geography.point.GfrPointDms;
import org.geoforge.lang.util.geography.point.GfrPointUtm;
import org.geoforge.lang.util.geography.projection.GfrPrjAbs;
import org.geoforge.lang.util.geography.projection.cylindrical.PrjCylAbs;
/**
*
* @author Amadeus.Sowerby
*
* email: Amadeus.Sowerby_AT_gmail.com
* ... please remove "_AT_" from the above string to get the right email address
*/
public class GfrUtilConversionCylindricalUtm
{
public static String s_getUtmZoneLetter(double phi)
{
if (phi > 90 || phi < -90)
return null;
if (phi > 84)
return "Z";
if (phi >= 72)
return "X";
if (phi >= 64)
return "W";
if (phi >= 56)
return "V";
if (phi >= 48)
return "U";
if (phi >= 40)
return "T";
if (phi >= 32)
return "S";
if (phi >= 24)
return "R";
if (phi >= 16)
return "Q";
if (phi >= 8)
return "P";
if (phi >= 0)
return "N";
if (phi >= -8)
return "M";
if (phi >= -16)
return "L";
if (phi >= -24)
return "K";
if (phi >= -32)
return "J";
if (phi >= -40)
return "H";
if (phi >= -48)
return "G";
if (phi >= -56)
return "F";
if (phi >= -64)
return "E";
if (phi >= -72)
return "D";
if (phi >= -80)
return "C";
return "Z";
}
public static int s_getUtmZoneNumber(double lambda)
{
if (lambda > 0)
{
return ((int) lambda) / ((int) 6) + 31;
}
return ((int) lambda) / ((int) 6) + 30;
}
public static GfrPointUtm s_getUtmCoordinateFromDms(
GfrPointDms pnt,
PrjCylAbs prj)
{
double phiRad = Math.toRadians(pnt.getLatitude());
double lambdaDeg = pnt.getLongitude();
double lambda0 = Math.toRadians(GfrUtilConversionCylindricalUtm._s_getLambda0_(lambdaDeg));
double lambdaRad = Math.toRadians(lambdaDeg);
double a = prj.getEllipsoid().getDblRadius();
double e_2 = prj.getEllipsoid().getDblSquareExcentricity();
double sin_phi = Math.sin(phiRad);
double tan_phi = Math.tan(phiRad);
double sin_2_phi = sin_phi * sin_phi;
double cos_phi = Math.cos(phiRad);
double cos_2_phi = cos_phi * cos_phi;
double nu = 1D / (Math.sqrt(1D - e_2 * sin_2_phi));
double A = (lambdaRad - lambda0) * cos_phi;
double A_2 = A * A;
double A_3 = A_2 * A;
double A_4 = A_3 * A;
double A_5 = A_4 * A;
double A_6 = A_5 * A;
double T = tan_phi * tan_phi;
double T_2 = T * T;
double C = (e_2 / (1 - e_2)) * cos_2_phi;
double k0 = 0.9996;
double E0 = 500000;
double N0 = GfrUtilConversionCylindricalUtm._s_getN0_(phiRad);
double s = GfrUtilConversionCylindrical.s_getMeridionalArc_precise(phiRad, e_2, a) / a;
double E = E0 + k0 * a * nu * (A + (1D - T + C) * (A_3 / 6D) + (5D - 18D * T + T_2) * (A_5 / 120D));
double N = N0 + k0 * a * (s + nu * tan_phi * (A_2 / 2D + (5D - T + 9D * C + 4 * C * C) * (A_4 / 24D) + (61D - 58D * T + T_2) * (A_6 / 720D)));
double prec = 100;
E = Math.round(E * prec) / prec;
N = Math.round(N * prec) / prec;
String strZoneLetter = GfrUtilConversionCylindricalUtm.s_getUtmZoneLetter(Math.toDegrees(phiRad));
int zone = GfrUtilConversionCylindricalUtm.s_getUtmZoneNumber(lambdaDeg);
return new GfrPointUtm(N, E, zone, strZoneLetter);
}
public static double s_getHemisphereFromLetterZoneUtm(String str) throws Exception
{
int intAscii = str.hashCode();
if (intAscii >= 67 && intAscii <= 88) //between C and X
{
if (intAscii < 78)// (A to M)
return -1;
return 1;
}
if (intAscii >= 99 && intAscii <= 120) //between c and x
{
if (intAscii < 110)// (a to m)
return -1;
return 1;
}
String strMessage = "Bad UTM zone letter : " + str;
throw new Exception(strMessage);
}
private static double _s_getN0_(GfrPointUtm pnt)
{
String strZoneLetter = pnt.getZoneLetter();
try
{
double dblHemisphere =
GfrUtilConversionCylindricalUtm.s_getHemisphereFromLetterZoneUtm(strZoneLetter);
if (dblHemisphere == -1)
return 10000000;
}
catch (Exception ex)
{
ex.printStackTrace();
}
return 0;
}
private static double _s_getN0_(double phi)
{
if (phi < 0)
return 10000000;
else
return 0;
}
private static double _s_getLambda0_(double lambda)
{
if (lambda > 0)
return (double) ((((int) lambda) / 6) * 6D + 3);
return (double) ((((int) lambda) / 6) * 6D - 3);
}
private static double _s_getLambda0_(int intZone)
{
double lambda0 = -180 + intZone * 6 - 3;
return (double) lambda0;
}
private static GfrPointDms _s_getDmsCoordinateFromUtmImpl_(
GfrPointUtm pnt,
GfrPrjAbs prj)// throws Exception
{
double x = pnt.getEasting();
double y = pnt.getNorthing();
x -= 500000;
y -= GfrUtilConversionCylindricalUtm._s_getN0_(pnt);
double a = prj.getEllipsoid().getDblRadius();
double e_2 = prj.getEllipsoid().getDblSquareExcentricity();
double e_4 = e_2 * e_2;
double e_6 = e_4 * e_2;
double k0 = 0.9996;
double M = y / k0;
double mu = M / (a * (1 - e_2 / 4 - 3 * e_4 / 64 - 5 * e_6 / 256));
double sin_2mu = Math.sin(2 * mu);
double sin_4mu = Math.sin(4 * mu);
double sin_6mu = Math.sin(6 * mu);
double sin_8mu = Math.sin(8 * mu);
double e1 = (1 - Math.sqrt(1 - e_2)) / (1 + Math.sqrt(1 - e_2));
double e1_2 = e1 * e1;
double e1_3 = e1_2 * e1;
double e1_4 = e1_3 * e1;
double J1 = (3 * e1 / 2 - 27 * e1_3 / 32);
double J2 = (21 * e1_2 / 16 - 55 * e1_4 / 32);
double J3 = (151 * e1_3 / 96);
double J4 = (1097 * e1_4 / 512);
double fp = mu + J1 * sin_2mu + J2 * sin_4mu + J3 * sin_6mu + J4 * sin_8mu;
double sin_fp = Math.sin(fp);
double cos_fp = Math.cos(fp);
double tan_fp = Math.tan(fp);
double cos_2_fp = cos_fp * cos_fp;
double sin_2_fp = sin_fp * sin_fp;
double e_prime_2 = e_2 / (1 - e_2);
double C1 = e_prime_2 * cos_2_fp;
double C1_2 = C1 * C1;
double T1 = tan_fp * tan_fp;
double T1_2 = T1 * T1;
double sqrt_exp = Math.sqrt(1 - e_2 * sin_2_fp);
double sqrt_3_exp = sqrt_exp * sqrt_exp * sqrt_exp;
double R1 = a * (1 - e_2) / (sqrt_3_exp);
double N1 = a / sqrt_exp;
double D = x / (N1 * k0);
double D_2 = D * D;
double D_3 = D_2 * D;
double D_4 = D_2 * D_2;
double D_5 = D_4 * D;
double D_6 = D_4 * D_2;
double Q1 = N1 * tan_fp / R1;
double Q2 = D_2 / 2;
double Q3 = (5 + 3 * T1 + 10 * C1 * 4 * C1_2 - 9 * e_prime_2) * D_4 / 24;
double Q4 = (61 + 90 * T1 + 298 * C1 + 45 * T1_2 - 3 * C1_2 - 252 * e_prime_2) * D_6 / 720;
double lat = fp - Q1 * (Q2 - Q3 + Q4);
int intZone = pnt.getZoneNumber();
double long0 = Math.toRadians(GfrUtilConversionCylindricalUtm._s_getLambda0_(intZone));
double Q5 = D;
double Q6 = (1 + 2 * T1 + C1) * D_3 / 6;
double Q7 = (5 - 2 * C1 + 28 * T1 - 3 * C1_2 + 8 * e_prime_2 + 24 * T1_2) * D_5 / 120;
double lon = long0 + (Q5 - Q6 + Q7) / cos_fp;
lat = Math.toDegrees(lat);
lon = Math.toDegrees(lon);
double prec = 0;
if (Math.abs(lat) < 85)
prec = 1000000D;
else if (Math.abs(lat) <= 89)
prec = 100000D;
else
prec = 10000D;
return new GfrPointDms(
Math.round(lon * prec) / prec,
Math.round(lat * prec) / prec);
}
public static GfrPointDms s_getDmsCoordinateFromUtm(GfrPointUtm pnt, GfrPrjAbs prj)// throws Exception
{
return _s_getDmsCoordinateFromUtmImpl_(pnt, prj);
}
}