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

# 安装 PROJ.4 C库 (以Ubuntu为例)
sudo apt-get update
sudo apt-get install libproj-dev proj-bin

# 安装 go-proj-4 Go模块
go get github.com/pebbe/go-proj-4/proj

package main

import (
    "fmt"
    "log"

    "github.com/pebbe/go-proj-4/proj"
)

func main() {
    // 定义源坐标系统：OSGB36 (EPSG:27700)
    sourceCRS := "EPSG:27700"
    // 定义目标坐标系统：WGS84 经纬度 (EPSG:4326)
    targetCRS := "EPSG:4326"

    // 创建一个坐标转换器
    transformer, err := proj.NewCRSTransformer(sourceCRS, targetCRS)
    if err != nil {
        log.Fatalf("创建CRS转换器失败: %v", err)
    }
    defer transformer.Close() // 确保在函数结束时关闭转换器资源

    // 示例OSGB36坐标 (伦敦市中心附近)
    // Easting: 530000.0, Northing: 180000.0
    easting := 530000.0
    northing := 180000.0

    // 执行转换
    // Transform函数接收 (x, y, z)，对于2D坐标，z可以设为0
    // 返回值是 (longitude, latitude, z)
    lon, lat, err := transformer.Transform(easting, northing, 0)
    if err != nil {
        log.Fatalf("坐标转换失败: %v", err)
    }

    fmt.Printf("原始OSGB36坐标: 东向(Easting) %.2f, 北向(Northing) %.2f\n", easting, northing)
    fmt.Printf("转换后的WGS84坐标: 纬度(Latitude) %.6f, 经度(Longitude) %.6f\n", lat, lon)

    // --- 关于问题描述中示例的说明 ---
    fmt.Println("\n--- 关于问题描述中示例的说明 ---")
    fmt.Println("请注意，问题描述中提供的示例输入 (348356, 862582) 和输出 (41.40338, 2.17403) 似乎不属于OSGB36系统。")
    fmt.Println("输出的经纬度 (41.40338, 2.17403) 对应西班牙巴塞罗那地区，而OSGB36坐标系统主要覆盖英国。")
    fmt.Println("如果将 (348356, 862582) 作为OSGB36坐标进行转换，其结果将落在北海区域，而非巴塞罗那。")
    fmt.Println("因此，本教程的示例代码以真实的OSGB36坐标进行演示，以确保结果的准确性。")
}

原始OSGB36坐标: 东向(Easting) 530000.00, 北向(Northing) 180000.00
转换后的WGS84坐标: 纬度(Latitude) 51.513418, 经度(Longitude) -0.088629

package main

import (
    "fmt"
    "math"
)

// 定义椭球体参数
type Ellipsoid struct {
    A float64 // 长半轴
    B float64 // 短半轴
    F float64 // 扁率
    E2 float64 // 第一偏心率的平方
}

// OSGB36 (Airy 1830) 椭球体参数
var airy1830 = Ellipsoid{
    A: 6377563.396,
    B: 6356256.910,
    F: 1 / 299.3249646,
}

// WGS84 椭球体参数
var wgs84 = Ellipsoid{
    A: 6378137.0,
    B: 6356752.314245,
    F: 1 / 298.257223563,
}

func init() {
    airy1830.E2 = (airy1830.A*airy1830.A - airy1830.B*airy1830.B) / (airy1830.A*airy1830.A)
    wgs84.E2 = (wgs84.A*wgs84.A - wgs84.B*wgs84.B) / (wgs84.A*wgs84.A)
}

// 定义Helmert变换参数 (OSGB36 到 WGS84)
// 这些参数需要精确定义，通常来自权威机构
type HelmertParams struct {
    Tx, Ty, Tz float64 // 平移 (m)
    Rx, Ry, Rz float64 // 旋转 (弧度)
    S          float64 // 尺度因子
}

// 这是一个简化的Helmert参数示例，实际值需要查阅资料
var osgb36ToWgs84Helmert = HelmertParams{
    Tx: -446.448, Ty: 125.157, Tz: -542.060,
    Rx: -0.1502 * math.Pi / (180 * 3600), // 转换为弧度
    Ry: -0.2470 * math.Pi / (180 * 3600),
    Rz: -0.8421 * math.Pi / (180 * 3600),
    S