作为一个 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/trans…
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_2716…）
这里曾经思考到环状多边形的状况。
实用场景：判断点是否在面内
/**
* 点和面关系
* @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/xiaoz…
/**
* 线和面关系
* @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/），不得用于商业目标，基于本文批改后的作品务必以雷同的许可公布。