[PHP] 计算两个坐标点之间的现实距离 (m)

29 阅读4分钟

此方法亦可实现最简单的圆型电子围栏

另外附上多边形判断

<?php

namespace utils;

class GpsUtils 
{

    private const a = 6378137; // 长半径a=6378137
    private const b = 6356752.314245; // 短半径b=6356752.314245
    private const f = 1 / 298.257223563; // 扁率f=1/298.257223563

    /**
     *  计算两个坐标之间间距离 Vincenty公式
     *  扩展思路: 当前坐标距离上一点距离, 以及圆心距离较远时, 再去判断是否超出围栏 (其实距离圆心超出围栏半径, 即可判断为超出围栏)
     * @return float
     */
    public static function getDistance($ponits1, $ponits2) 
    {

        $ponits1Arr = explode(',', $ponits1);
        $ponits2Arr = explode(',', $ponits2);
        $lat_one = $ponits1Arr[1];
        $lon_one = $ponits1Arr[0];
        $lat_two = $ponits2Arr[1];
        $lon_two = $ponits2Arr[0];

        $L = self::toRadians($lon_one - $lon_two);
        $U1 = atan((1 - self::f) * tan(self::toRadians($lat_one)));
        $U2 = atan((1 - self::f) * tan(self::toRadians($lat_two)));
        $sinU1 = sin($U1);
        $cosU1 = cos($U1);
        $sinU2 = sin($U2);
        $cosU2 = cos($U2);
        $lambda = $L;
        $lambdaP = pi();
        $cosSqAlpha = 0;
        $sinSigma = 0;
        $cos2SigmaM = 0;
        $cosSigma = 0;
        $sigma = 0;
        $circleCount = 40;

        while (abs($lambda - $lambdaP) > 1e-12 && --$circleCount > 0) {
            $sinLambda = sin($lambda);
            $cosLambda = cos($lambda);
            $sinSigma = sqrt(($cosU2 * $sinLambda) * ($cosU2 * $sinLambda) +
                    ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda));
            if ($sinSigma == 0) {
                return 0;  // co-incident points
            }
            $cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
            $sigma = atan2($sinSigma, $cosSigma);
            $alpha = asin($cosU1 * $cosU2 * $sinLambda / $sinSigma);
            $cosSqAlpha = cos($alpha) * cos($alpha);
            $cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;
            $C = self::f / 16 * $cosSqAlpha * (4 + self::f * (4 - 3 * $cosSqAlpha));
            $lambdaP = $lambda;
            $lambda = $L + (1 - $C) * self::f * sin($alpha) *
                    ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM)));
        }

        if ($circleCount == 0) {
            return NAN;  // formula failed to converge
        }

        $uSq = $cosSqAlpha * (self::a * self::a - self::b * self::b) / (self::b * self::b);
        $A = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
        $B = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
        $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) -
                $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma) * (-3 + 4 * $cos2SigmaM * $cos2SigmaM)));

        $result = self::b * $A * ($sigma - $deltaSigma);
        $distance = round($result, 5); // 保留5位小数
        return $distance;
    }

    private static function toRadians($angle) {
        return $angle * pi() / 180;
    }

    /**
     * 判断点是否在多边形内
     * @return bool
     */
    public static function containsPoint($vertices, $ponits) 
    {
        $ponitsArr = explode(',', $ponits);
        $lat = $ponitsArr[1];
        $lng = $ponitsArr[0];


        $crossings = 0;
        $count = count($vertices);

        for ($i = 0; $i < $count; $i++) {
            $a = $vertices[$i];
            $b = $vertices[($i + 1) % $count];

            if ((($a[1] <= $lng) && ($lng < $b[1])) || (($b[1] <= $lng) && ($lng < $a[1]))) {
                $vt = ($lng - $a[1]) / ($b[1] - $a[1]);
                if ($lat < $a[0] + $vt * ($b[0] - $a[0])) {
                    $crossings++;
                }
            }
        }
        // 奇数次交点:点在多边形内部;偶数次交点:点在多边形外部
        return ($crossings % 2 != 0);
    }

    public function java()
    {
        $eof = <<<EOF

       
EOF;
    }
    
}

java

package com.util.geo;

import cn.hutool.core.collection.CollectionUtil;
import cn.hutool.core.util.StrUtil;
import com.util.geo.core.*;

import java.util.Arrays;
import java.util.List;
import java.util.ArrayList;


public class GeoUtils {



    /**
     * 经纬度逗号隔开 重新处理
     *
     * @param location 经度,纬度
     * @return 。
     */
    public static Point parseLocation(String location) {
        List<String> points = StrUtil.split(location, ",");
        return new Point(Double.parseDouble(points.get(0)), Double.parseDouble(points.get(1)));
    }
//    public static Point parseLocation(String location) {
//        List<String> points = Arrays.asList(StrUtil.split(location, ","));
//        return new Point(Double.parseDouble(points.get(0)), Double.parseDouble(points.get(1)));
//    }




    /**
     * 计算顺序标点的距离,必须是顺序标表否则计算不准确
     *
     * @param points 标点
     * @return 。
     */
    public static double calculateGeoDistanceSphere(List<Point> points) {
        //顺序计算总距离
        if (CollectionUtil.isEmpty(points) || points.size() == 1) {
            return 0.00;
        }
        Point tempPoint = points.get(0);
        double distance = 0.00;
        for (int i = points.size() - 1; i >= 1; i--) {
            Point targetPoint = points.get(i);
            distance += calculateGeoDistanceSphere(tempPoint, targetPoint);
            tempPoint = targetPoint;
        }
        return distance;
    }

