GPS坐标区域判定2
1.两点之间距离
通过手动算法计算：
public static double getEarthDistance(double lat1,double lon1,
double lat2, double lon2){
double φ1 = Math.toRadians(lat1);
double φ2 = Math.toRadians(lat2);
double Δφ = Math.toRadians(lat2-lat1);
double Δλ = Math.toRadians(lon2-lon1);
double a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
Math.cos(φ1) * Math.cos(φ2) *
Math.sin(Δλ/2) * Math.sin(Δλ/2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
double d = EARTH_RADIUS_KM * c;
return d;
}
通过android中现成的方法:
private static void computeDistanceAndBearing(double lat1, double lon1,
double lat2, double lon2, float[] results) {
// Based on http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
// using the "Inverse Formula" (section 4)
int MAXITERS = 20;
// Convert lat/long to radians
lat1 *= Math.PI / 180.0;
lat2 *= Math.PI / 180.0;
lon1 *= Math.PI / 180.0;
lon2 *= Math.PI / 180.0;
double a = 6378137.0; // WGS84 major axis
double b = 6356752.3142; // WGS84 semi-major axis
double f = (a - b) / a;
double aSqMinusBSqOverBSq = (a * a - b * b) / (b * b);
double L = lon2 - lon1;
double A = 0.0;
double U1 = Math.atan((1.0 - f) * Math.tan(lat1));
double U2 = Math.atan((1.0 - f) * Math.tan(lat2));
double cosU1 = Math.cos(U1);
double cosU2 = Math.cos(U2);
double sinU1 = Math.sin(U1);
double sinU2 = Math.sin(U2);
double cosU1cosU2 = cosU1 * cosU2;
double sinU1sinU2 = sinU1 * sinU2;
double sigma = 0.0;
double deltaSigma = 0.0;
double cosSqAlpha = 0.0;
double cos2SM = 0.0;
double cosSigma = 0.0;
double sinSigma = 0.0;
double cosLambda = 0.0;
double sinLambda = 0.0;
double lambda = L; // initial guess
for (int iter = 0; iter < MAXITERS; iter++) {
double lambdaOrig = lambda;
cosLambda = Math.cos(lambda);
sinLambda = Math.sin(lambda);
double t1 = cosU2 * sinLambda;
double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
double sinSqSigma = t1 * t1 + t2 * t2; // (14)
sinSigma = Math.sqrt(sinSqSigma);
cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
sigma = Math.atan2(sinSigma, cosSigma); // (16)
double sinAlpha = (sinSigma == 0) ? 0.0 :
cosU1cosU2 * sinLambda / sinSigma; // (17)
cosSqAlpha = 1.0 - sinAlpha * sinAlpha;
cos2SM = (cosSqAlpha == 0) ? 0.0 :
cosSigma - 2.0 * sinU1sinU2 / cosSqAlpha; // (18)
double uSquared = cosSqAlpha * aSqMinusBSqOverBSq; // defn
A = 1 + (uSquared / 16384.0) * // (3)
(4096.0 + uSquared *
(-768 + uSquared * (320.0 - 175.0 * uSquared)));
double B = (uSquared / 1024.0) * // (4)
(256.0 + uSquared *
(-128.0 + uSquared * (74.0 - 47.0 * uSquared)));
double C = (f / 16.0) *
cosSqAlpha *
(4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
double cos2SMSq = cos2SM * cos2SM;
deltaSigma = B * sinSigma * // (6)
(cos2SM + (B / 4.0) *
(cosSigma * (-1.0 + 2.0 * cos2SMSq) -
(B / 6.0) * cos2SM *
(-3.0 + 4.0 * sinSigma * sinSigma) *
(-3.0 + 4.0 * cos2SMSq)));
lambda = L +
(1.0 - C) * f * sinAlpha *
(sigma + C * sinSigma *
(cos2SM + C * cosSigma *
(-1.0 + 2.0 * cos2SM * cos2SM))); // (11)
double delta = (lambda - lambdaOrig) / lambda;
if (Math.abs(delta) < 1.0e-12) {
break;
}
}
float distance = (float) (b * A * (sigma - deltaSigma));
results[0] = distance;
if (results.length > 1) {
float initialBearing = (float) Math.atan2(cosU2 * sinLambda,
cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);
initialBearing *= 180.0 / Math.PI;
results[1] = initialBearing;
if (results.length > 2) {
float finalBearing = (float) Math.atan2(cosU1 * sinLambda,
-sinU1 * cosU2 + cosU1 * sinU2 * cosLambda);
finalBearing *= 180.0 / Math.PI;
results[2] = finalBearing;
}
}
}
2.判断一点是否在不规则多边形中
判断一点是否在多边形中，常用的是射线法，个人觉得也是最简单的一种，因为比较角度和面积都会牵扯到正负的判断，而射线法则不需要。规则：
以给定的点为起点，向任意方向做一条射线。如果射线与给定的多边形的交点为奇数则该点在多边形内部（点不再多边形边上且不与顶点重合），为偶数则为点在多边形外（0是偶数）。
单单这样看的话，也并不简单，首先应该先建立坐标系，然后取一条过点的射线，列出函数表达式，然后因为已知多边形各各点的坐标，求出每条边的函数，遍历每条边，计算与射线的交点并计算交点是否在多边形的边上。如果在，算作一个交点。最后算出总的交点个数即可。
特殊情况：当所选取的射线与一条边重合时，如果点不在重合线段上，算作两个交点，如果在，算作一个交点，并且排除相邻两边的计算，排除以上情况后，如果射线过多边形的一个顶点，那么依旧算作两个交点，并且排除过该顶点的另一条边的计算。其余的按照正常计算即可。
具体问题具体分析，gps坐标是地理坐标，是在一个球体上表示，如果求一个点是否在指定图形区域内，如果要在笛卡尔坐标系中计算，就需要将坐标投影到平面上，可以用web墨卡托坐标系，在平面坐标系中求得多边形距离给定点最近的一点的坐标，将此点再次转化为gps坐标，然后计算给定点与该点的距离即可。
之所以做两次转化是因为平面坐标系计算距离不准确，而且维度不同误差很大，所以平面坐标系将点到多边形转化为点到点的计算，然后可以在球面坐标系中用现成的方法来计算
但是，以上方法还是太麻烦了，因为是应用计算，完全不必要，只需要取一个特殊的射线即可，我取的是 给定点指向北极点（将极点看作无穷远）的射线。那么只需要3步骤即可：
- 1.判断射线是否与多边形边重合
- 2.判断给定点是否与多边形顶点重合
- 3.判断多边形每条边是否与射线相交（遍历，以下是对多边形每条线段的处理）
- 3.1.线段两个端点的纬度如果大于给定点，经度如果一个大于一个小于给定点则相交，记为一个交点。
- 3.2.如果线段两端点经度一个大于一个小于给定点，同时纬度一个大于一个小于给定点，那么就需要通过平行相交的比例关系判断与射线所在直线的交点是否是在射线上。如果在，记为一个交点。
具体算法:
/**
* 判断点是否在图形区域内
* 原理:射线法，取点所在的经线作为射线来判断
*
* 特殊情况：
* 1.射线与n条边重合
* 2.射线与经过n个顶点
* @param point 点的经纬度，0/1 latitude/ longitude 纬度/经度
* @param polygon 多边形的点集合，0/1 latitude/ longitude 纬度/经度
* */
public static boolean isPointInPolygon(double[] point,double[][] polygon){
/**参考射线取点所在经线自该点向南极发射的一条射线*/
int node = 0;
for (int i = 0; i < polygon.length; i++) {
int next = i+1;
if (i==polygon.length-1)
next = 0;
/**首先判断边是否与参考经线重合*/
if(point[1]==polygon[i][1]&&polygon[i][1]==polygon[next][1]){
if((point[0]>=polygon[i][0]&&point[0]<=polygon[next][0])||(point[0]<=polygon[i][0]&&point[0]>=polygon[next][0])){
/**在线段内*/
return true;
}else if(point[0]<polygon[i][0]&&point[0]<polygon[next][0]){
/**在线段外的射线上*/
System.out.println("~~点在线段延长线上:"+i);
node+=2;
}else{
/**在线段外，并且不在射线上*/
}
}else if(point[1]==polygon[i][1]){
/**判断顶点是否在参考经线之上*/
if(point[0]==polygon[i][0]){
/**点与图形顶点重合*/
return true;
}else if(point[0]<polygon[i][0]){
node+=2;
}
}else{
/**判断边是否与射线相交*/
/**这种计算的先决条件是线段两点必须分布在当前位置点所在经线的两端*/
if((polygon[i][1]>point[1]&&polygon[next][1]>point[1])||(polygon[i][1]<point[1]&&polygon[next][1]<point[1])){
/**参考经线不会经过该线段*/
continue;
}else{
/**判断纬度是否需要计算*/
if(polygon[i][0]>point[0]&& polygon[next][0]>point[0]){
System.out.println("~~纬度判断:"+i);
node++;
}else if(polygon[i][0]<point[0]&& polygon[next][0]>point[0]){
/***/
double tempLength1 = getDistanceByAndroid(polygon[i][0],polygon[i][1],polygon[i][0],point[1]);
double tempLength2 = getDistanceByAndroid(polygon[next][0],polygon[next][1],polygon[next][0],point[1]);
/**参考距离*/
double tempLength4 = getDistanceByAndroid(point[0],point[1],polygon[next][0],point[1]);
double tempLength5 = getDistanceByAndroid(polygon[i][0],point[1],polygon[next][0],point[1]);
if((tempLength4/tempLength5)>(tempLength2/(tempLength1+tempLength2))){
/**经过一个交点*/
System.out.println("~~夹点判断1:"+i);
node++;
}
}else if(polygon[i][0]>point[0]&& polygon[next][0]<point[0]){
/***/
double tempLength1 = getDistanceByAndroid(polygon[i][0],polygon[i][1],polygon[i][0],point[1]);
double tempLength2 = getDistanceByAndroid(polygon[next][0],polygon[next][1],polygon[next][0],point[1]);
/**参考距离*/
double tempLength3 = getDistanceByAndroid(point[0],point[1],polygon[i][0],point[1]);
double tempLength5 = getDistanceByAndroid(polygon[i][0],point[1],polygon[next][0],point[1]);
if((tempLength3/tempLength5)>(tempLength1/(tempLength1+tempLength2))){
/**经过一个交点*/
System.out.println("~~夹点判断2:"+i);
node++;
}
}
/**
* polygon[i][0]<point[0]&& polygon[next][0]<point[0]
相交于射线的反向延长线上(最后一种可能，基于性能不做判断)
* */
}
}
}
System.out.println("~~node:"+node);
if(node%2==0) return false;
else return true;
}
3.计算点到多边形的距离
直接计算我没有好的方法，我计算的是配合上面判断点是否在多边形内部，然后计算点到各个多边形边的最短距离
- 1.遍历多边形，得出每条边
- 2.已知给定点到线段端点的长度，以及线段长度，求出点到线段的距离
算法1，已知三边长，求高：
/** * 计算三角上a边上的高
*@param a a边长
*@param b b边长
*@param c c边长
* */
public static double getTrigCh(double b,double c,double a){
double p1 = a+c+b;
double p2 = a+c-b;
double p3 = b+a-c;
double p4 = b-a+c;
double r1 = p1*p2*p3*p4;
double r2 = Math.sqrt(r1);
return r2/(2*a);
}
算法2，计算点到多边形的距离：
/**计算点到多边形之间的距离*/
public static double calculateDistance(double[] point,double[][] polygon){
double[] distances = new double[polygon.length];
double[] sideLengths = new double[polygon.length];
for (int i = 0; i < polygon.length; i++) {
int next = i+1;
if (i==polygon.length-1)
next = 0;
distances[i] = getDistance(point[0],point[1],polygon[i][0],polygon[i][1]);
sideLengths[i] = getDistance(polygon[i][0],polygon[i][1],polygon[next][0],polygon[next][1]);
System.out.println("~distances["+i+"]:"+distances[i]);
System.out.println("~sideLengths["+i+"]:"+sideLengths[i]);
}
double minDistance = Double.MAX_VALUE;
for (int i = 0; i < polygon.length; i++) {
int next = i+1;
if (i==polygon.length-1)
next = 0;
double tampDistance;
if(distances[i]>distances[next]&&distances[i]>sideLengths[i]){
tampDistance = distances[next];
}else if(distances[next]>distances[i]&&distances[next]>sideLengths[i]){
tampDistance = distances[i];
}else{
tampDistance = getTrigCh(distances[i],distances[next],sideLengths[i]);
}
System.out.println("~tampDistance:"+tampDistance);
minDistance = (minDistance>tampDistance)?tampDistance:minDistance;
}
return minDistance;
}