通过 OpenLayers 加载CAD导出位图 和 math.js 构造的仿射变换实现地理坐标系到任意CAD图上像素坐标系的互转

WebGIS开发过程中会遇到这样一种情况:需要使用OpenLayers加载一个未校准的CAD导出的位图;并且还需要通过经纬度坐标数据在这个位图上做一些标记,还需要能通过在OpenLayers取得的图上要素的像素坐标获知实际的经纬度。

总结起来就是两个需求:

加载位图
经纬度坐标与像素坐标互转

分析:

由于从CAD导出的位图并不带有定位信息,所以需要通过仿射变换将图上的像素坐标转换到地理坐标。

即:(左边为某厂的卫星地图,右边为该厂的CAD导出位图,最终实现效果就是用OpenLayers加载位图,并实现坐标转换)
在这里插入图片描述

关于求解仿射变换的过程请见这里。主要的算法思想如下:

Arcgis中就带有了仿射变换的计算模块,OpenLayers没有仿射变换计算的能力,所以使用math.js这个数学库来进行实现。
代码:

算法在上面的截图已经有了,直接用相应的API实现就好:

        //定义仿射变换函数
        function affineTransform(point, from, to) {
            if (from.length != to.length) return;
            //根据参数构造仿射变换所需的矩阵            
            var X = [];
            var Y = [];
            var I = [];
            var U = [];
            var V = [];
            from.forEach((item, index) => {
                X.push(item[0]);
                Y.push(item[1]);
                I.push(1);
                U.push([to[index][0]])
                V.push([to[index][1]])
            })
            //开始最小二乘法的计算过程
            var XYIt = [X, Y, I];
            var resultINV = math.inv(math.multiply(XYIt, math.transpose(XYIt)))
            var resultMulti = math.multiply(resultINV, XYIt);
            var vec1 = math.multiply(resultMulti, U)
            var vec2 = math.multiply(resultMulti, V)
            //使用vec1和vec2计算转换后的坐标
            return [vec1[0][0] * point[0] + point[1] * vec1[1][0] + vec1[2][0], vec2[0][0] * point[0] + point[1] * vec2[1][0] + vec2[2][0]]
        }

CAD导出的位图直接使用ImageStatic加载,并自定义一个像素坐标系:

        //定义地图的像素坐标四至
        var extent = [0, 0, 4000, 2000];
 
        //定义地图的投影坐标系,像素坐标
        var projection = new ol.proj.Projection({
            code: 'factory-image',
            units: 'pixels',
            extent: extent
        });
 
        //初始化地图
        var map = new ol.Map({
            layers: [
                new ol.layer.Image({
                    source: new ol.source.ImageStatic({
                        url: './data/10-9.png',
                        projection: projection,
                        imageExtent: extent
                    })
                })
            ],
            target: 'map',
            view: new ol.View({
                projection: projection,
                center: ol.extent.getCenter(extent),
                zoom: 2,
                maxZoom: 8
            })
        });

概略分别获取图上厂区四角的坐标,图片像素坐标是用potoshop量取的,经纬度坐标是在google地图上拾取的:

        var upperLeft = [119.071450, 39.309006];
        var lowerLeft = [119.074536, 39.305893];
        var upperRight = [119.075858, 39.311641];
        var lowerRight = [119.078934, 39.308527];
 
        var upperLeftPixel = [959, 1897];
        var lowerLeftPixel = [959, 112];
        var upperRightPixel = [2924, 1897];
        var lowerRightPixel = [2924, 112];

包括使用旗杆坐标打点测试的完整代码:

<!DOCTYPE html>
<html>
 
<head>
    <title>厂区地图计算</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v3.20.1/css/ol.css" type="text/css">
    <script src="https://unpkg.com/mathjs@6.2.3/dist/math.js"></script>
    <script src="https://openlayers.org/en/v3.20.1/build/ol.js"></script>
 
</head>
<style>
 
</style>
 
<body>
    <div id="map" class="map"></div>
    <script>
 
        //定义仿射变换函数
        function affineTransform(point, from, to) {
            if (from.length != to.length) return;
            var X = [];
            var Y = [];
            var I = [];
            var U = [];
            var V = [];
            from.forEach((item, index) => {
                X.push(item[0]);
                Y.push(item[1]);
                I.push(1);
                U.push([to[index][0]])
                V.push([to[index][1]])
            })
            var XYIt = [X, Y, I];
            var resultINV = math.inv(math.multiply(XYIt, math.transpose(XYIt)))
            var resultMulti = math.multiply(resultINV, XYIt);
            var vec1 = math.multiply(resultMulti, U)
            var vec2 = math.multiply(resultMulti, V)
            return [vec1[0][0] * point[0] + point[1] * vec1[1][0] + vec1[2][0], vec2[0][0] * point[0] + point[1] * vec2[1][0] + vec2[2][0]]
        }
 
        //Google坐标
 
        var upperLeft = [119.071450, 39.309006];
        var lowerLeft = [119.074536, 39.305893];
        var upperRight = [119.075858, 39.311641];
        var lowerRight = [119.078934, 39.308527];
 
        var upperLeftPixel = [959, 1897];
        var lowerLeftPixel = [959, 112];
        var upperRightPixel = [2924, 1897];
        var lowerRightPixel = [2924, 112];
 
 
 
        //定义地图的像素坐标四至
        var extent = [0, 0, 4000, 2000];
 
        //定义地图的投影坐标系,像素坐标
        var projection = new ol.proj.Projection({
            code: 'factory-image',
            units: 'pixels',
            extent: extent
        });
 
        //初始化地图
        var map = new ol.Map({
            layers: [
                new ol.layer.Image({
                    source: new ol.source.ImageStatic({
                        url: './data/10-9.png',
                        projection: projection,
                        imageExtent: extent
                    })
                })
            ],
            target: 'map',
            view: new ol.View({
                projection: projection,
                center: ol.extent.getCenter(extent),
                zoom: 2,
                maxZoom: 8
            })
        });
 
        //这里用旗杆的坐标演示坐标转换的使用
        var flagPole = [119.077710, 39.309195];
 
        var flagPolePixel = affineTransform(
            flagPole, 
        [upperLeft, lowerLeft, upperRight, lowerRight], 
        [upperLeftPixel, lowerLeftPixel, upperRightPixel, lowerRightPixel]
        )
 
        var p=affineTransform(
            flagPolePixel, [upperLeftPixel, lowerLeftPixel, upperRightPixel, lowerRightPixel],
        [upperLeft, lowerLeft, upperRight, lowerRight]
        
        )
        console.log(flagPole);
        console.log(p);
        var f = new ol.Feature(new ol.geom.Point(flagPolePixel));
 
        var vSource = new ol.source.Vector();
        var vLayer = new ol.layer.Vector({
            source: vSource
        })
 
        vSource.addFeature(f);
 
        f.setStyle(
            new ol.style.Style({
                image: new ol.style.Icon({
                    src: './data/icon.png',
                    anchor: [0.5, 1],
                    scale: 0.3
 
                }),
 
            })
        )
 
        map.addLayer(vLayer);
        map.getView().fit(extent, map.getSize())
        map.render()
 
 
    </script>
</body>
 
</html>

这里使用旗杆的坐标进行了转换和逆转换,并在console里输出,结果如下:
在这里插入图片描述

使用图上像素坐标转换的经纬度,在Google地图上标记旗杆的位置如下:
在这里插入图片描述

经过多次计算和实地对比,精度差距大约在1米以内,比较符合实际需要。