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

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

using UnityEngine;

public class GeoUtils : MonoBehaviour
{
    // WGS84椭球体参数
    private const double WGS84_A = 6378137.0; // 长半轴
    private const double WGS84_B = 6356752.314245; // 短半轴
    private const double WGS84_F = 1 / 298.257223563; // 扁率
    private const double WGS84_E2 = 0.00669437999014; // 第一偏心率的平方
    private const double WGS84_E2C = 0.00673949674228; // 第二偏心率的平方

    // UTM 坐标系参数
    private const double UTM_K0 = 0.9996; // UTM 投影的比例因子
    private const double UTM_FE = 500000; // UTM 起始点的东偏移量
    private const double UTM_FN_N = 0; // 北半球的起始点纬度
    private const double UTM_FN_S = 10000000; // 南半球的起始点纬度

    // 经纬度转 UTM 坐标
    public static Vector2 LatLonToUTM(double lat, double lon)
    {
        double a = WGS84_A;
        double b = WGS84_B;
        double f = WGS84_F;
        double e2 = WGS84_E2;
        double e2c = WGS84_E2C;
        double k0 = UTM_K0;
        double fe = UTM_FE;
        double fn = lat < 0 ? UTM_FN_S : UTM_FN_N;

        double lonRad = lon * Mathf.Deg2Rad;
        double latRad = lat * Mathf.Deg2Rad;

        double utmZone = Mathf.Floor((float)((lon + 180.0) / 6.0)) + 1.0;

        double lonOrigin = (utmZone - 1) * 6 - 180 + 3;
        double lonOriginRad = lonOrigin * Mathf.Deg2Rad;

        double e2p = (Mathf.Pow(a, 2) - Mathf.Pow(b, 2)) / Mathf.Pow(b, 2);
        double n = a / Mathf.Sqrt(1 - e2 * Mathf.Pow(Mathf.Sin(latRad), 2));
        double t = Mathf.Tan(latRad) * Mathf.Tan(latRad);
        double c = e2p * Mathf.Pow(Mathf.Cos(latRad), 2);
        double a1 = Mathf.Cos(latRad) * (lonRad - lonOriginRad);
        double a2 = Mathf.Cos(latRad) * Mathf.Pow(a1, 3) * (1 / 6.0);
        double a3 = Mathf.Cos(latRad) * Mathf.Pow(a1, 5) * (1 / 120.0);
        double a4 = Mathf.Cos(latRad) * Mathf.Pow(a1, 7) * (1 / 5040.0);
        double m = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latRad -
            (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Mathf.Sin(2 * latRad) +
            (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Mathf.Sin(4 * latRad) -
            (35 * e2 * e2 * e2 / 3072) * Mathf.Sin(6 * latRad));

        double utmEasting = k0 * n * (a1 + (1 - t + c) * a2 + (5 - 18 * t + t * t + 72 * c - 58 * e2p) * a3 + (61 - 479 * t + 179 * Mathf.Pow(t, 2) - Mathf.Pow(t, 3)) * a4) + fe;
        double utmNorthing = k0 * (m + n * Mathf.Tan(latRad) * ((a1 * a1 * 0.5) + (1 - t + 2 * c) * Mathf.Pow(a1, 4) * (1 / 24.0) + (5 - 4 * t + 14 * c - 28 * e2p) * Mathf.Pow(a1, 6) * (1 / 720.0))) + fn;

        return new Vector2((float)utmEasting, (float)utmNorthing);
    }

    // UTM 坐标转经纬度
    public static Vector2 UTMToLatLon(double easting, double northing, int zone, bool isSouth)
    {
        double a = WGS84_A;
        double b = WGS84_B;
        double f = WGS84_F;
        double e2 = WGS84_E2;
        double e2c = WGS84_E2C;
        double k0 = UTM_K0;
        double fe = UTM_FE;
        double fn = isSouth ? UTM_FN_S : UTM_FN_N;

        double e1 = (1 - Mathf.Sqrt(1 - e2)) / (1 + Mathf.Sqrt(1 - e2));
        double x = easting - fe;
        double y = northing - fn;
        double m = y / k0;
        double sm = (a + b) / 2 * (1 + Mathf.Pow(e1, 2) / 4 + Mathf.Pow(e1, 4) / 64) * Mathf.Rad2Deg;
        double n = (a - b) / (a + b);
        double nu = a / Mathf.Sqrt(1 - e2 * Mathf.Pow(Mathf.Sin(sm * Mathf.Deg2Rad), 2));
        double p = x / (nu * k0);
        double p2 = p * p;
        double p3 = p2 * p;
        double p4 = p3 * p;
        double p5 = p4 * p;
        double p6 = p5 * p;
        double m1 = m / nu;
        double m2 = m1 * m1;
        double m3 = m2 * m1;
        double m4 = m3 * m1;
        double m5 = m4 * m1;
        double m6 = m5 * m1;
        double s1 = nu * e2 * Mathf.Sin(sm * Mathf.Deg2Rad) * Mathf.Cos(sm * Mathf.Deg2Rad) * Mathf.Pow(p, 2) * 0.5;
        double s2 = nu * e2 * Mathf.Sin(sm * Mathf.Deg2Rad) * Mathf.Cos(sm * Mathf.Deg2Rad) * Mathf.Pow(p, 4) * (1 / 24.0);
        double s3 = nu * e2 * Mathf.Sin(sm * Mathf.Deg2Rad) * Mathf.Cos(sm * Mathf.Deg2Rad) * Mathf.Pow(p, 6) * (1 / 720.0);
        double s4 = nu * e2 * Mathf.Sin(sm * Mathf.Deg2Rad) * Mathf.Cos(sm * Mathf.Deg2Rad) * Mathf.Pow(p, 8) * (1 / 40320.0);
        double s5 = nu * e2 * Mathf.Sin(sm * Mathf.Deg2Rad) * Mathf.Cos(sm * Mathf.Deg2Rad) * Mathf.Pow(p, 10) * (1 / 3628800.0);
        double s6 = nu * e2 * Mathf.Sin(sm * Mathf.Deg2Rad) * Mathf.Cos(sm * Mathf.Deg2Rad) * Mathf.Pow(p, 12) * (1 / 479001600.0);
        double lat = sm - m1 + (m2 - s1) * 2 + (m3 - s2) * 4 + (m4 - s3) * 6 + (m5 - s4) * 8 + (m6 - s5) * 10;
        double l1 = p - (1 + 2 * p2 + p4) * Mathf.Pow(n, 2) * m1 * 0.5 + (5 + 28 * p2 + 24 * p4 + 6 * p6 + 8 * n * n * p2) * Mathf.Pow(n, 4) * m1 * (1 / 24.0) -
            (61 + 662 * p2 + 1320 * p4 + 720 * p6) * Mathf.Pow(n, 6) * m1 * (1 / 720.0);
        double l2 = (1 + 2 * p2 + p4) * Mathf.Pow(n, 2) * l1 * 0.5 - (1 + 28 * p2 + 24 * p4 + 6 * p6 + 8 * n * n * p2) * Mathf.Pow(n, 4) * l1 * (1 / 24.0) +
            (61 + 662 * p2 + 1320 * p4 + 720 * p6) * Mathf.Pow(n, 6) * l1 * (1 / 720.0);
        double lon = zone * 6 - 183 + l2 / Mathf.Cos(sm * Mathf.Deg2Rad);

        return new Vector2((float)lat, (float)lon);
    }
}

其中，`LatLonToUTM` 方法接收一个纬度和经度，返回对应的 UTM 坐标系下的坐标。`UTMToLatLon` 方法接收一个 UTM 坐标系下的坐标和所在区域，返回对应的经纬度。需要注意的是，UTM 坐标系下的南半球起始点纬度为 10000000。