我有一个包含 4 个字段的英国邮政编码数据库表:

postcode,
east,
north,
pqi 

示例值:ST1 6BQ, 388605, 349057,10

主键是 postcodes 和 east 和 north 都是 int(11)` 字段。

有没有人有一个高效的 MySQL 查询，可以返回给定帖子的最近 5 个邮政编码？

我见过很多使用long 和lat 的示例，但没有使用northing 和easting。

最佳答案

将您的北距/东距转换为纬度/经度，请记住操作系统网格基于 OSGB36 而不是 WGS84。

我使用以下类:

<?php

namespace OSGB36;

use \Geodetic\Datum;

class Converter
{
    private $_osRef;
    private $_fromDatum;
    private $_toDatum;

    public function __construct()
    {
        $this->_osRef = new OSRef();
        $this->_fromDatum = new Datum(Datum::OSGB36);
        $this->_toDatum = new Datum(Datum::WGS84);
    }

    /**
     * Converts easting/northing into lat/long
     * 
     * @param integer $eastings
     * @param integer $northings
     * @return \Geodetic\LatLong
     */
    public function calculateLatLong($eastings, $northings)
    {
        $this->_osRef->setNorthings($northings)
                     ->setEastings($eastings);
        $OSGB36LatLong = $this->_osRef->toLatLong(
            $this->_fromDatum->getReferenceEllipsoid()
        );

        $ecef = $OSGB36LatLong->toECEF($this->_fromDatum);
        $ecef->toWGS84($this->_fromDatum);
        $WGS84LatLong = $ecef->toLatLong($this->_toDatum);

        return $WGS84LatLong;
    }
}

和

<?php

namespace OSGB36;

use \Geodetic\LatLong\CoordinateValues;

class OSRef
{
    private $_northings;
    private $_eastings;

    public function __construct($northings = NULL, $eastings = NULL)
    {
        $this->_northings = $northings;
        $this->_eastings  = $eastings;
    }

    public function setNorthings($northings)
    {
        $this->_northings = $northings;

        return $this;
    }

    public function setEastings($eastings)
    {
        $this->_eastings = $eastings;

        return $this;
    }


    private function _sinSquared($x) {
        return sin($x) * sin($x);
    }

    private function _tanSquared($x) {
        return tan($x) * tan($x);
    }

    private function _secant($x) {
        return 1.0 / cos($x);
    }

    private function _cosecant($x) {
        return 1.0 / sin($x);
    }

    private function _cotangent($x) {
        return 1.0 / tan($x);
    }

    public function toLatLong(\Geodetic\ReferenceEllipsoid $airy1830) {
        $OSGB_F0       = 0.9996012717;    //  Central Meridan Scale factor
        $N0            = -100000.0;       //  True origin Northing
        $E0            = 400000.0;        //  True origin Easting
        $phi0          = deg2rad(49.0);   //  True origin Latitude
        $lambda0       = deg2rad(-2.0);   //  True origin Longitude
        $semiMajorAxis = $airy1830->getSemiMajorAxis();
        $semiMinorAxis = $airy1830->getSemiMinorAxis();
        $eSquared      = $airy1830->getFirstEccentricitySquared();
        $easting       = $this->_eastings - $E0;
        $northing      = $this->_northings - $N0;
        $n             = ($semiMajorAxis - $semiMinorAxis) / ($semiMajorAxis + $semiMinorAxis);
        $M             = 0.0;
        $phiPrime      = ($northing / ($semiMajorAxis * $OSGB_F0)) + $phi0;

        do {
            $M = ($semiMinorAxis * $OSGB_F0) *
                (((1 + $n + ((5.0 / 4.0) * $n * $n) + ((5.0 / 4.0) * $n * $n * $n)) *
                ($phiPrime - $phi0)) -
                (((3 * $n) + (3 * $n * $n) + ((21.0 / 8.0) * $n * $n * $n)) *
                sin($phiPrime - $phi0) *
                cos($phiPrime + $phi0)) +
                ((((15.0 / 8.0) * $n * $n) + ((15.0 / 8.0) * $n * $n * $n)) *
                sin(2.0 * ($phiPrime - $phi0)) *
                cos(2.0 * ($phiPrime + $phi0))) -
                (((35.0 / 24.0) * $n * $n * $n) *
                sin(3.0 * ($phiPrime - $phi0)) *
                cos(3.0 * ($phiPrime + $phi0))));
            $phiPrime += ($northing - $M) / ($semiMajorAxis * $OSGB_F0);
        } while (($northing - $M) >= 0.001);

        $v = $semiMajorAxis * $OSGB_F0 * pow(1.0 - $eSquared * $this->_sinSquared($phiPrime), -0.5);
        $rho = $semiMajorAxis * $OSGB_F0 * (1.0 - $eSquared) *
            pow(1.0 - $eSquared * $this->_sinSquared($phiPrime), -1.5);
        $etaSquared = ($v / $rho) - 1.0;
        $VII = tan($phiPrime) / (2 * $rho * $v);
        $VIII = (tan($phiPrime) / (24.0 * $rho * pow($v, 3.0))) *
            (5.0 + (3.0 * $this->_tanSquared($phiPrime)) + $etaSquared - (9.0 * $this->_tanSquared($phiPrime) * $etaSquared));
        $IX = (tan($phiPrime) / (720.0 * $rho * pow($v, 5.0))) *
            (61.0 + (90.0 * $this->_tanSquared($phiPrime)) + (45.0 * $this->_tanSquared($phiPrime) * $this->_tanSquared($phiPrime)));
        $X = $this->_secant($phiPrime) / $v;
        $XI = ($this->_secant($phiPrime) / (6.0 * $v * $v * $v)) * (($v / $rho) + (2 * $this->_tanSquared($phiPrime)));
        $XII =  ($this->_secant($phiPrime) / (120.0 * pow($v, 5.0))) *
            (5.0 + (28.0 * $this->_tanSquared($phiPrime)) + (24.0 * $this->_tanSquared($phiPrime) * $this->_tanSquared($phiPrime)));
        $XIIA = ($this->_secant($phiPrime) / (5040.0 * pow($v, 7.0))) *
            (61.0 + (662.0 * $this->_tanSquared($phiPrime)) + (1320.0 * $this->_tanSquared($phiPrime) * $this->_tanSquared($phiPrime)) +
            (720.0 * $this->_tanSquared($phiPrime) * $this->_tanSquared($phiPrime) * $this->_tanSquared($phiPrime)));
        $phi = $phiPrime - ($VII * pow($easting, 2.0)) + ($VIII * pow($easting, 4.0)) - ($IX * pow($easting, 6.0));
        $lambda = $lambda0 +
            ($X * $easting) - ($XI * pow($easting, 3.0)) + ($XII * pow($easting, 5.0)) - ($XIIA * pow($easting, 7.0));

        $latLongCoordinates = new CoordinateValues(
            $phi,
            $lambda,
            \Geodetic\Angle::RADIANS,
            0.0,
            \Geodetic\Distance::METRES
        );
        return new \Geodetic\LatLong($latLongCoordinates);
    }

    function toGridRef()
    {
        $hundredkmE = floor($this->_eastings / 100000);
        $hundredkmN = floor($this->_northings / 100000);
        $firstLetter = "";
        if ($hundredkmN < 5) {
            if ($hundredkmE < 5) {
                $firstLetter = "S";
            } else {
                $firstLetter = "T";
            }
        } else if ($hundredkmN < 10) {
            if ($hundredkmE < 5) {
                $firstLetter = "N";
            } else {
                $firstLetter = "O";
            }
        } else {
            $firstLetter = "H";
        }

        $secondLetter = "";
        $index = 65 + ((4 - ($hundredkmN % 5)) * 5) + ($hundredkmE % 5);
        $ti = $index;
        if ($index >= 73)
        {
            $index++;
        }
        $secondLetter = chr($index);

        $e = round(($this->_eastings - (100000 * $hundredkmE)) / 100);
        $n = round(($this->_northings - (100000 * $hundredkmN)) / 100);

        return sprintf("%s%s%03d%03d", $firstLetter, $secondLetter, $e, $n);
    }

    public static function createOSRefFromGridRef($gridRef)
    {
        $char1 = substr($gridRef, 0, 1);
        $char2 = substr($gridRef, 1, 1);
        $east  = substr($gridRef, 2, 3) * 100;
        $north = substr($gridRef, 5, 3) * 100;
        if ($char1 == 'H') {
            $north += 1000000;
        } else if ($char1 == 'N') {
            $north += 500000;
        } else if ($char1 == 'O') {
            $north += 500000;
            $east  += 500000;
        } else if ($char1 == 'T') {
            $east += 500000;
        }
        $char2ord = ord($char2);
        if ($char2ord > 73)
        {
            $char2ord--; // Adjust for no I
        }
        $nx = (($char2ord - 65) % 5) * 100000;
        $ny = (4 - floor(($char2ord - 65) / 5)) * 100000;
        return new OSRef($north + $ny, $east + $nx);
    }
}

结合我的Geodetic library .

通过打开整个 CodePoint 给我一个邮政编码详细信息数据库(包括用于 OpenStreetMap 的 WGS84 纬度/经度坐标)的数据库非常快。

关于php - 使用军械测量东向和北向系统查找最接近给定邮政编码的邮政编码，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/16667121/