最近看到 deck.gl 在 medium 上的一篇文章，介绍了解决高缩放等级下「抖动」问题的思路。
How (sometimes) assuming the Earth is “flat” helps speed up rendering in deck.gl
我写了一个简单的 DEMO 演示这种现象：
GLSL 浮点数精度
出现抖动的原因是什么呢？
deck.gl 不同于 Mapbox，墨卡托坐标转换是在 vertex shader 中完成的。在 shader 中进行墨卡托投影，我们在本专栏上一篇「绘制地图」中已经介绍过了：
// deck.gl/shaders/project.glslvec2 project_mercator_(vec2 lnglat) {return vec2(radians(lnglat.x) + PI,PI - log(tan(PI * 0.25 + radians(lnglat.y) * 0.5)));}
因此我们需要将经纬度坐标传入 shader 中，然而 GLSL 中使用的是 32-bit float，不同于 JS 默认的 64-bit 双精度浮点数。在转换过程中就会丢失精度，为此 JS 提供了 Math.fround 方法从 64-bit 转换到 32-bit ：
Math.fround(-122.4000588);// -122.40006256103516
但是在转换过程中就会出现明显的精度丢失，在高缩放等级下就会造成肉眼可见的明显偏移，即「抖动」现象。并且随着缩放等级提升，误差将越来越大。下图来自 deck.gl medium 原文。
为了解决这个问题，deck.gl 早期在 GLSL 中使用两个 float 拼接成 64-bit 解决，但是对于 shader 编译和解析性能都有影响，需要在 JS 和 GLSL 中频繁 encode/decode，同时也会增加 CPU 向 GPU 传递的数据带宽。
而 v6.2 版本后 deck.gl使用了一种以屏幕中心作为动态坐标原点的「Offset Coords」方案，仅仅使用 32-bit 就解决了这个问题。
偏移坐标系
很自然想到，相近的两个坐标之差正好可以将高位抹去，只需要使用 32-bit 来存储差值，精度就完全足够了。
因此我们需要选取一个固定点，用来与其他顶点取差值。那么如何选择这个固定点呢？使用视口中心点是个很好的选择，而计算固定点和每个顶点坐标之差应该在 shader 中完成。那么能不能在 CPU 中算好差值传入呢？GPU 对于这种简单运算效率较高，而如果在 CPU 中，每一帧视口中心点的经纬度坐标都（可能）在改变，需要将原始数据都计算一遍，显然不行：
// deck.gl/shaders/project.glslvec4 project_position(vec4 position) {// 处理经纬度 offsetif (project_uCoordinateSystem == COORDINATE_SYSTEM_LNGLAT_AUTO_OFFSET) {// 与视口中心点的偏移，在经纬度坐标下因此可以保留低位，32-bit float 足够用float X = position.x - project_coordinate_origin.x;float Y = position.y - project_coordinate_origin.y;// 后续着重介绍return project_offset_(vec4(X, Y, position.z, position.w));} else {// 省略正常处理经纬度 -> 世界坐标return vec4(project_mercator(position.xy) * WORLD_SCALE * u_project_scale,project_scale(position.z),position.w);}}
因此 deck.gl 依据当前缩放等级，在正常和 Offset 两种坐标系之间切换，如果缩放等级大于 12，则需要计算出视口中心在经纬度坐标系下的坐标：
const LNGLAT_AUTO_OFFSET_ZOOM_THRESHOLD = 12;if (coordinateZoom < LNGLAT_AUTO_OFFSET_ZOOM_THRESHOLD) {} else {// 使用 Offset 坐标，传入经纬度坐标系下的视口中心点const lng = Math.fround(viewport.longitude);const lat = Math.fround(viewport.latitude);shaderCoordinateOrigin = [lng, lat];}
回到 vertex shader 中，最终像素空间的计算结果可以拆解成两部分：世界坐标系偏移进行正常矩阵运算以及视口中心点的最终投影结果。
// 处理 offset 和正常经纬度到世界坐标系转换vec4 project_pos = project_position(vec4(a_pos, 0.0, 1.0));gl_Position = u_mvp_matrix * project_pos + u_viewport_center_projection;
视口中心点最终的投影结果可以在 CPU 中每一个渲染帧计算，例如在 mapbox custom layer 中可以使用相机参数创建视图和投影矩阵。值得一提的是 uber 也开源了 viewport-mercator-project ，能够很方便地完成视图和投影矩阵的计算：
// 从 mapbox 上下文中获取相机参数const center = this.map.getCenter();const currentZoomLevel = this.map.getZoom();const pitch = this.map.getPitch();const center = this.map.getCenter();// /uber-common/viewport-mercator-projectconst viewport = new WebMercatorViewport({width: gl.drawingBufferWidth / 2,height: gl.drawingBufferHeight / 2,longitude: center.lng,latitude: center.lat,zoom: currentZoomLevel,pitch,bearing,});// 当前缩放等级下墨卡托坐标 -> 世界坐标，在 JS 中使用 64-bit 浮点数计算const positionPixels = viewport.projectFlat([ Math.fround(center.lng), Math.fround(center.lat) ],Math.pow(2, currentZoomLevel));// 视图投影后的最终结果const projectionCenter = vec4.transformMat4([],[positionPixels[0], positionPixels[1], 0.0, 1.0],viewProjectionMatrix);// 将视口中心点的最终投影结果传入 uniformdrawParams['u_viewport_center_projection'] = projectionCenter;
现在只剩下了最后一部分，我们已经得到了偏移坐标系下精确的经纬度差，如何使用这个差值来估计世界坐标系下的偏移量呢？
泰勒级数展开
为了使用墨卡托坐标系下的差值估计世界坐标系下的差值，可以使用泰勒级数展开。关于泰勒级数，可以阅读「怎样更好地理解并记忆泰勒展开式？」
有了
处的值，结合导数就可以对 在任意点进行估计。在偏移坐标系场景下， 就是动态的屏幕中心点的坐标，每个顶点与中心点的差就是 ，而 就是偏移坐标系到世界坐标系的转换公式。
由于 X 轴方向是线性的，对于 x 就可以应用线性估计：
而 Y 轴方向的变化是非线性的，可以使用 2 级泰勒展开：
对于经纬度坐标下的偏移转换到世界坐标，其中 ：
在 GLSL 中使用向量运算可以快速实现上述转换公式：其中 u_pixels_per_degree 对应 [K11, K21] ，而 u_pixels_per_degree2 对应 [0, K22]：
// offset:[delta lng, delta lat]vec4 project_offset(vec4 offset) {float dy = offset.y;dy = clamp(dy, -1., 1.);vec3 pixels_per_unit = u_pixels_per_degree + u_pixels_per_degree2 * dy;return vec4(offset.xyz * pixels_per_unit, offset.w);}
对于u_pixels_per_degree 和 u_pixels_per_degree2 的计算过程如下 。注释十分详细的介绍了求导过程：
export function getDistanceScales({latitude, longitude, zoom, scale, highPrecision = false}) {// Calculate scale from zoom if not providedscale = scale !== undefined ? scale : zoomToScale(zoom);const result = {};const worldSize = TILE_SIZE * scale;const latCosine = Math.cos(latitude * DEGREES_TO_RADIANS);/*** Number of pixels occupied by one degree longitude around current lat/lon:pixelsPerDegreeX = d(lngLatToWorld([lng, lat])[0])/d(lng)= scale * TILE_SIZE * DEGREES_TO_RADIANS / (2 * PI)pixelsPerDegreeY = d(lngLatToWorld([lng, lat])[1])/d(lat)= -scale * TILE_SIZE * DEGREES_TO_RADIANS / cos(lat * DEGREES_TO_RADIANS) / (2 * PI)*/const pixelsPerDegreeX = worldSize / 360;const pixelsPerDegreeY = pixelsPerDegreeX / latCosine;/*** LngLat: longitude -> east and latitude -> north (bottom left)* UTM meter offset: x -> east and y -> north (bottom left)* World space: x -> east and y -> south (top left)** Y needs to be flipped when converting delta degree/meter to delta pixels*/result.pixelsPerDegree = [pixelsPerDegreeX, -pixelsPerDegreeY, altPixelsPerMeter];/*** Taylor series 2nd order for 1/latCosinef'(a) * (x - a)= d(1/cos(lat * DEGREES_TO_RADIANS))/d(lat) * dLat= DEGREES_TO_RADIANS * tan(lat * DEGREES_TO_RADIANS) / cos(lat * DEGREES_TO_RADIANS) * dLat*/if (highPrecision) {const latCosine2 = DEGREES_TO_RADIANS * Math.tan(latitude * DEGREES_TO_RADIANS) / latCosine;const pixelsPerDegreeY2 = pixelsPerDegreeX * latCosine2 / 2;result.pixelsPerDegree2 = [0, -pixelsPerDegreeY2, altPixelsPerDegree2];}// Main results, used for converting meters to latlng deltas and scaling offsetsreturn result;}
至此对于偏移坐标系下的处理就完成了，高缩放等级下也不会抖动了，效果如下DEMO：
总结
Mapbox 是在 CPU 中进行墨卡托投影变换，因此并不会出现这个问题，可以查看官方的 DEMO。个人认为 Mapbox 对于数据进行了分片裁剪，有效控制了在 CPU 中进行的运算量，因此可以这么做，而 deck.gl 每次都需要渲染全量数据，所以将需要一次 log 和 tan 的投影运算放到 GPU 中就很有必要了。
另外关于 deck.gl，除了官方文档，知乎上 @hijiangtao 的系列文章 hijiangtao：从零开始学习时空数据可视化（序）也很值得关注。
本文使用到 DEMO 地址：
xiaoiver/custom-mapbox-layer
参考资料
「How (sometimes) assuming the Earth is “flat” helps speed up rendering in deck.gl」
「RFC: Improved 32-bit LNGLAT projection mode」
「Introduction to Taylor's theorem for multivariable functions」
