個(gè)人中心
退出
通過 OpenLayers 加載CAD導(dǎo)出位圖 和 math.js 構(gòu)造的仿射變換實(shí)現(xiàn)地理坐標(biāo)系到任意CAD圖上像素坐標(biāo)系的互轉(zhuǎn)
WebGIS開發(fā)過程中會(huì)遇到這樣一種情況：需要使用OpenLayers加載一個(gè)未校準(zhǔn)的CAD導(dǎo)出的位圖；并且還需要通過經(jīng)緯度坐標(biāo)數(shù)據(jù)在這個(gè)位圖上做一些標(biāo)記，還需要能通過在OpenLayers取得的圖上要素的像素坐標(biāo)獲知實(shí)際的經(jīng)緯度。
總結(jié)起來就是兩個(gè)需求：
加載位圖
經(jīng)緯度坐標(biāo)與像素坐標(biāo)互轉(zhuǎn)
分析：
由于從CAD導(dǎo)出的位圖并不帶有定位信息，所以需要通過仿射變換將圖上的像素坐標(biāo)轉(zhuǎn)換到地理坐標(biāo)。
即：（左邊為某廠的衛(wèi)星地圖，右邊為該廠的CAD導(dǎo)出位圖，最終實(shí)現(xiàn)效果就是用OpenLayers加載位圖，并實(shí)現(xiàn)坐標(biāo)轉(zhuǎn)換）
關(guān)于求解仿射變換的過程請(qǐng)見這里。主要的算法思想如下：
Arcgis中就帶有了仿射變換的計(jì)算模塊，OpenLayers沒有仿射變換計(jì)算的能力，所以使用math.js這個(gè)數(shù)學(xué)庫來進(jìn)行實(shí)現(xiàn)。
代碼：
算法在上面的截圖已經(jīng)有了，直接用相應(yīng)的API實(shí)現(xiàn)就好：
//定義仿射變換函數(shù)
function affineTransform(point, from, to) {
if (from.length != to.length) return;
//根據(jù)參數(shù)構(gòu)造仿射變換所需的矩陣
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]])
})
//開始最小二乘法的計(jì)算過程
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計(jì)算轉(zhuǎn)換后的坐標(biāo)
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導(dǎo)出的位圖直接使用ImageStatic加載，并自定義一個(gè)像素坐標(biāo)系：
//定義地圖的像素坐標(biāo)四至
var extent = [0, 0, 4000, 2000];
//定義地圖的投影坐標(biāo)系，像素坐標(biāo)
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
})
});
概略分別獲取圖上廠區(qū)四角的坐標(biāo)，圖片像素坐標(biāo)是用potoshop量取的，經(jīng)緯度坐標(biāo)是在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];
包括使用旗桿坐標(biāo)打點(diǎn)測試的完整代碼：
<!DOCTYPE html>
<html>
<head>
<title>廠區(qū)地圖計(jì)算</title>
<link rel="stylesheet" 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>
//定義仿射變換函數(shù)
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坐標(biāo)
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];
//定義地圖的像素坐標(biāo)四至
var extent = [0, 0, 4000, 2000];
//定義地圖的投影坐標(biāo)系，像素坐標(biāo)
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
})
});
//這里用旗桿的坐標(biāo)演示坐標(biāo)轉(zhuǎn)換的使用
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>
這里使用旗桿的坐標(biāo)進(jìn)行了轉(zhuǎn)換和逆轉(zhuǎn)換，并在console里輸出，結(jié)果如下：
使用圖上像素坐標(biāo)轉(zhuǎn)換的經(jīng)緯度，在Google地圖上標(biāo)記旗桿的位置如下：
經(jīng)過多次計(jì)算和實(shí)地對(duì)比，精度差距大約在1米以內(nèi)，比較符合實(shí)際需要。
分類導(dǎo)航