Web墨卡托投影坐标转换的Golang实现

Web墨卡托投影(又称球形墨卡托投影)是Web地图(如Google Maps、OpenStreetMap等)常用的投影方式。下面我将介绍如何使用Golang实现Web墨卡托投影的坐标转换。

基本概念

Web墨卡托投影有以下特点：

• 将地球视为完美的球体
• 经度直接映射为X坐标
• 纬度使用墨卡托公式转换为Y坐标
• 坐标范围：经度[-180,180]，纬度[-85.05112878, 85.05112878]

核心实现代码

package mercator

import (
	"math"
)

const (
	earthRadius = 6378137.0 // WGS84椭球体长半轴
	maxLatitude = 85.05112877980659 // 最大纬度(弧度)
)

// LonLatToMeters 将经纬度转换为Web墨卡托投影坐标(单位：米)
func LonLatToMeters(lon, lat float64) (x, y float64) {
	// 确保纬度在有效范围内
	if lat > maxLatitude {
		lat = maxLatitude
	} else if lat < -maxLatitude {
		lat = -maxLatitude
	}

	x = lon * earthRadius * math.Pi / 180.0
	y = math.Log(math.Tan((90.0+lat)*math.Pi/360.0)) * earthRadius

	return x, y
}

// MetersToLonLat 将Web墨卡托投影坐标转换回经纬度
func MetersToLonLat(x, y float64) (lon, lat float64) {
	lon = (x / earthRadius) * 180.0 / math.Pi
	lat = (math.Atan(math.Exp(y/earthRadius)) * 360.0 / math.Pi) - 90.0

	return lon, lat
}

// LonLatToPixels 将经纬度转换为像素坐标(基于缩放级别和瓦片大小)
func LonLatToPixels(lon, lat float64, zoom int, tileSize int) (x, y float64) {
	xMeters, yMeters := LonLatToMeters(lon, lat)
	resolution := (2 * math.Pi * earthRadius) / (float64(tileSize) * math.Exp2(float64(zoom)))
	
	x = (xMeters + (math.Pi * earthRadius)) / resolution
	y = ((math.Pi * earthRadius) - yMeters) / resolution
	
	return x, y
}

// PixelsToLonLat 将像素坐标转换回经纬度
func PixelsToLonLat(x, y float64, zoom int, tileSize int) (lon, lat float64) {
	resolution := (2 * math.Pi * earthRadius) / (float64(tileSize) * math.Exp2(float64(zoom)))
	xMeters := x*resolution - (math.Pi * earthRadius)
	yMeters := (math.Pi * earthRadius) - y*resolution
	
	return MetersToLonLat(xMeters, yMeters)
}

使用示例

package main

import (
	"fmt"
	"yourmodulepath/mercator"
)

func main() {
	// 经纬度坐标(天安门)
	lon, lat := 116.391275, 39.906217

	// 转换为墨卡托坐标(米)
	x, y := mercator.LonLatToMeters(lon, lat)
	fmt.Printf("墨卡托坐标(米): x=%.2f, y=%.2f\n", x, y)

	// 转换回经纬度
	lon2, lat2 := mercator.MetersToLonLat(x, y)
	fmt.Printf("转换回经纬度: lon=%.6f, lat=%.6f\n", lon2, lat2)

	// 转换为像素坐标(缩放级别12, 瓦片大小256)
	zoom := 12
	tileSize := 256
	px, py := mercator.LonLatToPixels(lon, lat, zoom, tileSize)
	fmt.Printf("像素坐标: x=%.2f, y=%.2f\n", px, py)

	// 从像素坐标转换回经纬度
	lon3, lat3 := mercator.PixelsToLonLat(px, py, zoom, tileSize)
	fmt.Printf("从像素转换回经纬度: lon=%.6f, lat=%.6f\n", lon3, lat3)
}

扩展功能

如果需要更完整的功能，可以考虑添加以下方法：

// TileToLonLat 计算瓦片的左上角经纬度
func TileToLonLat(x, y, zoom int) (lon, lat float64) {
	n := math.Pi - 2.0*math.Pi*float64(y)/math.Exp2(float64(zoom))
	lon = float64(x)/math.Exp2(float64(zoom))*360.0 - 180.0
	lat = 180.0 / math.Pi * math.Atan(0.5*(math.Exp(n)-math.Exp(-n)))
	return lon, lat
}

// LonLatToTile 计算包含给定经纬度的瓦片坐标
func LonLatToTile(lon, lat float64, zoom int) (x, y int) {
	x = int(math.Floor((lon + 180.0) / 360.0 * math.Exp2(float64(zoom))))
	y = int(math.Floor((1.0 - math.Log(math.Tan(lat*math.Pi/180.0)+1.0/math.Cos(lat*math.Pi/180.0))/math.Pi) / 2.0 * math.Exp2(float64(zoom))))
	return x, y
}

注意事项

1. Web墨卡托投影在高纬度地区会有较大变形，因此纬度范围被限制在约±85.05度
2. 对于需要更高精度的应用，可以考虑使用椭球体模型而非球体模型
3. 在实际地图应用中，通常还需要处理瓦片坐标系统和像素坐标系统

这个实现提供了Web墨卡托投影的核心功能，可以满足大多数Web地图应用的需求。如果需要更复杂的功能，可以考虑在此基础上进行扩展。