作為一個(gè)GISer，在日常WebGIS開(kāi)發(fā)中，會(huì)常用到的turf.js，這是一個(gè)地理空間分析的JAVAScript庫(kù)，經(jīng)常搭配各種GIS JS API使用，如leaflet、mapboxgl、openlayers等；在后臺(tái)Java開(kāi)發(fā)中，也有個(gè)比較強(qiáng)大的GIS庫(kù)，geotools，里面包含構(gòu)建一個(gè)完整的地理信息系統(tǒng)所需要的全部工具類；數(shù)據(jù)庫(kù)端常用是postgis擴(kuò)展，需要在postgres庫(kù)中引入使用。

然而在開(kāi)發(fā)某一些業(yè)務(wù)系統(tǒng)的時(shí)候，有些需求只需要調(diào)用某一個(gè)GIS算法，簡(jiǎn)單的幾行代碼即可完成，沒(méi)有必要去引用一個(gè)GIS類庫(kù)。

而且有些算法在這些常用的GIS類庫(kù)中沒(méi)有對(duì)應(yīng)接口，就比如在下文記錄的這幾種常用算法中，求垂足、判斷線和面的關(guān)系，在turf.js就沒(méi)有對(duì)應(yīng)接口。

下面文章中是我總結(jié)的一些常用GIS算法，這里統(tǒng)一用JavaScript語(yǔ)言實(shí)現(xiàn)，因?yàn)镴S代碼相對(duì)比較簡(jiǎn)潔，方便理解其中算法邏輯，也方便在瀏覽器下預(yù)覽效果。在具體應(yīng)用時(shí)可以根據(jù)具體需求，翻譯成Java、C#、Python等語(yǔ)言來(lái)使用。

文中代碼大部分為之前遇到需求時(shí)在網(wǎng)上搜索得到，然后自己根據(jù)具體需要做了優(yōu)化修改，通過(guò)這篇文章做個(gè)總結(jié)收集，也方便后續(xù)使用時(shí)查找。

1、常用算法

以下方法中傳參的點(diǎn)、線、面都是對(duì)應(yīng)geojson格式中coordinates，方便統(tǒng)一調(diào)用。geojson標(biāo)準(zhǔn)參考：
https://www.oschina.net/translate/geojson-spec

1.1、計(jì)算兩經(jīng)緯度點(diǎn)之間的距離

適用場(chǎng)景：測(cè)量

