unity 实现经纬度和utm坐标转换代码

以下是 Unity 中实现经纬度和 UTM 坐标转换的代码：

using UnityEngine;

public class GeoCoordinateConverter : MonoBehaviour
{
    private const double EquatorialRadius = 6378137.0; // The equatorial radius of the earth in meters
    private const double PolarRadius = 6356752.3; // The polar radius of the earth in meters
    private const double Flattening = 0.00335281066474748; // The flattening of the earth

    // Convert latitude and longitude (in degrees) to UTM coordinates
    public static Vector2 LatLonToUTM(double latitude, double longitude)
    {
        double falseNorthing = latitude < 0 ? 10000000 : 0; // UTM false northing
        double zone = Mathf.Floor((float)(longitude + 180) / 6.0f) + 1; // UTM zone

        double latRad = latitude * Mathf.Deg2Rad;
        double lonRad = longitude * Mathf.Deg2Rad;

        double lonOrigin = (zone - 1) * 6 - 180 + 3; // UTM central meridian
        double lonOriginRad = lonOrigin * Mathf.Deg2Rad;

        double eccPrimeSquared = Flattening / (1 - Flattening); // Eccentricity prime squared
        double N = EquatorialRadius / Mathf.Sqrt(1 - Flattening * Mathf.Pow(Mathf.Sin((float)latRad), 2)); // Radius of curvature of the earth's surface
        double T = Mathf.Pow(Mathf.Tan((float)latRad), 2); // Tangent of latitude
        double C = eccPrimeSquared * Mathf.Pow(Mathf.Cos((float)latRad), 2); // Coefficient

        double A = (lonRad - lonOriginRad) * Mathf.Cos((float)latRad);
        double M = EquatorialRadius * ((1 - Flattening / 4 - 3 * Mathf.Pow(Flattening, 2) / 64 - 5 * Mathf.Pow(Flattening, 3) / 256) * latRad
            - (3 * Flattening / 8 + 3 * Mathf.Pow(Flattening, 2) / 32 + 45 * Mathf.Pow(Flattening, 3) / 1024) * Mathf.Sin(2 * (float)latRad)
            + (15 * Mathf.Pow(Flattening, 2) / 256 + 45 * Mathf.Pow(Flattening, 3) / 1024) * Mathf.Sin(4 * (float)latRad)
            - (35 * Mathf.Pow(Flattening, 3) / 3072) * Mathf.Sin(6 * (float)latRad));

        double easting = 500000 + (0.9996 * N * (A + (1 - T + C) * Mathf.Pow(A, 3) / 6 + (5 - 18 * T + Mathf.Pow(T, 2) + 72 * C - 58 * eccPrimeSquared) * Mathf.Pow(A, 5) / 120));
        double northing = falseNorthing + (0.9996 * (M + N * Mathf.Tan((float)latRad) * (Mathf.Pow(A, 2) / 2 + (5 - T + 9 * C + 4 * Mathf.Pow(C, 2)) * Mathf.Pow(A, 4) / 24 + (61 - 58 * T + Mathf.Pow(T, 2) + 600 * C - 330 * eccPrimeSquared) * Mathf.Pow(A, 6) / 720)));

        return new Vector2((float)easting, (float)northing);
    }

    // Convert UTM coordinates to latitude and longitude (in degrees)
    public static Vector2 UTMToLatLon(double easting, double northing, int zone, bool isNorth)
    {
        double falseNorthing = isNorth ? 0 : 10000000; // UTM false northing

        double a = EquatorialRadius;
        double eSquared = Flattening * (2 - Flattening);
        double ePrimeSquared = eSquared / (1 - eSquared);
        double e1 = (1 - Mathf.Sqrt(1 - (float)eSquared)) / (1 + Mathf.Sqrt(1 - (float)eSquared));

        double x = easting - 500000; // UTM easting
        double y = northing - falseNorthing; // UTM northing

        double lonOrigin = (zone - 1) * 6 - 180 + 3; // UTM central meridian
        double lonOriginRad = lonOrigin * Mathf.Deg2Rad;

        double M = y / 0.9996;
        double mu = M / (a * (1 - eSquared / 4 - 3 * Mathf.Pow(eSquared, 2) / 64 - 5 * Mathf.Pow(eSquared, 3) / 256));

        double phi1Rad = mu + (3 * e1 / 2 - 27 * Mathf.Pow(e1, 3) / 32) * Mathf.Sin(2 * mu)
            + (21 * Mathf.Pow(e1, 2) / 16 - 55 * Mathf.Pow(e1, 4) / 32) * Mathf.Sin(4 * mu)
            + (151 * Mathf.Pow(e1, 3) / 96) * Mathf.Sin(6 * mu)
            + (1097 * Mathf.Pow(e1, 4) / 512) * Mathf.Sin(8 * mu);

        double phi1 = phi1Rad * Mathf.Rad2Deg;

        double N1 = a / Mathf.Sqrt(1 - eSquared * Mathf.Pow(Mathf.Sin((float)phi1Rad), 2));
        double T1 = Mathf.Pow(Mathf.Tan((float)phi1Rad), 2);
        double C1 = ePrimeSquared * Mathf.Pow(Mathf.Cos((float)phi1Rad), 2);
        double R1 = a * (1 - eSquared) / Mathf.Pow(1 - eSquared * Mathf.Pow(Mathf.Sin((float)phi1Rad), 2), 1.5);
        double D = x / (N1 * 0.9996);

        double lat = phi1Rad - (N1 * Mathf.Tan((float)phi1Rad) / R1) *
            (Mathf.Pow(D, 2) / 2
            - (5 + 3 * T1 + 10 * C1 - 4 * Mathf.Pow(C1, 2) - 9 * ePrimeSquared) * Mathf.Pow(D, 4) / 24
            + (61 + 90 * T1 + 298 * C1 + 45 * Mathf.Pow(T1, 2) - 252 * ePrimeSquared - 3 * Mathf.Pow(C1, 2)) * Mathf.Pow(D, 6) / 720);

        double lon = lonOriginRad + (D -
            (1 + 2 * T1 + C1) * Mathf.Pow(D, 3) / 6
            + (5 - 2 * C1 + 28 * T1 - 3 * Mathf.Pow(C1, 2) + 8 * ePrimeSquared + 24 * Mathf.Pow(T1, 2)) * Mathf.Pow(D, 5) / 120) / Mathf.Cos((float)phi1Rad);

        return new Vector2((float)lat * Mathf.Rad2Deg, (float)lon * Mathf.Rad2Deg);
    }
}

使用示例：

// Convert latitude and longitude to UTM coordinates
Vector2 utmCoord = GeoCoordinateConverter.LatLonToUTM(39.914142, 116.403959);
Debug.Log("UTM coordinates: " + utmCoord.x + ", " + utmCoord.y);

// Convert UTM coordinates to latitude and longitude
Vector2 latLonCoord = GeoCoordinateConverter.UTMToLatLon(455207.0339, 4405140.4426, 50, true);
Debug.Log("Latitude and longitude: " + latLonCoord.x + ", " + latLonCoord.y);

注意，需要使用正确的地球半径和扁率来进行计算，以上代码中使用的是 WGS84 椭球体的参数。此外，UTM 转换需要指定正确的区域和北半球/南半球。