目錄
作為一個GISer，在日常WebGIS開發中，會常用到的turf.js，這是一個地理空間分析的JavaScript庫，經常搭配各種GIS JS API使用，如leaflet、mapboxgl、openlayers等；在後臺Java開發中，也有個比較強大的GIS庫，geotools，裡面包含構建一個完整的地理資訊系統所需要的全部工具類；資料庫端常用是postgis擴充套件，需要在postgres庫中引入使用。
然而在開發某一些業務系統的時候，有些需求只需要呼叫某一個GIS演算法，簡單的幾行程式碼即可完成，沒有必要去引用一個GIS類庫。
而且有些演算法在這些常用的GIS類庫中沒有對應介面，就比如在下文記錄的這幾種常用演算法中，求垂足、判斷線和麵的關係，在turf.js就沒有對應介面。
下面文章中是我總結的一些常用GIS演算法，這裡統一用JavaScript語言實現，因為JS程式碼相對比較簡潔，方便理解其中演算法邏輯，也方便在瀏覽器下預覽效果。在具體應用時可以根據具體需求，翻譯成Java、C#、Python等語言來使用。
文中程式碼大部分為之前遇到需求時在網上搜尋得到，然後自己根據具體需要做了優化修改，通過這篇文章做個總結收集，也方便後續使用時查詢。
1、常用演算法
以下方法中傳參的點、線、面都是對應geojson格式中coordinates，方便統一呼叫。geojson標準參考：https://www.oschina.net/translate/geojson-spec
1.1、計算兩經緯度點之間的距離
適用場景：測量
/**
* 計算兩經緯度點之間的距離(單位：米)
* @param p1 起點的座標；[經度,緯度]；例：[116.35,40.08]
* @param p2 終點的座標；[經度,緯度]；例：[116.72,40.18]
*
* @return d 返回距離
*/
function getDistance(p1, p2) {
var rlat1 = p1[1] * Math.PI / 180.0;
var rlat2 = p2[1] * Math.PI / 180.0;
var a = rlat1 - rlat2;
var b = p1[0] * Math.PI / 180.0 - p2[0] * Math.PI / 180.0;
var d = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2) + Math.cos(rlat1) * Math.cos(rlat2) * Math.pow(Math.sin(b / 2), 2)));
d = d * 6378.137;
d = Math.round(d * 10000) / 10;
return d
}
1.2、根據已知線段以及到起點距離，求目標點座標
適用場景：封閉管段定位問題點
/**
* 根據已知線段以及到起點距離（單位：米），求目標點座標
* @param line 線段；[[經度,緯度],[經度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param dis 到起點距離（米）；Number；例：500
*
* @return point 返回座標
*/
function getLinePoint(line, dis) {
var p1 = line[0]
var p2 = line[1]
var d = getDistance(p1, p2) // 計算兩經緯度點之間的距離(單位：米)
var dx = p2[0] - p1[0]
var dy = p2[1] - p1[1]
return [p1[0] + dx * (dis / d), p1[1] + dy * (dis / d)]
}
1.3、已知點、線段，求垂足
垂足可能線上段上，也可能線上段延長線上。
適用場景：求垂足
/**
* 已知點、線段，求垂足
* @param line 線段；[[經度,緯度],[經度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param p 點；[經度,緯度]；例：[116.35,40.08]
*
* @return point 返回垂足座標
*/
function getFootPoint(line, p) {
var p1 = line[0]
var p2 = line[1]
var dx = p2[0] - p1[0];
var dy = p2[1] - p1[1];
var cross = dx * (p[0] - p1[0]) + dy * (p[1] - p1[1])
var d2 = dx * dx + dy * dy
var u = cross / d2
return [(p1[0] + u * dx), (p1[1] + u * dy)]
}
1.4、線段上距離目標點最近的點
不同於上面求垂足方法，該方法求出的點肯定線上段上。
如果垂足線上段上，則最近的點就是垂足，如果垂足線上段延長線上，則最近的點就是線段某一個端點。
適用場景：根據求出最近的點計算點到線段的最短距離
/**
* 線段上距離目標點最近的點
* @param line 線段；[[經度,緯度],[經度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param p 點；[經度,緯度]；例：[116.35,40.08]
*
* @return point 最近的點座標
*/
function getShortestPointInLine(line, p) {
var p1 = line[0]
var p2 = line[1]
var dx = p2[0] - p1[0];
var dy = p2[1] - p1[1];
var cross = dx * (p[0] - p1[0]) + dy * (p[1] - p1[1])
if (cross <= 0) {
return p1
}
var d2 = dx * dx + dy * dy
if (cross >= d2) {
return p2
}
// 垂足
var u = cross / d2
return [(p1[0] + u * dx), (p1[1] + u * dy)]
}
1.5、點緩衝
這裡緩衝屬於測地線方法，由於這裡並沒有嚴格的投影轉換體系，所以與標準的測地線緩衝還有些許誤差，不過經測試，半徑100KM內，誤差基本可以忽略。具體緩衝型別可看下之前的文章你真的會用PostGIS中的buffer緩衝嗎？
適用場景：根據點和半徑畫圓
/**
* 點緩衝
* @param center 中心點；[經度,緯度]；例：[116.35,40.08]
* @param radius 半徑（米）；Number；例：5000
* @param vertices 返回圓麵點的個數；預設64；Number；例：32
*
* @return coords 面的座標
*/
function bufferPoint(center, radius, vertices) {
if (!vertices) vertices = 64;
var coords = []
// 111319.55：在赤道上1經度差對應的距離，111133.33：在經線上1緯度差對應的距離
var distanceX = radius / (111319.55 * Math.cos(center[1] * Math.PI / 180));
var distanceY = radius / 111133.33;
var theta, x, y;
for (var i = 0; i < vertices; i++) {
theta = (i / vertices) * (2 * Math.PI);
x = distanceX * Math.cos(theta);
y = distanceY * Math.sin(theta);
coords.push([center[0] + x, center[1] + y]);
}
return [coords]
}
1.6、點和麵關係
該方法採用射線法思路實現。（瞭解射線法可參考：https://blog.csdn.net/qq_27161673/article/details/52973866）
這裡已經考慮到環狀多邊形的情況。
適用場景：判斷點是否在面內
/**
* 點和麵關係
* @param point 點；[經度,緯度]；例：[116.353455, 40.080173]
* @param polygon 面；geojson格式中的coordinates；例：[[[116.1,39.5],[116.1,40.5],[116.9,40.5],[116.9,39.5]],[[116.3,39.7],[116.3,40.3],[116.7,40.3],[116.7,39.7]]]
*
* @return inside 點和麵關係；0:多邊形外，1：多邊形內，2：多邊形邊上
*/
function pointInPolygon(point, polygon) {
var isInNum = 0;
for (var i = 0; i < polygon.length; i++) {
var inside = pointInRing(point, polygon[i])
if (inside === 2) {
return 2;
} else if (inside === 1) {
isInNum++;
}
}
if (isInNum % 2 == 0) {
return 0;
} else if (isInNum % 2 == 1) {
return 1;
}
}
/**
* 點和麵關係
* @param point 點
* @param ring 單個閉合面的座標
*
* @return inside 點和麵關係；0:多邊形外，1：多邊形內，2：多邊形邊上
*/
function pointInRing(point, ring) {
var inside = false,
x = point[0],
y = point[1],
intersects, i, j;
for (i = 0, j = ring.length - 1; i < ring.length; j = i++) {
var xi = ring[i][0],
yi = ring[i][1],
xj = ring[j][0],
yj = ring[j][1];
if (xi == xj && yi == yj) {
continue
}
// 判斷點與線段的相對位置，0為線上段上，>0 點在左側，<0 點在右側
if (isLeft(point, [ring[i], ring[j]]) === 0) {
return 2; // 點在多邊形邊上
} else {
if ((yi > y) !== (yj > y)) { // 垂直方向目標點在yi、yj之間
// 求目標點在當前線段上的x座標。 由於JS小數運算後會轉換為精確15位的float，因此需要去一下精度
var xx = Number(((xj - xi) * (y - yi) / (yj - yi) + xi).toFixed(10))
if (x <= xx) { // 目標點水平射線與當前線段有交點
inside = !inside;
}
}
}
}
return Number(inside);
}
/**
* 判斷點與線段的相對位置
* @param point 目標點
* @param line 線段
*
* @return isLeft，點與線段的相對位置，0為線上段上，>0 p在左側，<0 p在右側
*/
function isLeft(point, line) {
var isLeft = ((line[0][0] - point[0]) * (line[1][1] - point[1]) - (line[1][0] - point[0]) * (line[0][1] - point[1]))
// 由於JS小數運算後會轉換為精確15位的float，因此需要去一下精度
return Number(isLeft.toFixed(10))
}
1.7、線段與線段的關係
適用場景：判斷線和線的關係
/**
* 線段與線段的關係
* @param line1 線段；[[經度,緯度],[經度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param line2 線段；[[經度,緯度],[經度,緯度]]；例：[[116.33,40.21],[116.36,39.76]]
*
* @return intersect 線段與線段的關係；0:相離，1：相交，2：相切
*/
function intersectLineAndLine(line1, line2) {
var x1 = line1[0][0],
y1 = line1[0][1],
x2 = line1[1][0],
y2 = line1[1][1],
x3 = line2[0][0],
y3 = line2[0][1],
x4 = line2[1][0],
y4 = line2[1][1]
//快速排斥：
//兩個線段為對角線組成的矩形，如果這兩個矩形沒有重疊的部分，那麼兩條線段是不可能出現重疊的
//這裡的確如此，這一步是判定兩矩形是否相交
//1.線段ab的低點低於cd的最高點（可能重合）
//2.cd的最左端小於ab的最右端（可能重合）
//3.cd的最低點低於ab的最高點（加上條件1，兩線段在豎直方向上重合）
//4.ab的最左端小於cd的最右端（加上條件2，兩直線在水平方向上重合）
//綜上4個條件，兩條線段組成的矩形是重合的
//特別要注意一個矩形含於另一個矩形之內的情況
if (!(Math.min(x1, x2) <= Math.max(x3, x4) && Math.min(y3, y4) <= Math.max(y1, y2) &&
Math.min(x3, x4) <= Math.max(x1, x2) && Math.min(y1, y2) <= Math.max(y3, y4))) {
return 0
}
// 判斷點與線段的相對位置，0為線上段上，>0 點在左側，<0 點在右側
if (isLeft(line1[0], line2) === 0 || isLeft(line1[1], line2) === 0) {
return 2
}
//跨立實驗：
//如果兩條線段相交，那麼必須跨立，就是以一條線段為標準，另一條線段的兩端點一定在這條線段的兩段
//也就是說a b兩點線上段cd的兩端，c d兩點線上段ab的兩端
var kuaili1 = ((x3 - x1) * (y2 - y1) - (x2 - x1) * (y3 - y1)) * ((x4 - x1) * (y2 - y1) - (x2 - x1) * (y4 - y1))
var kuaili2 = ((x1 - x3) * (y4 - y3) - (x4 - x3) * (y1 - y3)) * ((x2 - x3) * (y4 - y3) - (x4 - x3) * (y2 - y3))
return Number(Number(kuaili1.toFixed(10)) <= 0 && Number(kuaili2.toFixed(10)) <= 0)
}
1.8、線和麵關係
適用場景：判斷線與面的關係
該方法考慮到環狀多邊形的情況，且把相切情況分為了內切和外切。
參考連結：https://www.cnblogs.com/xiaozhi_5638/p/4165353.html
/**
* 線和麵關係
* @param line 線段；[[經度,緯度],[經度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param polygon 面；geojson格式中的coordinates；例：[[[116.1,39.5],[116.1,40.5],[116.9,40.5],[116.9,39.5]],[[116.3,39.7],[116.3,40.3],[116.7,40.3],[116.7,39.7]]]
*
* @return intersect 線和麵關係；0:相離，1：相交，2：包含，3：內切，4：外切
*/
function intersectLineAndPolygon(line, polygon) {
var isTangent = false
var isInNum = 0
var intersect = 0
for (var i = 0; i < polygon.length; i++) {
// 線和麵關係；0:相離，1：相交，2：包含，3：內切，4：外切
intersect = intersectLineAndRing(line, polygon[i])
if (intersect === 1) {
return 1
} else if (intersect === 2) {
isInNum++
} else if (intersect === 3) {
isInNum++
isTangent = true
} else if (intersect === 4) {
isTangent = true
}
}
if (isInNum % 2 == 0) {
if (isTangent) {
return 4 // 外切
} else {
return 0 // 相離
}
} else if (isInNum % 2 == 1) {
if (isTangent) {
return 3 // 內切
} else {
return 2 // 包含
}
}
}
/**
* 線和麵關係
* @param line 線段
* @param ring 單面
*
* @return intersect 線和麵關係；0:相離，1：相交，2：包含，3：內切，4：外切
*/
function intersectLineAndRing(line, ring) {
var inserset = 0
var isTangent = false
var inserset1 = pointInRing(line[0], ring) // 點和麵關係；0:多邊形外，1：多邊形內，2：多邊形邊上
var inserset2 = pointInRing(line[1], ring) // 點和麵關係；0:多邊形外，1：多邊形內，2：多邊形邊上
if (inserset1 === inserset2 === 0) {
inserset = 0
} else if ((inserset1 * inserset2) === 1) {
inserset = 2
} else if ((inserset1 * inserset2) === 2) {
inserset = 3
} else if ((inserset1 === 2 || inserset2 === 2) && (inserset1 === 0 || inserset2 === 0)) {
inserset = 4
} else if ((inserset1 === 1 || inserset2 === 1) && (inserset1 === 0 || inserset2 === 0)) {
return 1 // 相交
}
for (var i = 0, j = ring.length - 1; i < ring.length; j = i++) {
var line2 = [ring[j], ring[i]]
// 目標線段與當前線段的關係；0:相離，1：相交，2：相切
var intersectLine = intersectLineAndLine(line, line2)
if (intersectLine == 1) {
return 1 // 相交
}
}
return inserset
}
1.9、geojson 面轉線
適用場景：只有geojson面資料，獲取線的邊界
/**
* 面轉線
* @param geojson 面geojson
*
* @return geojson 線geojson
*/
function convertPolygonToPolyline(polygonGeoJson) {
var polylineGeoJson = JSON.parse(JSON.stringify(polygonGeoJson))
for (var i = 0; i < polylineGeoJson.features.length; i++) {
var MultiLineString = []
if (polylineGeoJson.features[i].geometry.type === 'Polygon') {
var Polygon = polylineGeoJson.features[i].geometry.coordinates
Polygon.forEach(LinearRing => {
var LineString = LinearRing
MultiLineString.push(LineString)
})
} else if (polylineGeoJson.features[i].geometry.type === 'MultiPolygon') {
var MultiPolygon = polylineGeoJson.features[i].geometry.coordinates
MultiPolygon.forEach(Polygon => {
Polygon.forEach(LinearRing => {
var LineString = LinearRing
MultiLineString.push(LineString)
})
})
} else {
console.error('請確認輸入引數為geojson格式面資料！')
return null
}
polylineGeoJson.features[i].geometry.type = 'MultiLineString' //面轉線
polylineGeoJson.features[i].geometry.coordinates = MultiLineString
}
return polylineGeoJson
}
2、線上示例
線上示例：http://gisarmory.xyz/blog/index.html?demo=GISAlgorithm
程式碼地址：http://gisarmory.xyz/blog/index.html?source=GISAlgorithm
原文地址：http://gisarmory.xyz/blog/index.html?blog=GISAlgorithm
關注《GIS兵器庫》， 只給你網上搜不到的GIS知識技能。
本文章採用 知識共享署名-非商業性使用-相同方式共享 4.0 國際許可協議 進行許可。歡迎轉載、使用、重新發布，但務必保留文章署名《GIS兵器庫》（包含連結： http://gisarmory.xyz/blog/），不得用於商業目的，基於本文修改後的作品務必以相同的許可釋出。