    /**
     * 计算经纬度距离
     *
     * @param point       经纬度
     * @param targetPoint 目标经纬度
     * @return 距离 单位M
     */
    public static double calculateGeoDistanceSphere(Point point, Point targetPoint) {
        GlobalCoordinates source = new GlobalCoordinates(point.getLat(), point.getLng());
        GlobalCoordinates target = new GlobalCoordinates(targetPoint.getLat(), targetPoint.getLng());
        GeodeticCurve getS = new GeodeticCalculator().calculateGeodeticCurve(Ellipsoid.Sphere, source, target);
        return getS.getEllipsoidalDistance();
    }

    public static double calculateGeoDistanceSphere2(Point point, Point targetPoint) {
        GlobalCoordinates source = new GlobalCoordinates(point.getLat(), point.getLng());
        GlobalCoordinates target = new GlobalCoordinates(targetPoint.getLat(), targetPoint.getLng());
        GeodeticCurve getS = new GeodeticCalculator().calculateGeodeticCurve(Ellipsoid.Sphere, source, target);
        return getS.getEllipsoidalDistance();
    }

    /**
     * 计算是否在圆内
     *
     * @param radius      半径(单位/米)
     * @param centerPoint 圆心坐标
     * @param targetPoint 判断点坐标
     * @return boolean true:在圆内,false:在圆外
     */
    public static boolean isInCircle(double radius, Point centerPoint, Point targetPoint) {
        //先求圆心到目标的距离
        double distance = calculateGeoDistanceSphere(centerPoint, targetPoint);
        return distance <= radius;
    }


    /**
     * 判断点是否在多边形内部
     * <p>
     * 将多边形顶点坐标按照顺序存储在数组中;
     * 创建一个方法,参数为需要检测的点的坐标(x,y)和多边形的顶点数组verticesvertices;
     * 定义一个变量counter,用于记录交点数量初始为0;
     * 循环遍历vertices数组中的每一条边:假设当前遍历到的是第i条边,它连接(vertices[i].x, vertices[i].y)和(vertices[i+1].x, vertices[i+1].y)两个点;
     * 计算射线与当前边(vertices[i],vertices[i+1])是否相交,计算方法如下:
     * 如果射线位置在多边形某个凸角上,忽略该点;
     * 如果射线位置恰好在当前边的下端点(vertices[i+1].y==point.y),忽略该点;
     * 如果射线位置在当前边的上端点(vertices[i].y==point.y),将交点数量加1;
     * 如果当前边不水平,则计算射线y坐标等于此时两侧端点高度投影相交点横坐标,若此值小于x,则交点数量加1;
     * 根据交点数量的奇偶性来判断该点是否在围栏内。如果交点数量为奇数,则该点在围栏内,否则在围栏外。
     * 如果该点在围栏内,返回true;否则返回false。
     * 需要注意的是,在遍历多边形的边的时候,vertices[i+1]指的是下一个顶点,如果当前边已经是最后一条边,那么vertices[i+1]应该是多边形的第一个顶点。
     */
    public boolean containsPoint(List<Double[]> vertices, double lat, double lng) {
        int crossings = 0;
        for (int i = 0; i < vertices.size(); i++) {
            Double[] a = vertices.get(i);
            Double[] b = vertices.get((i + 1) % vertices.size());
            if (((a[1] <= lng) && (lng < b[1])) ||
                    ((b[1] <= lng) && (lng < a[1]))) {
                double vt = (lng - a[1]) / (b[1] - a[1]);
                if (lat < a[0] + vt * (b[0] - a[0])) {
                    crossings++;
                }
            }
        }
        // 奇数次交点:点在多边形内部;偶数次交点:点在多边形外部
        return (crossings % 2 != 0);
    }

    public boolean containsPoint(List<Point> vertices, Point point) {
        int crossings = 0;
        for (int i = 0; i < vertices.size(); i++) {
            Point a = vertices.get(i);
            Point b = vertices.get((i + 1) % vertices.size());

            if (((a.getLng() <= point.getLng()) && (point.getLng() < b.getLng())) ||
                    ((b.getLng() <= point.getLng()) && (point.getLng() < a.getLng()))) {
                double vt = (point.getLng() - a.getLng()) / (b.getLng() - a.getLng());
                if (point.getLat() < a.getLat() + vt * (b.getLat() - a.getLat())) {
                    crossings++;
                }
            }
        }
        // 奇数次交点:点在多边形内部;偶数次交点:点在多边形外部
        return (crossings % 2 != 0);
    }


}
package com.util.geo;

import com.util.geo.GeoUtils;
import com.util.geo.core.Point;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;

import java.util.ArrayList;
import java.util.List;

public class GeoUtilsTest {

    public static void main(String[] args) {
        GeoUtilsTest geoUtilsTest = new GeoUtilsTest();
        geoUtilsTest.testContainsPoint();
    }

    @Test
    public void testContainsPoint() {
        List<Double[]> vertices = new ArrayList<>();
        vertices.add(new Double[]{116.282962, 39.979764});
        vertices.add(new Double[]{116.474983, 39.978437});
        vertices.add(new Double[]{116.484757, 39.846056});
        vertices.add(new Double[]{116.295610, 39.839851});

        double testLat = 39.929322;
        double testLng = 116.395645;

        boolean isInside = new GeoUtils().containsPoint(vertices, testLat, testLng);

        System.out.println("Point is inside polygon: " + (isInside ? "Yes" : "No"));
    }
}