如何基于WMS的BBOX参数计算CRS:25832下的地图瓦片四角点？

我想用PHP编写一个小型WMS服务，从MariaDB数据库查询多边形数据并配合Leaflet使用。Leaflet.js结合proj4leaflet.js会向该服务传递BBOX参数，提供西北点（X1,Y1）和东南点（X2,Y2）的坐标，WMS服务需要据此计算请求覆盖的区域形状。

在常规直角坐标系中，可通过以下四个点计算区域：

X1,Y1  X2,Y1  X2,Y2  X1,Y2

但这种方法无法适配CRS:25832的“倾斜”特性，导致网格出现恼人的缝隙。

我尝试用数学方法求解东北点，伪代码如下：

DIST = DISTANCE(POINT(X1,Y1), POINT(X2,Y2)) / SQRT(2);
NE = POINT(X1 + DIST, Y2 + DIST)

（西南点也用相同方法计算）。生成的网格缝隙已大幅缩小，但仍不够完美。

请问如何计算地图瓦片菱形区域的实际四角点？理想情况下可完全通过SQL虚拟列实现，PHP方案也能接受。

补充说明 #

方块1、2、3和4本应完美对齐，但实际并未对齐：

仔细观察会发现，1号方块的东南点与4号方块的西北点（粉色）因Leaflet的标识而对齐：

2号方块的西南点（紫色）与3号方块的东北点（蓝色）本应完全重合，但实际并未重合。

核心逻辑 #

CRS:25832（UTM Zone 32N）本身是直角坐标系，但Leaflet结合proj4leaflet渲染时，瓦片投影转换的「网格对齐」要求，导致常规矩形计算的区域和实际覆盖的菱形区域不匹配。正确做法是基于UTM网格规则，将BBOX边界对齐到网格线，确保相邻瓦片边缘完全重合。

方法1：MariaDB SQL实现 #

可直接在查询中计算，或创建虚拟列复用逻辑：

-- 传入参数：@x1（西北X）、@y1（西北Y）、@x2（东南X）、@y2（东南Y）
-- 1. 设置当前瓦片层级的UTM分辨率（示例为Zoom 18，德国北部区域参考值）
SET @resolution = 0.598;

-- 2. 计算对齐到网格线的东北、西南点坐标
SET @ne_x = FLOOR(@x2 / @resolution) * @resolution;
SET @sw_y = FLOOR(@y1 / @resolution) * @resolution;

-- 3. 获取四个顶点
SELECT
  ST_Point(@x1, @y1) AS nw,
  ST_Point(@ne_x, @y1) AS ne,
  ST_Point(@x2, @y2) AS se,
  ST_Point(@x1, @sw_y) AS sw;

-- 4. 生成可用于查询的多边形
SELECT ST_PolygonFromText(CONCAT(
  'POLYGON((',
  @x1,' ',@y1,', ',
  @ne_x,' ',@y1,', ',
  @x2,' ',@y2,', ',
  @x1,' ',@sw_y,', ',
  @x1,' ',@y1,
  '))'
)) AS tile_polygon;

方法2：PHP实现 #

在服务端直接计算，逻辑与SQL一致：

<?php
// 从GET参数获取BBOX：格式为 x1,y1,x2,y2
list($x1, $y1, $x2, $y2) = explode(',', $_GET['bbox']);

// 根据瓦片层级设置UTM分辨率，可通过公式动态计算：
// resolution = 156543.03392 * cos(deg2rad(中心纬度)) / (2^zoom)
$zoom = 18;
$center_lat = 53; // 示例：德国北部纬度
$resolution = 156543.03392 * cos(deg2rad($center_lat)) / pow(2, $zoom);

// 对齐到UTM网格线
$ne_x = floor($x2 / $resolution) * $resolution;
$sw_y = floor($y1 / $resolution) * $resolution;

// 四个顶点坐标
$vertices = [
  [$x1, $y1], // 西北
  [$ne_x, $y1], // 东北
  [$x2, $y2], // 东南
  [$x1, $sw_y], // 西南
  [$x1, $y1] // 闭合多边形
];

// 生成WKT格式多边形
$wkt_parts = array_map(function($point) {
  return sprintf('%f %f', $point[0], $point[1]);
}, $vertices);
$wkt = 'POLYGON((' . implode(', ', $wkt_parts) . '))';
?>

关键注意事项 #

1. 分辨率动态计算：不要硬编码分辨率，根据当前Zoom层级和区域中心纬度用公式计算，确保不同层级、不同区域的网格对齐准确。
2. 空间索引优化：在MariaDB中给多边形字段创建SPATIAL INDEX，可大幅提升基于瓦片区域的查询效率。
3. 边界校验：计算后可额外校验顶点顺序，确保多边形为闭合的有效形状。

内容的提问来源于stack exchange，提问作者LukasKroess