/**
* 計(jì)算兩經(jīng)緯度點(diǎn)之間的距離(單位：米)
* @param p1 起點(diǎn)的坐標(biāo)；[經(jīng)度,緯度]；例：[116.35,40.08]
* @param p2 終點(diǎn)的坐標(biāo)；[經(jīng)度,緯度]；例：[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、根據(jù)已知線段以及到起點(diǎn)距離，求目標(biāo)點(diǎn)坐標(biāo)

適用場(chǎng)景：封閉管段定位問(wèn)題點(diǎn)

/**
* 根據(jù)已知線段以及到起點(diǎn)距離（單位：米），求目標(biāo)點(diǎn)坐標(biāo)
* @param line 線段；[[經(jīng)度,緯度],[經(jīng)度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param dis 到起點(diǎn)距離（米）；Number；例：500
*
* @return point 返回坐標(biāo)
*/
function getLinePoint(line, dis) {
  var p1 = line[0]
  var p2 = line[1]
  var d = getDistance(p1, p2) // 計(jì)算兩經(jīng)緯度點(diǎn)之間的距離(單位：米)
  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、已知點(diǎn)、線段，求垂足

垂足可能在線段上，也可能在線段延長(zhǎng)線上。

適用場(chǎng)景：求垂足

/**
* 已知點(diǎn)、線段，求垂足
* @param line 線段；[[經(jīng)度,緯度],[經(jīng)度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param p 點(diǎn)；[經(jīng)度,緯度]；例：[116.35,40.08]
*
* @return point 返回垂足坐標(biāo)
*/
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、線段上距離目標(biāo)點(diǎn)最近的點(diǎn)

不同于上面求垂足方法，該方法求出的點(diǎn)肯定在線段上。

如果垂足在線段上，則最近的點(diǎn)就是垂足，如果垂足在線段延長(zhǎng)線上，則最近的點(diǎn)就是線段某一個(gè)端點(diǎn)。

適用場(chǎng)景：根據(jù)求出最近的點(diǎn)計(jì)算點(diǎn)到線段的最短距離

/**
* 線段上距離目標(biāo)點(diǎn)最近的點(diǎn)
* @param line 線段；[[經(jīng)度,緯度],[經(jīng)度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param p 點(diǎn)；[經(jīng)度,緯度]；例：[116.35,40.08]
*
* @return point 最近的點(diǎn)坐標(biāo)
*/
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、點(diǎn)緩沖

這里緩沖屬于測(cè)地線方法，由于這里并沒(méi)有嚴(yán)格的投影轉(zhuǎn)換體系，所以與標(biāo)準(zhǔn)的測(cè)地線緩沖還有些許誤差，不過(guò)經(jīng)測(cè)試，半徑100KM內(nèi)，誤差基本可以忽略。具體緩沖類型可看下之前的文章你真的會(huì)用PostGIS中的buffer緩沖嗎？

適用場(chǎng)景：根據(jù)點(diǎn)和半徑畫圓

/**
* 點(diǎn)緩沖
* @param center 中心點(diǎn)；[經(jīng)度,緯度]；例：[116.35,40.08]
* @param radius 半徑（米）；Number；例：5000
* @param vertices 返回圓面點(diǎn)的個(gè)數(shù)；默認(rèn)64；Number；例：32
*
* @return coords 面的坐標(biāo)
*/
function bufferPoint(center, radius, vertices) {
  if (!vertices) vertices = 64;
  var coords = []
  // 111319.55：在赤道上1經(jīng)度差對(duì)應(yīng)的距離，111133.33：在經(jīng)線上1緯度差對(duì)應(yīng)的距離
  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、點(diǎn)和面關(guān)系

該方法采用射線法思路實(shí)現(xiàn)。（了解射線法可參考：
https://blog.csdn.net/qq_27161673/article/details/52973866）

這里已經(jīng)考慮到環(huán)狀多邊形的情況。

適用場(chǎng)景：判斷點(diǎn)是否在面內(nèi)

/**
* 點(diǎn)和面關(guān)系
* @param point 點(diǎn)；[經(jīng)度,緯度]；例：[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 點(diǎn)和面關(guān)系；0:多邊形外，1：多邊形內(nèi)，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;
  }
}


/**
* 點(diǎn)和面關(guān)系
* @param point 點(diǎn)
* @param ring 單個(gè)閉合面的坐標(biāo)
*
* @return inside 點(diǎn)和面關(guān)系；0:多邊形外，1：多邊形內(nèi)，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
    }
    // 判斷點(diǎn)與線段的相對(duì)位置，0為在線段上，>0 點(diǎn)在左側(cè)，<0 點(diǎn)在右側(cè)
    if (isLeft(point, [ring[i], ring[j]]) === 0) {
      return 2; // 點(diǎn)在多邊形邊上
    } else {
      if ((yi > y) !== (yj > y)) { // 垂直方向目標(biāo)點(diǎn)在yi、yj之間
        // 求目標(biāo)點(diǎn)在當(dāng)前線段上的x坐標(biāo)。 由于JS小數(shù)運(yùn)算后會(huì)轉(zhuǎn)換為精確15位的float，因此需要去一下精度
        var xx = Number(((xj - xi) * (y - yi) / (yj - yi) + xi).toFixed(10))
        if (x <= xx) { // 目標(biāo)點(diǎn)水平射線與當(dāng)前線段有交點(diǎn)
          inside = !inside;
        }
      }
    }
  }
  return Number(inside);
}


/**
* 判斷點(diǎn)與線段的相對(duì)位置
* @param point 目標(biāo)點(diǎn)
* @param line 線段
*
* @return isLeft，點(diǎn)與線段的相對(duì)位置，0為在線段上，>0 p在左側(cè)，<0 p在右側(cè)
*/
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小數(shù)運(yùn)算后會(huì)轉(zhuǎn)換為精確15位的float，因此需要去一下精度
  return Number(isLeft.toFixed(10))
}

1.7、線段與線段的關(guān)系

適用場(chǎng)景：判斷線和線的關(guān)系

/**
* 線段與線段的關(guān)系
* @param line1 線段；[[經(jīng)度,緯度],[經(jīng)度,緯度]]；例：[[116.01,40.01],[116.52,40.01]]
* @param line2 線段；[[經(jīng)度,緯度],[經(jīng)度,緯度]]；例：[[116.33,40.21],[116.36,39.76]]
*
* @return intersect 線段與線段的關(guān)系；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]

  //快速排斥：
  //兩個(gè)線段為對(duì)角線組成的矩形，如果這兩個(gè)矩形沒(méi)有重疊的部分，那么兩條線段是不可能出現(xiàn)重疊的

  //這里的確如此，這一步是判定兩矩形是否相交
  //1.線段ab的低點(diǎn)低于cd的最高點(diǎn)（可能重合）
  //2.cd的最左端小于ab的最右端（可能重合）
  //3.cd的最低點(diǎn)低于ab的最高點(diǎn)（加上條件1，兩線段在豎直方向上重合）
  //4.ab的最左端小于cd的最右端（加上條件2，兩直線在水平方向上重合）
  //綜上4個(gè)條件，兩條線段組成的矩形是重合的
  //特別要注意一個(gè)矩形含于另一個(gè)矩形之內(nèi)的情況
  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
  }

  // 判斷點(diǎn)與線段的相對(duì)位置，0為在線段上，>0 點(diǎn)在左側(cè)，<0 點(diǎn)在右側(cè)
  if (isLeft(line1[0], line2) === 0 || isLeft(line1[1], line2) === 0) {
    return 2
  }

  //跨立實(shí)驗(yàn)：
  //如果兩條線段相交，那么必須跨立，就是以一條線段為標(biāo)準(zhǔn)，另一條線段的兩端點(diǎn)一定在這條線段的兩段
  //也就是說(shuō)a b兩點(diǎn)在線段cd的兩端，c d兩點(diǎn)在線段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、線和面關(guān)系

適用場(chǎng)景：判斷線與面的關(guān)系

該方法考慮到環(huán)狀多邊形的情況，且把相切情況分為了內(nèi)切和外切。

參考鏈接：
https://www.cnblogs.com/xiaozhi_5638/p/4165353.html

/**
* 線和面關(guān)系
* @param line 線段；[[經(jīng)度,緯度],[經(jīng)度,緯度]]；例：[[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 線和面關(guān)系；0:相離，1：相交，2：包含，3：內(nèi)切，4：外切
*/
function intersectLineAndPolygon(line, polygon) {
  var isTangent = false
  var isInNum = 0
  var intersect = 0
  for (var i = 0; i < polygon.length; i++) {
    // 線和面關(guān)系；0:相離，1：相交，2：包含，3：內(nèi)切，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 // 內(nèi)切
    } else {
      return 2 // 包含
    }
  }
}


