GIS算法总结(求经纬度,求距离,求方位角)
程序员文章站
2022-06-11 16:10:36
...
公司项目中用到的一些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));
上一篇: 对CIS的认识
下一篇: arcgis文档地址