Go语言中OSGB36东/北坐标到经纬度坐标的转换教程

package main

import (
    "fmt"
    "log"

    "github.com/pebbe/go-proj-4/proj" // 假设已安装此库
)

func main() {
    // 定义源投影：OSGB36国家格网 (EPSG:27700)
    // PROJ.4 字符串定义通常为 "+init=epsg:XXXX" 或详细参数
    // 这里使用EPSG代码，需要PROJ.4库支持
    sourceCRS := "+init=epsg:27700" // OSGB36 British National Grid
    // 定义目标投影：WGS84经纬度 (EPSG:4326)
    targetCRS := "+init=epsg:4326"  // WGS84 Latitude/Longitude

    // 创建转换上下文
    pj, err := proj.New(sourceCRS, targetCRS)
    if err != nil {
        log.Fatalf("创建投影转换器失败: %v", err)
    }
    defer pj.Close() // 确保在函数结束时关闭转换器

    // 输入的OSGB36东/北坐标 (米)
    easting := 348356.0
    northing := 862582.0

    // 执行转换
    // 注意：proj.Transform 返回的顺序通常是经度、纬度
    // 且单位是弧度，需要转换为度
    lonRad, latRad, err := pj.Transform(easting, northing)
    if err != nil {
        log.Fatalf("坐标转换失败: %v", err)
    }

    // 将弧度转换为十进制度
    lonDeg := proj.RadToDeg(lonRad)
    latDeg := proj.RadToDeg(latRad)

    fmt.Printf("OSGB36 Easting: %.0f, Northing: %.0f\n", easting, northing)
    fmt.Printf("转换为WGS84 Latitude: %.5f, Longitude: %.5f\n", latDeg, lonDeg)

    // 预期输出示例: Latitude: 41.40338, Longitude: 2.17403
    // 注意：这里的示例输出与原始问题中的示例输出不符，
    // 原始问题中的示例 (41.40338, 2.17403) 看起来更像是欧洲大陆的坐标，
    // 而OSGB36是英国格网。
    // 真实的OSGB36 (348356, 862582) 对应的WGS84坐标应在英国境内。
    // 例如：Easting: 348356, Northing: 862582 接近苏格兰高地。
    // 实际结果可能为 Lat: 57.77, Lon: -4.73 左右。
    // 请以实际运行结果为准。
}

package main

import (
    "fmt"
    "math"
)

// 定义坐标结构体
type LatLon struct {
    Latitude  float64 // 纬度，十进制度
    Longitude float64 // 经度，十进制度
}

type EastingNorthing struct {
    Easting  float64 // 东坐标，米
    Northing float64 // 北坐标，米
}

// OSGB36椭球体参数 (Airy 1830)
const (
    a_airy = 6377563.396 // 长半轴
    b_airy = 6356256.910 // 短半轴
    // 其他参数如扁率、偏心率等需要根据a_airy和b_airy计算
)

// WGS84椭球体参数
const (
    a_wgs84 = 6378137.0 // 长半轴
    b_wgs84 = 6356752.314245 // 短半轴
    // 其他参数
)

// Helmert七参数转换参数 (OSGB36 -> WGS84)
// 这些参数是标准值，需要准确获取
const (
    tx = -446.448 // X轴平移 (米)
    ty = 125.717  // Y轴平移 (米)
    tz = -542.008 // Z轴平移 (米)
    rx = -0.150   // X轴旋转 (弧度)
    ry = -0.247   // Y轴旋转 (弧度)
    rz = -0.842   // Z轴旋转 (弧度)
    s  = 20.489   // 尺度因子 (ppm)
)

// ConvertOSGB36ToWGS84 是主转换函数
// 它将OSGB36的东/北坐标转换为WGS84的经纬度坐标
func ConvertOSGB36ToWGS84(en EastingNorthing) (LatLon, error) {
    // 1. 将OSGB36格网坐标反投影到OSGB36大地坐标 (经纬度)
    // 这部分涉及高斯投影的逆运算，非常复杂，需要大量的数学公式
    // 例如，计算子午线弧长、子午线收敛角、投影带宽度等
    // 伪代码:
    // latOSGB36, lonOSGB36 := inverseGaussProjection(en.Easting, en.Northing, a_airy, b_airy, /* 其他投影参数 */)

    // 假设我们已经得到了OSGB36的大地坐标（这里用假值代替，实际需要计算）
    // 例如，对于输入 348356, 862582，对应的OSGB36大地坐标可能是：
    // latOSGB36 := 57.77 * math.Pi / 180 // 纬度，弧度
    // lonOSGB36 := -4.73 * math.Pi / 180 // 经度，弧度

    // 2. 将OSGB36大地坐标转换为地心直角坐标 (X, Y, Z)
    // x_osgb36, y_osgb36, z_osgb36 := geodeticToCartesian(latOSGB36, lonOSGB36, a_airy, b_airy)

    // 3. 应用Helmert七参数转换，将OSGB36地心坐标转换为WGS84地心坐标
    // x_wgs84, y_wgs84, z_wgs84 := helmertTransform(x_osgb36, y_osgb36, z_osgb36, tx, ty, tz, rx, ry, rz, s)

    // 4. 将WGS84地心坐标转换回WGS84大地坐标 (经纬度)
    // latWGS84, lonWGS84 := cartesianToGeodetic(x_wgs84, y_wgs84, z_wgs84, a_wgs84, b_wgs84)

    // 这里返回一个模拟结果，实际需要实现上述所有步骤
    // 原始问题中的示例输出 (41.40338, 2.17403) 并非英国坐标，
    // 而是巴塞罗那附近的坐标，这与OSGB36的地理范围不符。
    // 实际OSGB36 (348356, 862582) 对应的WGS84坐标大约在苏格兰高地。
    // 为了演示，这里假设转换成功并返回一个合理的值（但不是原始问题中的示例值）
    // 你需要根据实际移植的算法来计算。
    resultLatLon := LatLon{
        Latitude:  57.77, // 示例值，实际应通过计算获得
        Longitude: -4.73, // 示例值，实际应通过计算获得
    }

    return resultLatLon, nil
}

func main() {
    inputEN := EastingNorthing{
        Easting:  348356.0,
        Northing: 862582.0,
    }

    outputLL, err := ConvertOSGB36ToWGS84(inputEN)
    if err != nil {
        fmt.Printf("转换失败: %v\n", err)
        return
    }

    fmt.Printf("输入 OSGB36 Easting: %.0f, Northing: %.0f\n", inputEN.Easting, inputEN.Northing)
    fmt.Printf("输出 WGS84 Latitude: %.5f, Longitude: %.5f\n", outputLL.Latitude, outputLL.Longitude)
}

// 辅助函数（实际需要实现）
// inverseGaussProjection(easting, northing, a, b, ...) (lat, lon)
// geodeticToCartesian(lat, lon, a, b) (x, y, z)
// helmertTransform(x, y, z, tx, ty, tz, rx, ry, rz, s) (x_prime, y_prime, z_prime)
// cartesianToGeodetic(x, y, z, a, b) (lat, lon)