/**
* 線和面關(guān)系
* @param line 線段
* @param ring 單面
*
* @return intersect 線和面關(guān)系；0:相離，1：相交，2：包含，3：內(nèi)切，4：外切
*/
function intersectLineAndRing(line, ring) {
  var inserset = 0
  var isTangent = false
  var inserset1 = pointInRing(line[0], ring) // 點(diǎn)和面關(guān)系；0:多邊形外，1：多邊形內(nèi)，2：多邊形邊上
  var inserset2 = pointInRing(line[1], ring) // 點(diǎn)和面關(guān)系；0:多邊形外，1：多邊形內(nèi)，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]]
    // 目標(biāo)線段與當(dāng)前線段的關(guān)系；0:相離，1：相交，2：相切
    var intersectLine = intersectLineAndLine(line, line2)
    if (intersectLine == 1) {
      return 1 // 相交
    }
  }
  return inserset
}

1.9、geojson 面轉(zhuǎn)線

適用場(chǎng)景：只有g(shù)eojson面數(shù)據(jù)，獲取線的邊界

/**
* 面轉(zhuǎn)線
* @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('請(qǐng)確認(rèn)輸入?yún)?shù)為geojson格式面數(shù)據(jù)！')
      return null
    }
    polylineGeoJson.features[i].geometry.type = 'MultiLineString' //面轉(zhuǎn)線
    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

關(guān)注《GIS兵器庫(kù)》， 只給你網(wǎng)上搜不到的GIS知識(shí)技能。

本文章采用 知識(shí)共享署名-非商業(yè)性使用-相同方式共享 4.0 國(guó)際許可協(xié)議 進(jìn)行許可。歡迎轉(zhuǎn)載、使用、重新發(fā)布，但務(wù)必保留文章署名《GIS兵器庫(kù)》（包含鏈接：
http://gisarmory.xyz/blog/），不得用于商業(yè)目的，基于本文修改后的作品務(wù)必以相同的許可發(fā)布。