Python生成栅格瓦片叠加EPSG:3395海图的偏移问题排查

我正在开展一个基于EPSG:3395坐标系在海图上可视化船舶位置数据的项目。现有一批船舶的WGS84（EPSG:4326）经纬度坐标，需将其渲染为绿色菱形点生成瓦片，叠加至EPSG:3395海图。

执行步骤：

• 从JSON文件加载船舶位置数据；
• 将坐标从EPSG:4326转换至EPSG:3395；
• 生成不同缩放级别的瓦片，绘制绿色菱形标记船舶位置；
• 将生成的瓦片叠加到EPSG:3395海图上。

遇到的问题：叠加到EPSG:3395海图时存在明显偏移，但叠加到标准Web墨卡托投影的Web地图时对齐正常。使用Geoserver渲染瓦片则显示正常。

Python瓦片生成代码：

import os
import time
import cv2
import numpy as np
import orjson
from geo_tool import tile as tile_util
from pyproj import CRS, Transformer
from rtree import index

TILE_SIZE = 256
PADDING_PIXEL = 10

proj_wgs84 = CRS.from_epsg(4326)
proj_3395 = CRS.from_epsg(3395)
transformer_to_3395 = Transformer.from_crs(proj_wgs84, proj_3395, always_xy=True)
transformer_to_wgs84 = Transformer.from_crs(proj_3395, proj_wgs84, always_xy=True)


def draw_diamond(image, center, size, color):
    gap = 0
    points = np.array(
        [
            [center[0] + gap, center[1] - size + gap],
            [center[0] - size + gap, center[1] + gap],
            [center[0] + gap, center[1] + size + gap],
            [center[0] + size + gap, center[1] + gap],
        ]
    )
    cv2.fillPoly(image, [points], color)


def generate_tiles_for_zoom(zoom_level):
    max_tile = 2**zoom_level - 1
    return [
        (x, y, zoom_level) for x in range(max_tile + 1) for y in range(max_tile + 1)
    ]


def generate_tiles_for_zoom_range(zoom_level_start, zoom_level_end):
    return [
        generate_tiles_for_zoom(zoom_level)
        for zoom_level in range(zoom_level_start, zoom_level_end + 1)
    ]


def generate_map_tiles_cv(in_range_latlngs, x, y, z, bbox):
    tile_dir = "tiles"
    min_x, min_y, max_x, max_y = bbox
    tile_filename = f"{tile_dir}/{z}/{y}/{x}.png"
    extended_bbox = (
        min_x - (max_x - min_x) / TILE_SIZE * PADDING_PIXEL,
        min_y - (max_y - min_y) / TILE_SIZE * PADDING_PIXEL,
        max_x + (max_x - min_x) / TILE_SIZE * PADDING_PIXEL,
        max_y + (max_y - min_y) / TILE_SIZE * PADDING_PIXEL,
    )
    image = np.zeros(
        (TILE_SIZE + PADDING_PIXEL * 2, TILE_SIZE + PADDING_PIXEL * 2, 4),
        dtype=np.uint8,
    )
    positions = [
        (
            int(
                (latlng["x"] - extended_bbox[0])
                * (TILE_SIZE + PADDING_PIXEL * 2)
                / (extended_bbox[2] - extended_bbox[0])
            ),
            int(
                (extended_bbox[3] - latlng["y"])
                * (TILE_SIZE + PADDING_PIXEL * 2)
                / (extended_bbox[3] - extended_bbox[1])
            ),
        )
        for latlng in in_range_latlngs
    ]
    for px, py in positions:
        draw_diamond(image, (px, py), 2, (0, 255, 0, 255))
    if not os.path.exists(os.path.dirname(tile_filename)):
        os.makedirs(os.path.dirname(tile_filename))
    image = image[
        PADDING_PIXEL : TILE_SIZE + PADDING_PIXEL,
        PADDING_PIXEL : TILE_SIZE + PADDING_PIXEL,
    ]
    cv2.imwrite(tile_filename, image)


