I have a database table of UK postcodes that has 4 fields:
postcode,
east,
north,
pqi
Example values: ST1 6BQ, 388605, 349057,10
The primary key is postcodes and east and north are bothint(11)` fields.
Does anybody have an efficient MySQL query that will return the nearest 5 postcodes to a given post?
I have seen many examples using long and lat, but not northing and easting.
分享 Convert your northings/eastings to lat/long, remembering that The OS grid is based on OSGB36 rather than WGS84.
I use the following class:
<?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;
}
}
and
<?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);
}
}
combined with my Geodetic library.
It's pretty quick chundering through the entirety of CodePoint open to give me a database of postcode details including WGS84 lat/long coordinates for use with OpenStreetMap.
分享
