共计 8533 个字符,预计需要花费 22 分钟才能阅读完成。
作为一个 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/),不得用于商业目标,基于本文批改后的作品务必以雷同的许可公布。