/*
* Copyright (C) 2011 United States Government as represented by the Administrator of the
* National Aeronautics and Space Administration.
* All Rights Reserved.
*/
package gov.nasa.worldwind.geom.coords;
import gov.nasa.worldwind.avlist.AVKey;
import gov.nasa.worldwind.geom.*;
/**
* Converter used to translate UTM coordinates to and from geodetic latitude and longitude.
*
* @author Patrick Murris
* @version $Id$
* @see UTMCoord, TMCoordConverter
*/
/**
* Ported to Java from the NGA GeoTrans utm.c and utm.h
*
* @author Garrett Headley, Patrick Murris
*/
class UTMCoordConverter
{
public final static double CLARKE_A = 6378206.4;
public final static double CLARKE_B = 6356583.8;
public final static double CLARKE_F = 1 / 294.9786982;
public final static double WGS84_A = 6378137;
public final static double WGS84_F = 1 / 298.257223563;
public final static int UTM_NO_ERROR = 0x0000;
public final static int UTM_LAT_ERROR = 0x0001;
public final static int UTM_LON_ERROR = 0x0002;
public final static int UTM_EASTING_ERROR = 0x0004;
public final static int UTM_NORTHING_ERROR = 0x0008;
public final static int UTM_ZONE_ERROR = 0x0010;
public final static int UTM_HEMISPHERE_ERROR = 0x0020;
public final static int UTM_ZONE_OVERRIDE_ERROR = 0x0040;
public final static int UTM_A_ERROR = 0x0080;
public final static int UTM_INV_F_ERROR = 0x0100;
public final static int UTM_TM_ERROR = 0x0200;
private final static double PI = 3.14159265358979323;
//private final static double MIN_LAT = ((-80.5 * PI) / 180.0); /* -80.5 degrees in radians */
//private final static double MAX_LAT = ((84.5 * PI) / 180.0); /* 84.5 degrees in radians */
private final static double MIN_LAT = ((-82 * PI) / 180.0); /* -82 degrees in radians */
private final static double MAX_LAT = ((86 * PI) / 180.0); /* 86 degrees in radians */
private final static int MIN_EASTING = 100000;
private final static int MAX_EASTING = 900000;
private final static int MIN_NORTHING = 0;
private final static int MAX_NORTHING = 10000000;
private double UTM_a = 6378137.0; /* Semi-major axis of ellipsoid in meters */
private double UTM_f = 1 / 298.257223563; /* Flattening of ellipsoid */
private long UTM_Override = 0; /* Zone override flag */
private double Easting;
private double Northing;
private String Hemisphere;
private int Zone;
private double Latitude;
private double Longitude;
private double Central_Meridian;
UTMCoordConverter(){}
UTMCoordConverter(double a, double f)
{
setUTMParameters(a, f, 0);
}
/**
* The function Set_UTM_Parameters receives the ellipsoid parameters and UTM zone override parameter as inputs, and
* sets the corresponding state variables. If any errors occur, the error code(s) are returned by the function,
* otherwise UTM_NO_ERROR is returned.
*
* @param a Semi-major axis of ellipsoid, in meters
* @param f Flattening of ellipsoid
* @param override UTM override zone, zero indicates no override
*
* @return error code
*/
private long setUTMParameters(double a, double f, long override)
{
double inv_f = 1 / f;
long Error_Code = UTM_NO_ERROR;
if (a <= 0.0)
{ /* Semi-major axis must be greater than zero */
Error_Code |= UTM_A_ERROR;
}
if ((inv_f < 250) || (inv_f > 350))
{ /* Inverse flattening must be between 250 and 350 */
Error_Code |= UTM_INV_F_ERROR;
}
if ((override < 0) || (override > 60))
{
Error_Code |= UTM_ZONE_OVERRIDE_ERROR;
}
if (Error_Code == UTM_NO_ERROR)
{ /* no errors */
UTM_a = a;
UTM_f = f;
UTM_Override = override;
}
return (Error_Code);
}
/**
* The function Convert_Geodetic_To_UTM converts geodetic (latitude and longitude) coordinates to UTM projection
* (zone, hemisphere, easting and northing) coordinates according to the current ellipsoid and UTM zone override
* parameters. If any errors occur, the error code(s) are returned by the function, otherwise UTM_NO_ERROR is
* returned.
*
* @param Latitude Latitude in radians
* @param Longitude Longitude in radians
*
* @return error code
*/
public long convertGeodeticToUTM(double Latitude, double Longitude)
{
long Lat_Degrees;
long Long_Degrees;
long temp_zone;
long Error_Code = UTM_NO_ERROR;
double Origin_Latitude = 0;
double False_Easting = 500000;
double False_Northing = 0;
double Scale = 0.9996;
if ((Latitude < MIN_LAT) || (Latitude > MAX_LAT))
{ /* Latitude out of range */
Error_Code |= UTM_LAT_ERROR;
}
if ((Longitude < -PI) || (Longitude > (2 * PI)))
{ /* Longitude out of range */
Error_Code |= UTM_LON_ERROR;
}
if (Error_Code == UTM_NO_ERROR)
{ /* no errors */
if (Longitude < 0)
Longitude += (2 * PI) + 1.0e-10;
Lat_Degrees = (long) (Latitude * 180.0 / PI);
Long_Degrees = (long) (Longitude * 180.0 / PI);
if (Longitude < PI)
temp_zone = (long) (31 + ((Longitude * 180.0 / PI) / 6.0));
else
temp_zone = (long) (((Longitude * 180.0 / PI) / 6.0) - 29);
if (temp_zone > 60)
temp_zone = 1;
/* UTM special cases */
if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > -1) && (Long_Degrees < 3))
temp_zone = 31;
if ((Lat_Degrees > 55) && (Lat_Degrees < 64) && (Long_Degrees > 2) && (Long_Degrees < 12))
temp_zone = 32;
if ((Lat_Degrees > 71) && (Long_Degrees > -1) && (Long_Degrees < 9))
temp_zone = 31;
if ((Lat_Degrees > 71) && (Long_Degrees > 8) && (Long_Degrees < 21))
temp_zone = 33;
if ((Lat_Degrees > 71) && (Long_Degrees > 20) && (Long_Degrees < 33))
temp_zone = 35;
if ((Lat_Degrees > 71) && (Long_Degrees > 32) && (Long_Degrees < 42))
temp_zone = 37;
if (UTM_Override != 0)
{
if ((temp_zone == 1) && (UTM_Override == 60))
temp_zone = UTM_Override;
else if ((temp_zone == 60) && (UTM_Override == 1))
temp_zone = UTM_Override;
else if (((temp_zone - 1) <= UTM_Override) && (UTM_Override <= (temp_zone + 1)))
temp_zone = UTM_Override;
else
Error_Code = UTM_ZONE_OVERRIDE_ERROR;
}
if (Error_Code == UTM_NO_ERROR)
{
if (temp_zone >= 31)
Central_Meridian = (6 * temp_zone - 183) * PI / 180.0;
else
Central_Meridian = (6 * temp_zone + 177) * PI / 180.0;
Zone = (int) temp_zone;
if (Latitude < 0)
{
False_Northing = 10000000;
Hemisphere = AVKey.SOUTH;
}
else
Hemisphere = AVKey.NORTH;
try
{
TMCoord TM = TMCoord.fromLatLon(Angle.fromRadians(Latitude), Angle.fromRadians(Longitude),
this.UTM_a, this.UTM_f, Angle.fromRadians(Origin_Latitude),
Angle.fromRadians(Central_Meridian), False_Easting, False_Northing, Scale);
Easting = TM.getEasting();
Northing = TM.getNorthing();
if ((Easting < MIN_EASTING) || (Easting > MAX_EASTING))
Error_Code = UTM_EASTING_ERROR;
if ((Northing < MIN_NORTHING) || (Northing > MAX_NORTHING))
Error_Code |= UTM_NORTHING_ERROR;
}
catch (Exception e)
{
Error_Code = UTM_TM_ERROR;
}
}
}
return (Error_Code);
}
/** @return Easting (X) in meters */
public double getEasting()
{
return Easting;
}
/** @return Northing (Y) in meters */
public double getNorthing()
{
return Northing;
}
/**
* @return The coordinate hemisphere, either {@link gov.nasa.worldwind.avlist.AVKey#NORTH} or {@link
* gov.nasa.worldwind.avlist.AVKey#SOUTH}.
*/
public String getHemisphere()
{
return Hemisphere;
}
/** @return UTM zone */
public int getZone()
{
return Zone;
}
/**
* The function Convert_UTM_To_Geodetic converts UTM projection (zone, hemisphere, easting and northing) coordinates
* to geodetic(latitude and longitude) coordinates, according to the current ellipsoid parameters. If any errors
* occur, the error code(s) are returned by the function, otherwise UTM_NO_ERROR is returned.
*
* @param Zone UTM zone.
* @param Hemisphere The coordinate hemisphere, either {@link gov.nasa.worldwind.avlist.AVKey#NORTH} or {@link
* gov.nasa.worldwind.avlist.AVKey#SOUTH}.
* @param Easting Easting (X) in meters.
* @param Northing Northing (Y) in meters.
*
* @return error code.
*/
public long convertUTMToGeodetic(long Zone, String Hemisphere, double Easting, double Northing)
{
// TODO: arg checking
long Error_Code = UTM_NO_ERROR;
double Origin_Latitude = 0;
double False_Easting = 500000;
double False_Northing = 0;
double Scale = 0.9996;
if ((Zone < 1) || (Zone > 60))
Error_Code |= UTM_ZONE_ERROR;
if (!Hemisphere.equals(AVKey.SOUTH) && !Hemisphere.equals(AVKey.NORTH))
Error_Code |= UTM_HEMISPHERE_ERROR;
// if ((Easting < MIN_EASTING) || (Easting > MAX_EASTING)) //removed check to enable reprojecting images
// Error_Code |= UTM_EASTING_ERROR; //that extend into another zone
if ((Northing < MIN_NORTHING) || (Northing > MAX_NORTHING))
Error_Code |= UTM_NORTHING_ERROR;
if (Error_Code == UTM_NO_ERROR)
{ /* no errors */
if (Zone >= 31)
Central_Meridian = ((6 * Zone - 183) * PI / 180.0 /*+ 0.00000005*/);
else
Central_Meridian = ((6 * Zone + 177) * PI / 180.0 /*+ 0.00000005*/);
if (Hemisphere.equals(AVKey.SOUTH))
False_Northing = 10000000;
try
{
TMCoord TM = TMCoord.fromTM(Easting, Northing,
Angle.fromRadians(Origin_Latitude), Angle.fromRadians(Central_Meridian),
False_Easting, False_Northing, Scale);
Latitude = TM.getLatitude().radians;
Longitude = TM.getLongitude().radians;
if ((Latitude < MIN_LAT) || (Latitude > MAX_LAT))
{ /* Latitude out of range */
Error_Code |= UTM_NORTHING_ERROR;
}
}
catch (Exception e)
{
Error_Code = UTM_TM_ERROR;
}
}
return (Error_Code);
}
/** @return Latitude in radians. */
public double getLatitude()
{
return Latitude;
}
/** @return Longitude in radians. */
public double getLongitude()
{
return Longitude;
}
/** @return Central_Meridian in radians. */
public double getCentralMeridian()
{
return Central_Meridian;
}
public static LatLon convertWGS84ToNAD27(Angle latWGS, Angle lonWGS)
{
double deltaX = -12.0;
double deltaY = 130.0;
double deltaZ = 190.0;
double difA = WGS84_A - CLARKE_A;
double difF = WGS84_F - CLARKE_F;
double lat = latWGS.radians;
double lon = lonWGS.radians;
double f = 1 - CLARKE_B / CLARKE_A;
double e2 = 2 * f - Math.pow(f, 2);
double Rn = CLARKE_A / Math.sqrt(1 - e2 * Math.pow(Math.sin(lat), 2.0));
double Rm = (CLARKE_A * (1 - e2)) / Math.pow(1 - e2 * Math.pow(Math.sin(lat), 2.0), 1.5);
double errLon = (-1 * deltaX * Math.sin(lon) + deltaY * Math.cos(lon)) / (Rn * Math.cos(lat));
double errLat = (-1 * deltaX * Math.sin(lat) * Math.cos(lon) - deltaY * Math.sin(lat) * Math.sin(lon)
+ deltaZ * Math.cos(lat)
+ difA * (Rn * e2 * Math.sin(lat) * Math.cos(lat)) / CLARKE_A
+ difF * (Rm * CLARKE_A / CLARKE_B + Rn * CLARKE_B / CLARKE_A) * Math.sin(lat) * Math.cos(lat)) / Rm;
return LatLon.fromRadians(lat - errLat, lon - errLon);
}
}