def main():
    with open("./position_7.json", "rb") as f:
        data = [
            {"x": x, "y": y}
            for latlng in orjson.loads(f.read())
            for x, y in [
                transformer_to_3395.transform(
                    float(latlng["lon"]), float(latlng["lat"])
                )
            ]
        ]

    tiles = [
        item for sublist in generate_tiles_for_zoom_range(3, 10) for item in sublist
    ]

    idx = index.Index()
    for i, latlng in enumerate(data):
        idx.insert(i, (latlng["x"], latlng["y"], latlng["x"], latlng["y"]))

    for x, y, z in tiles:
        bbox_wgs84 = tile_util.tile_to_bbox(x, y, z)
        min_x, min_y = transformer_to_3395.transform(bbox_wgs84[0], bbox_wgs84[1])
        max_x, max_y = transformer_to_3395.transform(bbox_wgs84[2], bbox_wgs84[3])
        bbox = (min_x, min_y, max_x, max_y)

        extended_bbox = (
            min_x - (max_x - min_x) / TILE_SIZE * PADDING_PIXEL,
            min_y - (max_y - min_y) / TILE_SIZE * PADDING_PIXEL,
            max_x + (max_x - min_x) / TILE_SIZE * PADDING_PIXEL,
            max_y + (max_y - min_y) / TILE_SIZE * PADDING_PIXEL,
        )
        rst = [data[i] for i in idx.intersection(extended_bbox)]

        if rst:
            generate_map_tiles_cv(rst, x, y, z, bbox)


if __name__ == "__main__":
    main()

Leaflet渲染代码：

<!DOCTYPE html>
<html lang="en">
<head>
    <meta charset="UTF-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <title>Leaflet Map</title>
    <!-- Leaflet CSS -->
    <link rel="stylesheet" href="https://unpkg.com/leaflet/dist/leaflet.css" />
    <style>
        #app {
            height: 100vh;
        }
    </style>
</head>
<body>
    <div id="app"></div>

    <!-- Leaflet JS -->
    <script src="https://unpkg.com/leaflet/dist/leaflet.js"></script>
    <script>
        const map = L.map('app', {
            crs: L.CRS.EPSG3395,
            layers: [
                new L.TileLayer('https://m12.shipxy.com/tile.c?l=Na&m=o&x={x}&y={y}&z={z}&'),
            ],
            zoomAnimation: false,
        });

        map.setView([32.1103, 120.0924], 10);

        new L.TileLayer("http://127.0.0.1:8081/tiles/{z}/{y}/{x}.png").addTo(map);
    </script>
</body>
</html>

疑问：

1. 叠加EPSG:3395海图时出现偏移的原因是什么？
2. 我的坐标转换或像素位置计算是否存在问题？
3. 使用EPSG:3395投影处理海图时需注意哪些特殊事项？

解答 #

1. 偏移原因分析

核心问题是瓦片Y轴原点定义不匹配：

• 标准Web墨卡托（EPSG:3857）的Y轴原点在南极，而你使用的EPSG:3395海图瓦片采用北极作为Y轴原点，即瓦片Y号是倒序排列的。
• 你生成瓦片时沿用了Web墨卡托的Y轴逻辑，导致叠加后出现上下偏移。

2. 坐标转换与像素计算问题排查

• 坐标转换部分：Transformer.from_crs(proj_wgs84, proj_3395, always_xy=True) 是正确的，EPSG:3395为东-北坐标系，符合always_xy=True的参数要求。
• 像素位置计算问题：
  1. 瓦片Y号反向：生成瓦片时的y值是Web墨卡托逻辑，但海图的Y号应为2^z - 1 - y，需在生成或加载时反转。
  2. 瓦片边界计算偏差：tile_util.tile_to_bbox(x, y, z) 大概率是基于Web墨卡托实现的，先转WGS84再转EPSG:3395会因投影范围差异导致bbox不准确，应直接基于EPSG:3395投影计算瓦片边界。
  3. 像素Y轴计算：(extended_bbox[3] - latlng["y"]) 逻辑正确（图像Y轴向下，地理Y轴向上），但瓦片Y号反向会导致双重反转，加剧偏移。

3. EPSG:3395投影处理注意事项

• 瓦片Y轴方向：EPSG:3395瓦片服务存在两种Y轴定义（北极/南极原点），必须确认所用海图的Y轴规则，生成瓦片时严格匹配。
• 投影范围：EPSG:3395是全球墨卡托投影，但部分海图服务会裁剪高纬度区域，生成瓦片时需限定在海图有效范围内。
• 坐标单位：EPSG:3395单位为米，计算像素位置时要保证bbox数值精度，避免浮点运算误差累积。
• 边界计算方式：直接使用EPSG:3395的瓦片投影公式计算边界，不要通过WGS84中转，确保准确性。

内容的提问来源于stack exchange，提问作者J.Crax