经纬度相关公式及实现
需求
- 已知两点经纬度，求两点间的距离
- 已知两点经纬度，求两点间的方位角
- 已知一点经纬度为中心、r为半径，求出该经纬度区间
- 已知一点的经纬度、距离，方位角，求另一个点的经纬度
公式
-
已知两点经纬度，求两点间的距离:
AB=Rarccos[sin(wA)sin(wB)+cos(wA)cos(wB)cos(jA-jB)]
-
已知两点经纬度，求两点间的方位角
rarccos[cos(y1)cos(y2)cos(x1-x2)+sin(y1)sin(y2)]
这里得到的是弧度，我们需要将弧度转换成我们需要的方位角，这里需要考虑象限问题
纬度是横轴，经度时纵轴
B点在第一象限，Azimuth=A；
B点在第四象限，Azimuth=360+A;
B点在第二三象限，Azimuth=180-A -
已知一点经纬度为中心、r为半径，求出该经纬度区间
-
已知一点的经纬度、距离，方位角，求另一个点的经纬度
实现代码
import java.util.HashMap;
import java.util.Map;
public class DistanceUtil {
/*
* 大地坐标系资料WGS-84 长半径a=6378137 短半径b=6356752.3142 扁率f=1/298.2572236
*/
/** 长半径a=6378137 */
private static double EARTH_RADIUS = 6378137;
/** 短半径b=6356752.3142 */
private static double b = 6356752.3142;
/** 扁率f=1/298.2572236 */
private static double f = 1 / 298.2572236;
private static double rad(double d) {
return d * Math.PI / 180.0;
}
/**
* 弧度换成度
*
* @param x
*弧度
* @return 度
*/
private static double deg(double x) {
return x * 180 / Math.PI;
}
/**
* 根据两个位置的经纬度，来计算两地的距离（单位为KM）
* 参数为String类型
* @param lat1 用户经度
* @param lon1 用户纬度
* @param lat2 商家经度
* @param lon2 商家纬度
* @return
*/
public static double getDistance(int lat1, int lon1, int lat2, int lon2) {
double tmplat1 = lat1/1000000.0;
double tmplon1 = lon1/1000000.0;
double tmplat2 = lat2/1000000.0;
double tmplon2 = lon2/1000000.0;;
double radLat1 = rad(tmplat1);
double radLat2 = rad(tmplat2);
double difference = radLat1 - radLat2;
double mdifference = rad(tmplon1) - rad(tmplon2);
double distance = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(difference / 2), 2)
+ Math.cos(radLat1) * Math.cos(radLat2)
* Math.pow(Math.sin(mdifference / 2), 2)));
distance = distance * EARTH_RADIUS;
distance = Math.round(distance * 10000) / 10000;
return distance;
}
/**
* 获取当前用户一定距离以内的经纬度值,距离单位是米
* 单位米 return minLat
* 最小经度 minLon
* 最小纬度 maxLat
* 最大经度 maxLon
* 最大纬度 minLat
*/
public static Map<String,Integer> getAround(int lat, int lon, double raidus) {
Map<String,Integer> map = new HashMap();
double latitude = lat/1000000.0;// 传值给经度
double longitude = lon/1000000.0;// 传值给纬度
Double degree = (24901 * 1609) / 360.0; // 获取每度
double raidusMile = raidus;
Double mpdLon = Double.parseDouble((degree * Math.cos(latitude * (Math.PI / 180))+"").replace("-", ""));
Double dpmLon = 1 / mpdLon;
Double radiusLon = dpmLon * raidusMile;
//获取最小经度
double minLat = latitude - radiusLon;
// 获取最大经度
double maxLat = latitude + radiusLon;
Double dpmLat = 1 / degree;
Double radiusLat = dpmLat * raidusMile;
// 获取最小纬度
double minLon = longitude - radiusLat;
// 获取最大纬度
double maxLon = longitude + radiusLat;
map.put("minLat", new Double(minLat*1000000).intValue());
map.put("maxLat", new Double(maxLat*1000000).intValue());
map.put("minLon", new Double(minLon*1000000).intValue());
map.put("maxLon", new Double(maxLon*1000000).intValue());
return map;
}
/**
* 已知两个点的坐标，求两个的方位角
* @param lat1 起始点
* @param lon1
* @param lat2 结束点
* @param lon2
* @return
*/
public static double getAngle(int lat1, int lon1, int lat2, int lon2) {
double tmplat1 = lat1/1000000.0;
double tmplon1 = lon1/1000000.0;
double tmplat2 = lat2/1000000.0;
double tmplon2 = lon2/1000000.0;
double dRotateAngle = Math.atan2(Math.abs(tmplat1 - tmplat2), Math.abs(tmplon1 - tmplon2));
/*
纬度是横轴，经度时纵轴
*/
if (tmplat2 >= tmplat1){
if (tmplon2 >=tmplon1){
//第一象限
} else {
// 第四象限
dRotateAngle =2 * Math.PI - dRotateAngle;
}
} else{
if (tmplon2 >= tmplon1) {
//第二象限
dRotateAngle = Math.PI - dRotateAngle;
} else {
//第三象限
dRotateAngle = Math.PI + dRotateAngle;
}
}
dRotateAngle = dRotateAngle * 180 / Math.PI;
return dRotateAngle;
}
/**
* 求两经纬度方向角
*
*@param lat1
*第一点的纬度
* @param lon1
*第一点的经度
* @param lat2
*第二点的纬度
@param lon2
*第二点的经度
* @return 方位角，角度（单位：°）
* */
public static double getAzimuth(int lat1,int lon1,int lat2,int lon2) {
double tmplat1 = lat1/1000000.0;
double tmplon1 = lon1/1000000.0;
double tmplat2 = lat2/1000000.0;
double tmplon2 = lon2/1000000.0;
tmplat1 = rad(tmplat1);
tmplat2 = rad(tmplat2);
tmplon1 = rad(tmplon1);
tmplon2 = rad(tmplon2);
double azimuth = Math.sin(tmplat1) * Math.sin(tmplat2) + Math.cos(tmplat1)
* Math.cos(tmplat2) * Math.cos(tmplon2 - tmplon1);
azimuth = Math.sqrt(1 - azimuth * azimuth);
azimuth = Math.cos(tmplat2) * Math.sin(tmplon2 - tmplon1) / azimuth;
azimuth = Math.asin(azimuth) * 180 / Math.PI;
if (Double.isNaN(azimuth)) {
if (lon1 < lon2) {
azimuth = 90.0;
} else {
azimuth = 270.0;
}
}
/*
纬度是横轴，经度时纵轴
B点在第一象限，Azimuth=A；
B在第四象限，Azimuth=360+A;
B在第二三象限，Azimuth=180-A
*/
if (tmplat2 >= tmplat1){
if (tmplon2 >=tmplon1){
//第一象限
} else {
// 第四象限
azimuth =360 + azimuth;
}
} else{
if (tmplon2 >= tmplon1) {
//第二象限
azimuth = 180 - azimuth;
} else {
//第三象限
azimuth = 180 - azimuth;
}
}
return azimuth;
}
/**
* 计算另一点经纬度
*
* @param lon
*经度
* @param lat
*维度
* @param brng
*方位角
* @param dist
*距离（米）
*/
public static Map<String,Integer> computerThatLonLat(int lat,int lon, double brng, double dist) {
double tmplon = lon / 1000000.0;
double tmplat = lat / 1000000.0;
double alpha1 = rad(brng);
double sinAlpha1 = Math.sin(alpha1);
double cosAlpha1 = Math.cos(alpha1);
double tanU1 = (1 - f) * Math.tan(rad(tmplat));
double cosU1 = 1 / Math.sqrt((1 + tanU1 * tanU1));
double sinU1 = tanU1 * cosU1;
double sigma1 = Math.atan2(tanU1, cosAlpha1);
double sinAlpha = cosU1 * sinAlpha1;
double cosSqAlpha = 1 - sinAlpha * sinAlpha;
double uSq = cosSqAlpha * (EARTH_RADIUS * EARTH_RADIUS - b * b) / (b * b);
double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
double cos2SigmaM=0;
double sinSigma=0;
double cosSigma=0;
double sigma = dist / (b * A), sigmaP = 2 * Math.PI;
while (Math.abs(sigma - sigmaP) > 1e-12) {
cos2SigmaM = Math.cos(2 * sigma1 + sigma);
sinSigma = Math.sin(sigma);
cosSigma = Math.cos(sigma);
double 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 = dist / (b * A) + deltaSigma;
}
double tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1;
double lat2 = Math.atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1,
(1 - f) * Math.sqrt(sinAlpha * sinAlpha + tmp * tmp));
double lambda = Math.atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1);
double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
double L = lambda - (1 - C) * f * sinAlpha
* (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
double revAz = Math.atan2(sinAlpha, -tmp); // final bearing
Map<String,Integer> map = new HashMap<String, Integer>();
map.put("lon",+new Double((tmplon+deg(L))*1000000).intValue());
map.put("lat",+new Double(deg(lat2)*1000000).intValue());
return map;
}
/**
* 判断是否漂移
* @return
*/
/**
* 判断是否漂移，如果漂移则返回true,否则返回false
* @param lat1 单位是10^-6
* @param lon1 单位是10^-6
* @param raidus 半径单位是米
* @param lat2
* @param lon2
* @return
*/
public static boolean isDrift(int lat1,int lon1,int raidus,int lat2,int lon2){
//中国纬度范围:4°~53°
if(lat1<4000000 || lat1 >53000000 ){
return true;
}
//中国经度范围:73°~135°
if(lon1<73000000 || lon1 >135000000){
return true;
}
Map<String,Integer> map = getAround(lat1,lon1,raidus);
if(lat2>map.get("maxLat") || lat2<map.get("minLat")){
return true;
}
if(lon2>map.get("maxLon") || lon2<map.get("minLon")){
return true;
}
return false;
}
}
