您当前的位置:首页 > 电脑百科 > 程序开发 > 算法

GIS常用算法

时间:2021-09-13 16:28:18  来源:  作者:GIS兵器库

作为一个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

GIS常用算法

 

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/),不得用于商业目的,基于本文修改后的作品务必以相同的许可发布。



Tags:GIS   点击:()  评论:()
声明:本站部分内容及图片来自互联网,转载是出于传递更多信息之目的,内容观点仅代表作者本人,如有任何标注错误或版权侵犯请与我们联系(Email:2595517585@qq.com),我们将及时更正、删除,谢谢。
▌相关推荐
SimpleAI.人工智能、机器学习、深度学习还是遥不可及?来这里看看吧~ 从基本的概念、原理、公式,到用生动形象的例子去理解,到动手做实验去感知,到著名案例的学习,到用所学来实现...【详细内容】
2021-10-19  Tags: GIS  点击:(47)  评论:(0)  加入收藏
作为一个GISer,在日常WebGIS开发中,会常用到的turf.js,这是一个地理空间分析的JavaScript库,经常搭配各种GIS JS API使用,如leaflet、mapboxgl、openlayers等;在后台Java开发中,也...【详细内容】
2021-09-13  Tags: GIS  点击:(65)  评论:(0)  加入收藏
本文是ArcGIS完成空间插值的文字版教程和视频教程,其中文字版教程部分转载自TOOL尖尖的木木作者,视频版教程来自二狮兄发布在百度文库的视频教程。其中上集高清在线视频地址为...【详细内容】
2020-08-20  Tags: GIS  点击:(294)  评论:(0)  加入收藏
▌简易百科推荐
前言Kafka 中有很多延时操作,比如对于耗时的网络请求(比如 Produce 是等待 ISR 副本复制成功)会被封装成 DelayOperation 进行延迟处理操作,防止阻塞 Kafka请求处理线程。Kafka...【详细内容】
2021-12-27  Java技术那些事    Tags:时间轮   点击:(1)  评论:(0)  加入收藏
博雯 发自 凹非寺量子位 报道 | 公众号 QbitAI在炼丹过程中,为了减少训练所需资源,MLer有时会将大型复杂的大模型“蒸馏”为较小的模型,同时还要保证与压缩前相当的结果。这就...【详细内容】
2021-12-24  量子位    Tags:蒸馏法   点击:(9)  评论:(0)  加入收藏
分稀疏重建和稠密重建两类:稀疏重建:使用RGB相机SLAMOrb-slam,Orb-slam2,orb-slam3:工程地址在: http://webdiis.unizar.es/~raulmur/orbslam/ DSO(Direct Sparse Odometry)因为...【详细内容】
2021-12-23  老师明明可以靠颜值    Tags:算法   点击:(7)  评论:(0)  加入收藏
1. 基本概念希尔排序又叫递减增量排序算法,它是在直接插入排序算法的基础上进行改进而来的,综合来说它的效率肯定是要高于直接插入排序算法的;希尔排序是一种不稳定的排序算法...【详细内容】
2021-12-22  青石野草    Tags:希尔排序   点击:(6)  评论:(0)  加入收藏
ROP是一种技巧,我们对execve函数进行拼凑来进行system /bin/sh。栈迁移的特征是溢出0x10个字符,在本次getshell中,还碰到了如何利用printf函数来进行canary的泄露。ROP+栈迁移...【详细内容】
2021-12-15  星云博创    Tags:栈迁移   点击:(19)  评论:(0)  加入收藏
一、什么是冒泡排序1.1、文字描述冒泡排序是一种简单的排序算法。它重复地走访要排序的数列,一次比较两个元素,如果他们的顺序错误就把他们交换过来。走访数列的工作是重复地...【详细内容】
2021-12-15    晓掌柜丶韶华  Tags:排序算法   点击:(16)  评论:(0)  加入收藏
在了解golang的map之前,我们需要了解哈希这个概念。哈希表,又称散列表(Hash table),是根据键(key)而直接访问在内存储存位置的数据结构。也就是说,它通过计算出一个键值的函数,将...【详细内容】
2021-12-07  一棵梧桐木    Tags:哈希表   点击:(13)  评论:(0)  加入收藏
前面文章在谈论分布式唯一ID生成的时候,有提到雪花算法,这一次,我们详细点讲解,只讲它。SnowFlake算法据国家大气研究中心的查尔斯&middot;奈特称,一般的雪花大约由10^19个水分子...【详细内容】
2021-11-17  小心程序猿QAQ    Tags:雪花算法   点击:(24)  评论:(0)  加入收藏
导读:在大数据时代,对复杂数据结构中的各数据项进行有效的排序和查找的能力非常重要,因为很多现代算法都需要用到它。在为数据恰当选择排序和查找策略时,需要根据数据的规模和类型进行判断。尽管不同策略最终得到的结果完...【详细内容】
2021-11-04  华章科技    Tags:排序算法   点击:(37)  评论:(0)  加入收藏
这是我在网上找的资源的一个总结,会先给出一个我看了觉得还行的关于算法的讲解,再配上实现的代码: Original author: Bill_Hoo Original Address: http://blog.sina.com.cn/s/bl...【详细内容】
2021-11-04  有AI野心的电工和码农    Tags: KMP算法   点击:(36)  评论:(0)  加入收藏
最新更新
栏目热门
栏目头条