原创 freegiser Spatial Data 2021-03-18 17:30
收录于合集
#MapboxGL5
#气象7
#WebGIS19
#GIS17
一 背景
WebGIS中非常常见的一个功能，就是叠加地理相关的Image到底图上，众所周知，底图是有地理坐标系的，Image通常也是基于某个地理坐标系制图的产物，如果底图坐标系和Image坐标系不一致，那么他们叠加肯定就会发生位置的变形与偏移：
不同投影叠加形变
上图，地理底图一般是墨卡托投影(epsg:3857)，假设叠加的图片是生产单位按照WGS84经纬度(epsg:4326)制图输出的，并且已知这个图片真实地理范围是[84minx,84miny,84maxx,84maxy]，常规思路就是计算84的范围对应底图的墨卡托的范围[mkt_minx,mkt_miny,mkt_maxx,mkt_maxy]，然后工工整整把这个图像贴到底图这个范围区间里，效果如上图，非常明显的投影形变。
看到这很多朋友会问：有方法让与底图坐标系不同的Image正确叠加底图吗？答案只有一个：让image和底图的坐标系相同。
那么下个问题：如何能让image和底图坐标系相同？
一个是逐像素转换，这不是本篇话题，计算量比较大，性能不好，通常不被采用。
另一个方法是：将图像分解三角形进行投影转换（很熟悉的味道啊，图像？三角形？对应webgl的纹理和三角片？）。
image投影到与底图一致
二 技术资料
openlayers目前原生就已经支持image重投影了，demo如下：
https:///en/latest/examples/reprojection-image.html
具体实现技术原理参考如下博客：
https://www.jianshu.com/p/98da62d88e83Openlayers raster reprojection 栅格重投影原理分析
只不过ol目前的源码里是用canvas实现的，在leaflet的里也有大神提供了一个webgl实现的图层，网速慢要等图像加载完毕：
https://ivansanchez./Leaflet.ImageOverlay.Arrugator/demo.html
本质上ol和leaflet思路都是图像拆分三角形，三角形整体投影（相比逐像素投影效率高非常多），很像数学里的积分概念，拆分三角形越多，那么计算越精确，可以想象，三角形足够的的多，是不是等于非常精确，和逐像素一样了？
下面我们主要基于leaflet的这个webgl方案做下简单说明，开发者是基于arrugator库实现三角形切分的：
https:///IvanSanchez/arrugator
这个库使用很简单：
第一步：定义图像坐标系
proj4.defs( "EPSG:25833", "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs");// 定义投影转换规则，如从图像的25833转底图的3857坐标系const projector = proj4("EPSG:25833", "EPSG:3857").forward;第二步：定义图像的控制坐标点（图像对应的地理bbox）与对应纹理的uv坐标
//图像的地理bbox坐标const controlPoints = [ [ -183622.3, 7996344.0 ], // top-left [ -183622.3, 6396344.0 ], // bottom-left [ 1416377.7, 7996344.0 ], // top-right [ 1416377.7, 6396344.0 ], // bottom-right];//纹理uv坐标，注意top left这些方向与地理的一一对应const sourceUV = [ [ 0, 0 ], // top-left [ 0, 1 ], // bottom-left [ 1, 0 ], // top-right [ 1, 1 ], // bottom-right];
第三步：计算三角形
let arruga = new Arrugator( forward, controlPoints, sourceUV, [ [ 0, 1, 3 ], [ 0, 3, 2 ] ] // 这里写死不变);// 根据阈值切分三角形const epsilon = 0.0000000001; // 不同坐标系这个阈值是不同的，实际使用可以自行调试arruga.lowerEpsilon(epsilon);// 结果输出let arrugado = arruga.output();
这里有个问题，阈值设置多少合适？开发者可以根据实际的坐标系转换时候，调试以下代码，先保证这个条件成立，然后在这个条件成立的基础上，epsilon越小，三角形切分越多，那么耗时多点，结果正确点呗。
例如：peek().epsilon是0.1，那么你传阈值必须要小于0.1，大于的话就不切分了其实，那我可以设置0.01，0.001都可以的，自己调试效果和性能，取个折衷数值较好，就是不能设置太小追求精确度，但牺牲性能，中庸点好。
这个库最终输出的其实有三个东西：
拿到地理顶点，uv顶点，indices的绘制索引，再加一张image作为纹理，使用webgl渲染简直是水到渠成的事。
三 mapboxgl实践
由于mapboxgl本身不支持image的投影转换，我们使用arrugator来帮我们实现这个功能，选择openlayers这个例子来改写到mapbox上。
https:///en/latest/examples/reprojection-image.html
首先对Image三角形切分：
import test_map from '../../images/2000px-British_National_Grid.svg.png';import proj4 from 'proj4';import Arrugator from 'arrugator';
// 墨卡托投影的左上角坐标，对应mapbox左上角起始坐标[0,0]const origin=[ -20037508.342789244,20037508.342789244];proj4.defs( 'EPSG:27700', '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 ' + '+x_0=400000 +y_0=-100000 +ellps=airy ' + '+towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 ' + '+units=m +no_defs' );// 定义转换规则const projector = proj4("EPSG:27700", "EPSG:3857").forward;// 改写坐标转换函数// 原因是mapbox的墨卡托坐标是0-1区间，并且对应地理范围与标准3857不同。function forward(coors) { // 墨卡托坐标 const coor_3857 = projector(coors); // 墨卡托坐标转换到 0-1区间，origin对应mapbox 0 0点，所以这么写 const mapbox_coor1 = Math.abs((coor_3857[0]-origin[0])/(20037508.342789244*2)); const mapbox_coor2 = Math.abs((coor_3857[1]-origin[1])/(20037508.342789244*2)); return [mapbox_coor1,mapbox_coor2];}const epsilon = 0.00000000001;const controlPoints = [ [ 0, 1300000 ], // top-left [ 0, 0 ], // bottom-left [ 700000, 1300000 ], // top-right [ 700000, 0 ], // bottom-right];//纹理uv坐标const sourceUV = [ [ 0, 0 ], // top-left [ 0, 1 ], // bottom-left [ 1, 0 ], // top-right [ 1, 1 ], // bottom-right];let arruga = new Arrugator( forward, controlPoints, sourceUV, [ [ 0, 1, 3 ], [ 0, 3, 2 ] ] );arruga.lowerEpsilon(epsilon);// arrugado里存了顶点，索引 let arrugado = arruga.output();
然后构建一个自定义图层，部分webgl自封装的方法，就不展开说了：
class ImageLayer { // 传入图像和切分三角形顶点和索引 constructor(img, arrugado) { this.id = 'test_layer'; this.type = 'custom'; this.renderingMode = '2d'; this.img = img; this.arrugado = arrugado; this.map; this.gl; } onAdd(m, gl) { this.map = m; this.gl = gl; // 墨卡托坐标转mapbox的墨卡托坐标 let pos = this.arrugado.projected.flat(); // uv纹理 let uv = this.arrugado.uv.flat(); // 三角形index let trigs = arrugado.trigs.flat(); // 着色器代码 const vs = `#version 300 es layout(location=0) in vec2 a_position; layout(location=1) in vec2 a_uv; uniform mat4 u_worldViewProjection; out vec2 v_uv; void main() { gl_Position = u_worldViewProjection*vec4(a_position, 0.0,1.0); v_uv = a_uv; }`; const fs = `#version 300 es precision highp int; precision highp float; uniform sampler2D u_image; in vec2 v_uv; out vec4 outColor; void main() { // 根据纹理坐标获取颜色 vec4 texel= texture(u_image, v_uv); if (texel.r > 0.95 && texel.g > 0.95 && texel.b > 0.95) { texel = vec4(0.6, 0.0, 0.0, 1.0); } else { outColor= vec4(texel.r*0.4,texel.g*0.4,texel.b*0.4,0.4); } }`; // 定义渲染model this.drawModel = createModel(gl, vs, fs); //定义vao对象，绑定顶点，索引 this.vao = gl.createVertexArray(); gl.bindVertexArray(this.vao); const positionBuffer = createBuffer(gl, gl.ARRAY_BUFFER, new Float32Array(pos)); bindAttribute(gl, positionBuffer, 0, 2); const uvBuffer = createBuffer(gl, gl.ARRAY_BUFFER, new Float32Array(uv)); bindAttribute(gl, uvBuffer, 1, 2); const indexBuffer = createBuffer(gl, gl.ELEMENT_ARRAY_BUFFER, new Uint16Array(trigs)); this.position_count = trigs.length;
// 绑定结束 gl.bindBuffer(gl.ARRAY_BUFFER, null); gl.bindVertexArray(null); //构造纹理并绑定 const texture = createTexture2D(gl, { data: this.img, mipLevel: 0, internalFormat: gl.RGBA,//webgl中格式 srcFormat: gl.RGBA,//输入数据源格式 type: gl.UNSIGNED_BYTE, height: 2000, width: 3696, parameters: { [ gl.TEXTURE_MAG_FILTER ]: gl.LINEAR, [ gl.TEXTURE_MIN_FILTER ]: gl.LINEAR, [ gl.TEXTURE_WRAP_S ]: gl.CLAMP_TO_EDGE, [ gl.TEXTURE_WRAP_T ]: gl.CLAMP_TO_EDGE } }); //drawProgram绑定纹理，不能绑定单元0 gl.useProgram(this.drawModel.program); bindTexture2D(gl, texture, 3); gl.uniform1i(this.drawModel.u_image, 3); } render(gl, matrix) { gl.useProgram(this.drawModel.program); //设置unifrom gl.uniformMatrix4fv(this.drawModel.u_worldViewProjection, false, matrix); //绑定顶点vao gl.bindVertexArray(this.vao); gl.drawElements(gl.TRIANGLES, this.position_count, gl.UNSIGNED_SHORT, 0); //如果取消绑定，会报错GL_INVALID_OPERATION: Insufficient buffer size. gl.bindVertexArray(null); } onRemove(map, gl) { // 。。。 }}
第三步：新建图层，叠加地图
let img = new Image(); img.src = test_map; img.onload = function () { let layer = new ImageLayer(img, arrugado); var container = document.createElement('div'); const width = 1600, height = 900; container.style.width = `${width}px`; container.style.height = `${height}px`; document.body.appendChild(container); // 兼容webgl2 mapboxglExtension(mapboxgl); var map = new mapboxgl.Map({ container: container, style: 'mapbox://styles/mapbox/dark-v10', center: [ 0, 65 ], zoom: 4, antialias: true }); map.on('load', () => { map.addLayer(layer); }); }至于怎么样让mapboxgl支持webgl2，为什么纹理单元不能写0，为什么这么多手动绑定解绑webgl资源对象，参考下文《使用WebGL2进行Mapbox自定义图层开发》：
https://mp.weixin.qq.com/s?__biz=Mzg2OTUxMzM2MA==&mid=2247483684&idx=1&sn=cbec2c833fa0a2a30e3ee0d6063fbf0c&chksm=ce9aa0dbf9ed29cd1b65bebf5e773eb8b4005c3d347b033426d0c4b65807ab74934268d9df88&token=1580692165&lang=zh_CN#rd
效果如下图：
四 应用延申
这个image重投影一般在哪用咧？比如气象领域，很多产品出的栅格图像，想要叠加到地图上，坐标系可能图像用4326经纬度，底图是墨卡托，这时候叠加歪了，就可以使用这个库去转换，例如下图：
转换前
转换后
转换前，日本列岛的温度都叠加底图上的日本区域，都飘上面去了，越往北偏移越大其实，赤道那近似，这个就是投影变换的原理了，球面转平面的关系，这里讲不清楚，自行研究。转换后，可以看到叠加的基本严丝合缝了。
