GIS算法总结（求经纬度，求距离，求方位角） |
您所在的位置：网站首页 › 经纬度距离算法 › GIS算法总结（求经纬度，求距离，求方位角） |
|
公司项目中用到的一些GIS算法，记录一下。所有经纬度都是基于WGS-84标准，为Google Earth上的经纬度。 1.已知两个点的经纬度，求它们的直线距离 python代码实现 from math import * def getDistance(latA, lonA, latB, lonB): ra = 6378140 # radius of equator: meter rb = 6356755 # radius of polar: meter flatten = (ra - rb) / ra # Partial rate of the earth # change angle to radians radLatA = radians(latA) radLonA = radians(lonA) radLatB = radians(latB) radLonB = radians(lonB) pA = atan(rb / ra * tan(radLatA)) pB = atan(rb / ra * tan(radLatB)) x = acos(sin(pA) * sin(pB) + cos(pA) * cos(pB) * cos(radLonA - radLonB)) c1 = (sin(x) - x) * (sin(pA) + sin(pB))**2 / cos(x / 2)**2 c2 = (sin(x) + x) * (sin(pA) - sin(pB))**2 / sin(x / 2)**2 dr = flatten / 8 * (c1 - c2) distance = ra * (x + dr) return distance distance = getDistance(30.25482684590,120.01375138640,30.25582684590,120.01375138640) print(distance, '米')2.已知两个点的经纬度，求它们的方位角 python代码实现 def getDegree(latA, lonA, latB, lonB): """ Args: point p1(latA, lonA) point p2(latB, lonB) Returns: bearing between the two GPS points, default: the basis of heading direction is north """ radLatA = radians(latA) radLonA = radians(lonA) radLatB = radians(latB) radLonB = radians(lonB) dLon = radLonB - radLonA y = sin(dLon) * cos(radLatB) x = cos(radLatA) * sin(radLatB) - sin(radLatA) * cos(radLatB) * cos(dLon) brng = degrees(atan2(y, x)) brng = (brng + 360) % 360 return brng brng = getDegree(30.25482684590,120.01375138640,30.25582684590,120.01375138640) print(brng)3.已知一个点的经纬度，方位角，距离，求另一点的经纬度 python代码实现 from math import * def rad(d): return d * pi / 180.0 def deg(x): return x * 180 / pi def getLonAndLat(lat, lng, brng, dist): a = 6378137 b = 6356752.3142 f = 1.0 / 298.257223563 lon1 = lng * 1 lat1 = lat * 1 s = dist alpha1 = rad(brng) sinAlpha1 = sin(alpha1) cosAlpha1 = cos(alpha1) tanU1 = (1 - f) * tan(rad(lat1)) cosU1 = 1 / sqrt((1 + tanU1 * tanU1)) sinU1 = tanU1 * cosU1 sigma1 = atan2(tanU1, cosAlpha1) sinAlpha = cosU1 * sinAlpha1 cosSqAlpha = 1 - sinAlpha * sinAlpha uSq = cosSqAlpha * (a * a - b * b) / (b * b) A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))) B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))) sigma = s / (b * A) sigmaP = 2 * pi while abs(sigma - sigmaP) > 1e-12: cos2SigmaM = cos(2 * sigma1 + sigma) sinSigma = sin(sigma) cosSigma = cos(sigma) deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * ( -3 + 4 * cos2SigmaM * cos2SigmaM))) sigmaP = sigma sigma = s / (b * A) + deltaSigma tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1 lat2 = atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, (1 - f) * sqrt(sinAlpha * sinAlpha + tmp * tmp)) lambdA = atan2(sinSigma * sinAlpha1, cosU1*cosSigma - sinU1 * sinSigma * cosAlpha1) C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)) L = lambdA - (1-C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM+C * cosSigma * (-1+2 * cos2SigmaM * cos2SigmaM))) revAz = atan2(sinAlpha, -tmp) lngLatObj = {'latitude': deg(lat2), 'longitude': lon1 + deg(L), } return lngLatObj print(getLonAndLat(30.2548268459,120.0137513864,0,110.85670061685299)) JavaScript代码实现 /** * 度换成弧度 * @param {Float} d 度 * @return {[Float} 弧度 */ function rad (d){ return d * Math.PI / 180.0; }; /** * 弧度换成度 * @param {Float} x 弧度 * @return {Float} 度 */ function deg (x){ return x*180/Math.PI; }; /** * * @param {*} lng 经度 113.3960698 * @param {*} lat 纬度 22.941386 * @param {*} brng 方位角 45 ---- 正北方：000°或360° 正东方：090° 正南方：180° 正西方：270° * @param {*} dist 距离 9000 * */ function getLonAndLat (lat,lng,brng,dist){ //大地坐标系资料WGS-84 长半径a=6378137 短半径b=6356752.3142 扁率f=1/298.2572236 var a=6378137; var b=6356752.3142; var f=1/298.257223563; var lon1 = lng*1; var lat1 = lat*1; var s = dist; var alpha1 = rad(brng); var sinAlpha1 = Math.sin(alpha1); var cosAlpha1 = Math.cos(alpha1); var tanU1 = (1-f) * Math.tan(rad(lat1)); var cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1*cosU1; var sigma1 = Math.atan2(tanU1, cosAlpha1); var sinAlpha = cosU1 * sinAlpha1; var cosSqAlpha = 1 - sinAlpha*sinAlpha; var uSq = cosSqAlpha * (a*a - b*b) / (b*b); var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); var sigma = s / (b*A), sigmaP = 2*Math.PI; while (Math.abs(sigma-sigmaP) > 1e-12) { var cos2SigmaM = Math.cos(2*sigma1 + sigma); var sinSigma = Math.sin(sigma); var cosSigma = Math.cos(sigma); var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)- B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM))); sigmaP = sigma; sigma = s / (b*A) + deltaSigma; } var tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1; var lat2 = Math.atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1, (1-f)*Math.sqrt(sinAlpha*sinAlpha + tmp*tmp)); var lambda = Math.atan2(sinSigma*sinAlpha1, cosU1*cosSigma - sinU1*sinSigma*cosAlpha1); var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha)); var L = lambda - (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM))); var revAz = Math.atan2(sinAlpha, -tmp); // final bearing var lngLatObj = {latitude:deg(lat2),longitude:lon1+deg(L),} return lngLatObj; } //调用 console.log(getLonAndLat(30.2548268459,120.0137513864,0,110.85670061685299)); |
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |
| CopyRight 2018-2019 实验室设备网 版权所有 win10的实时保护怎么永久关闭 |