<?php
/*******************************************************
Grid Reference Utilities
Version 3.0 - Written by Mark Wilton-Jones 1-10/4/2010
Updated 6/5/2010 to allow methods to return text
strings and add support for ITM, UTM and UPS
Updated 9/5/2010 to add dd_format
Updated 16/9/2020 to correct Helmert height handling
Updated 19/10/2020 to correct DD detection in dms_to_dd
Updated 23/12/2020 to add support for shift grids,
geoids, polynomial transformation, geodesics, trucated
grid references, floating point grid references, high
precision GPS coordinates, and more input formats
********************************************************
Please see http://www.howtocreate.co.uk/php/ for details
Please see http://www.howtocreate.co.uk/php/gridref.php for demos and instructions
Please see http://www.howtocreate.co.uk/php/gridrefapi.php for API documentation
Please see http://www.howtocreate.co.uk/jslibs/termsOfUse.html for terms and conditions of use
Provides functions to convert between different NGR formats and latitude/longitude formats.
This script contains the coefficients for OSi/LPS Polynomial Transformation. Coefficients are
available under BSD License: (c) Copyright and database rights Ordnance Survey Limited 2016,
(c) Crown copyright and and database rights Land & Property Services 2016 and/or (c) Ordnance
Survey Ireland, 2016. All rights reserved.
_______________________________________________________________________________________________*/
//use a class to keep all of the method/property clutter out of the global scope
class Grid_Ref_Utils {
//class instantiation forcing use as a singleton - state can also be stored if needed (or can use separate instances)
//could be done as an abstract class with static methods, but would still need instantiation to prepare UK_grid_numbers
//(as methods/properties do not exist during initial class constructor parsing), and would restrict future flexibility
private static $only_instance;
public static function toolbox() {
if( !isset( self::$only_instance ) ) {
self::$only_instance = new Grid_Ref_Utils();
}
return self::$only_instance;
}
private function __construct() {
$this->UK_grid_numbers = $this->grid_to_hashmap($this->UK_grid_squares);
$this->TEXT_EUROPEAN = html_entity_decode( '°', ENT_QUOTES, 'ISO-8859-1' );
$this->TEXT_UNICODE = html_entity_decode( '°', ENT_QUOTES, 'UTF-8' );
$this->TEXT_ASIAN = chr(129).chr(139); //x818b, Shift_JIS html_entity_decode fails
$this->TEXT_DEFAULT = html_entity_decode( '°', ENT_QUOTES );
if( !defined('PHP_FLOAT_EPSILON') ) {
//used by geodesy functions to test for accuracy limits
$epsilon = 1.0;
do {
$epsilon /= 2.0;
} while( 1.0 + ( $epsilon / 2.0 ) != 1.0 );
define( 'PHP_FLOAT_EPSILON', $epsilon );
}
}
//define return types
public $DATA_ARRAY = false;
public $HTML = true;
//will be set up when the object is constructed
public $TEXT_EUROPEAN;
public $TEXT_UNICODE;
public $TEXT_ASIAN;
public $TEXT_DEFAULT;
//common return values
private function generic_xy_coords( $return_type, $precise, $easting, $northing, $height = null, $additional = null ) {
//avoid stupid -0 by adding 0
if( $return_type ) {
if( is_nan($easting) || is_nan($northing) ) {
return 'OUT_OF_GRID';
} elseif( $precise ) {
return sprintf( '%F', $easting ) . ',' . sprintf( '%F', $northing ) . ( is_numeric($height) ? ( ',' . sprintf( '%F', $height ) ) : '' );
} else {
return ( round($easting) + 0 ) . ',' . ( round($northing) + 0 ) . ( is_numeric($height) ? ( ',' . ( round($height) + 0 ) ) : '' );
}
} else {
$returnarray = array( $easting + 0, $northing + 0 );
if( is_numeric($height) ) {
$returnarray[] = $height + 0;
}
if( $additional !== null ) {
$returnarray[] = $additional;
}
return $returnarray;
}
}
private function generic_lat_long_coords( $return_type, $latitude, $longitude, $height = null ) {
//force the longitude between -180 and 180
$longitude = $this->wrap_around_longitude($longitude);
if( $return_type ) {
//sprintf to produce simple numbers instead of scientific notation
//10 decimal places DD ~= 8 decimal places DM ~= 6 decimal places DMS
$deg = is_string($return_type) ? $return_type : '°';
$longitude = sprintf( '%.10F', $longitude );
if( $longitude == '-180.0000000000' ) {
//rounding can make -179.9...9 end up as -180
$longitude = '180.0000000000';
}
return sprintf( '%.10F', $latitude ) . $deg . ', ' . $longitude . $deg . ( is_numeric($height) ? ( ', ' . sprintf( '%F', $height ) ) : '' );
} else {
//avoid stupid -0 by adding 0
$temparray = array( $latitude + 0, $longitude + 0 );
if( is_numeric($height) ) {
$temparray[] = $height + 0;
}
return $temparray;
}
}
//common inputs
function expand_params( $ar1, $from, $to, $ar2 ) {
//PHP doesn't support overloaded functions/methods
//this function makes it possible to extract an initial array, returning parameters in a standard order
$param = 0;
$handback = array();
if( is_array($ar1) ) {
foreach( $ar1 as $onecell ) {
if( $param == $to ) {
//ignore extra cells
break;
}
$handback[$param] = $onecell;
$param++;
}
//this forces attempts to use an array of one item to then continue without the next param(s)
while( $param < $from ) {
//prefill with null so that PHP doesn't get upset about missing cells
$handback[$param] = null;
$param++;
}
} else {
//not using an array
$handback[0] = $ar1;
$param = 1;
}
foreach( $ar2 as $onecell ) {
if( $param < $to && is_array($onecell) ) {
//the first parameter that is an array must at least be at the "to" index
//this allows "height" to be put either inside or outside an array, and still be optional
//additional parameters in between are preserved and the calling function will continue with them
while( $param < $to ) {
$handback[$param] = null;
$param++;
}
$param = $to;
}
$handback[$param] = $onecell;
$param++;
}
return $handback;
}
//flooring to a number of significate figures
private function floor_precision( $num, $precision ) {
//regular floor does not support the precision param
//usual method floor( $number * pow(10, $precision) ) / pow(10, $precision) gets floating point errors with ( 9.7, 2 )
//string manipulation fails with 1E-15
//this function is immune to both problems
$rounded = round( $num, $precision );
if( $rounded > $num ) {
//pow works during testing for negative integer exponents, but perhaps it will fail on some systems, since JS does
//$rounded -= pow( 10, -1 * $precision );
//this approach always works
$rounded -= '1E' . ( -1 * $precision );
}
return $rounded + 0;
}
//longitude must be within -180<num<=180
private function wrap_around_longitude( $num, $negative_half = false ) {
if( !is_numeric($num) ) {
return $num;
}
if( abs($num) > 1440 ) {
//this can end up with floating point errors much more easily than simple + or -
//use it only in cases where the error is very large, to avoid the loops taking too long
$num -= 360 * floor( $num / 360 );
}
while( ( !$negative_half && $num > 180 ) || ( $negative_half && $num >= 180 ) ) {
$newnum = $num - 360;
if( $newnum == $num ) {
//avoid infinite loops with numbers high enough to need scientific notation
return $num;
}
$num = $newnum;
}
while( ( !$negative_half && $num <= -180 ) || ( $negative_half && $num < -180 ) ) {
$newnum = $num + 360;
if( $newnum == $num ) {
//avoid infinite loops with numbers high enough to need scientific notation
return $num;
}
$num = $newnum;
}
return $num;
}
//azimuth must be within 0<=num<360
private function wrap_around_azimuth($num) {
if( !is_numeric($num) ) {
return $num;
}
if( abs($num) > 1440 ) {
//this can end up with floating point errors much more easily than simple + or -
//use it only in cases where the error is very large, to avoid the loops taking too long
$num -= 360 * floor( $num / 360 );
}
while( $num >= 360 ) {
$newnum = $num - 360;
if( $newnum == $num ) {
//avoid infinite loops with numbers high enough to need scientific notation
return $num;
}
$num = $newnum;
}
while( $num < 0 ) {
$newnum = $num + 360;
if( $newnum == $num ) {
//avoid infinite loops with numbers high enough to need scientific notation
return $num;
}
$num = $newnum;
}
return $num;
}
//character grids used by map systems
private function grid_to_hashmap($grid_2d_array) {
//make a hashmap of arrays giving the x,y values for grid letters
$hashmap = array();
foreach( $grid_2d_array as $grid_row_index => $grid_row_array ) {
foreach( $grid_row_array as $grid_col_index => $grid_letter ) {
$hashmap[$grid_letter] = array($grid_col_index,$grid_row_index);
}
}
return $hashmap;
}
private $UK_grid_squares = array(
//the order of grid square letters in the UK NGR system - note that in the array they start from the bottom, like grid references
//there is no I in team
array('V','W','X','Y','Z'),
array('Q','R','S','T','U'),
array('L','M','N','O','P'),
array('F','G','H','J','K'),
array('A','B','C','D','E')
);
private $UK_grid_numbers; //will be set up when the object is constructed
//conversions between 12345,67890 and SS234789 grid reference formats
public function get_UK_grid_ref( $orig_x, $orig_y = false, $figures = false, $return_type = false, $deny_bad_reference = false, $use_rounding = false, $is_irish = false ) {
//no need to shuffle the $is_irish parameter over, since it will always be at the end
list( $orig_x, $orig_y, $figures, $return_type, $deny_bad_reference, $use_rounding ) =
$this->expand_params( $orig_x, 2, 2, array( $orig_y, $figures, $return_type, $deny_bad_reference, $use_rounding ) );
if( is_bool($figures) || $figures === null ) {
$figures = 4;
} else {
//enforce integer 1-25
$figures = max( min( floor( $figures ), 25 ), 1 );
}
//prepare factors used for enforcing a number of grid ref figures
$insig_figures = $figures - 5;
if( $use_rounding ) {
//round off unwanted detail
$x = round( $orig_x, $insig_figures );
$y = round( $orig_y, $insig_figures );
} else {
$x = $this->floor_precision( $orig_x, $insig_figures );
$y = $this->floor_precision( $orig_y, $insig_figures );
}
$errorletters = 'OUT_OF_GRID';
if( $is_irish ) {
//the Irish grid system uses the same letter layout as the UK system, but it only has one letter, with origin at V
$arY = floor( $y / 100000 );
$arX = floor( $x / 100000 );
if( $arX < 0 || $arX > 4 || $arY < 0 || $arY > 4 ) {
//out of grid
if( $deny_bad_reference ) {
return false;
}
$letters = $errorletters;
} else {
$letters = $this->UK_grid_squares[ $arY ][ $arX ];
}
} else {
//origin is at SV, not VV - offset it to VV instead
$x += 1000000;
$y += 500000;
$ar1Y = floor( $y / 500000 );
$ar1X = floor( $x / 500000 );
$ar2Y = floor( ( $y % 500000 ) / 100000 );
$ar2X = floor( ( $x % 500000 ) / 100000 );
if( is_nan($ar1X) || is_nan($ar1Y) || is_nan($ar2X) || is_nan($ar2Y) || $ar1X < 0 || $ar1X > 4 || $ar1Y < 0 || $ar1Y > 4 ) {
//out of grid - don't need to test $ar2Y and $ar2X since they cannot be more than 4, and ar1_ will also be less than 0 if ar2_ is
if( $deny_bad_reference ) {
return false;
}
$letters = $errorletters;
} else {
//first grid letter is for the 500km x 500km squares
$letters = $this->UK_grid_squares[ $ar1Y ][ $ar1X ];
//second grid letter is for the 100km x 100km squares
$letters .= $this->UK_grid_squares[ $ar2Y ][ $ar2X ];
}
}
//% casts to integer, so need to use floor instead
//$x %= 100000;
//$y %= 100000;
//floating point errors can reappear here if using numbers after the decimal point
//this approach also makes negative become positive remainder, but that is wanted anyway
$x -= floor( $x / 100000 ) * 100000;
$y -= floor( $y / 100000 ) * 100000;
//avoid stupid -0 by adding 0
if( $figures <= 5 ) {
$figure_factor = pow( 10, 5 - $figures );
$x = str_pad( ( ( $x % 100000 ) / $figure_factor ) + 0, $figures, '0', STR_PAD_LEFT );
$y = str_pad( ( ( $y % 100000 ) / $figure_factor ) + 0, $figures, '0', STR_PAD_LEFT );
} else {
//pad only the left side of the decimal point
//use sprintf to remove floating point errors introduced by -= and avoid stupid -0
$x = explode( '.', sprintf( '%.'.$insig_figures.'F', $x ) );
$x[0] = str_pad( $x[0], 5, '0', STR_PAD_LEFT );
$x = implode( '.', $x );
$y = explode( '.', sprintf( '%.'.$insig_figures.'F', $y ) );
$y[0] = str_pad( $y[0], 5, '0', STR_PAD_LEFT );
$y = implode( '.', $y );
}
if( $return_type && is_string($return_type) ) {
return $letters . ' ' . $x . ' ' . $y;
} elseif( $return_type ) {
return '<var class="grid">' . $letters . '</var><var>' . $x . '</var><var>' . $y . '</var>';
} else {
return array( $letters, $x, $y );
}
}
public function get_UK_grid_nums( $letters, $x = false, $y = false, $return_type = false, $deny_bad_reference = false, $use_rounding = false, $precise = false, $is_irish = false ) {
//no need to shuffle the $is_irish parameter over, since it will always be at the end
list( $letters, $x, $y, $return_type, $deny_bad_reference, $use_rounding, $precise ) =
$this->expand_params( $letters, 3, 3, array( $x, $y, $return_type, $deny_bad_reference, $use_rounding, $precise ) );
if( !is_string( $x ) || !is_string( $y ) ) {
//a single string 'X[Y]12345678' or 'X[Y] 1234 5678'
//manually move params
$precise = $deny_bad_reference;
$use_rounding = $return_type;
$deny_bad_reference = $y;
$return_type = $x;
//split into parts
//captured whitespace hack makes sure it fails reliably (UK/Irish grid refs must not match each other), while giving consistent matches length
if( !preg_match( "/^\s*([A-HJ-Z])".($is_irish?"(\s*)":"([A-HJ-Z])\s*")."([\d\.]+)\s*([\d\.]*)\s*$/", strtoupper($letters), $letters ) || ( !$letters[4] && strlen($letters[3]) < 2 ) || ( !$letters[4] && strpos( $letters[3], '.' ) !== false ) ) {
//invalid format
if( $deny_bad_reference ) {
return false;
}
//assume 0,0
$letters = array( '', $is_irish ? 'V' : 'S', 'V', '0', '0' );
$use_rounding = true;
}
if( !$letters[4] ) {
//a single string 'X[Y]12345678', break the numbers in half
$halflen = ceil( strlen($letters[3]) / 2 );
$letters[4] = substr( $letters[3], $halflen );
$letters[3] = substr( $letters[3], 0, $halflen );
}
$x = $letters[3];
$y = $letters[4];
} else {
//st becomes ['','S','T']
if( strlen($letters) != $is_irish ? 1 : 2 ) {
if( $deny_bad_reference ) {
return false;
}
}
$letters = str_split(strtoupper($letters));
array_unshift($letters,'');
}
$xdot = strpos( $x, '.' );
$ydot = strpos( $y, '.' );
if( is_numeric($x) ) {
if( !$use_rounding ) {
if( ( strlen($x) == 5 || $ydot !== false ) && $xdot === false ) {
$x .= '.5';
} else {
$x .= '5';
}
}
} else {
if( $deny_bad_reference ) {
return false;
}
$x = '0';
}
if( is_numeric($y) ) {
if( !$use_rounding ) {
if( ( strlen($y) == 5 || $xdot !== false ) && $ydot === false ) {
$y .= '.5';
} else {
$y .= '5';
}
}
} else {
if( $deny_bad_reference ) {
return false;
}
$y = '0';
}
//need 1m x 1m squares, but if either of them started with a dot, assume they are both in metres already
if( $ydot === false && strpos( $x, '.' ) === false ) {
$x = str_pad( $x, 5, '0' );
}
if( $xdot === false && strpos( $y, '.' ) === false ) {
$y = str_pad( $y, 5, '0' );
}
if( $is_irish ) {
if( isset( $this->UK_grid_numbers[$letters[1]] ) ) {
$x += $this->UK_grid_numbers[$letters[1]][0] * 100000;
$y += $this->UK_grid_numbers[$letters[1]][1] * 100000;
} elseif( $deny_bad_reference ) {
return false;
} else {
$x = $y = 0;
}
} else {
if( isset( $this->UK_grid_numbers[$letters[1]] ) && isset( $this->UK_grid_numbers[$letters[2]] ) ) {
//remove offset from VV to put origin back at SV
$x += ( $this->UK_grid_numbers[$letters[1]][0] * 500000 ) + ( $this->UK_grid_numbers[$letters[2]][0] * 100000 ) - 1000000;
$y += ( $this->UK_grid_numbers[$letters[1]][1] * 500000 ) + ( $this->UK_grid_numbers[$letters[2]][1] * 100000 ) - 500000;
} elseif( $deny_bad_reference ) {
return false;
} else {
$x = $y = 0;
}
}
if( $return_type && !$use_rounding ) {
//generic_xy_coords uses rounding
if( $precise ) {
$x = $this->floor_precision( $x, 6 );
$y = $this->floor_precision( $y, 6 );
} else {
$x = floor($x);
$y = floor($y);
}
}
return $this->generic_xy_coords( $return_type, $precise, $x, $y );
}
public function get_Irish_grid_ref( $x, $y = false, $figures = false, $return_type = false, $deny_bad_reference = false, $use_rounding = false ) {
return $this->get_UK_grid_ref( $x, $y, $figures, $return_type, $deny_bad_reference, $use_rounding, true );
}
public function get_Irish_grid_nums( $letters, $x = false, $y = false, $return_type = false, $deny_bad_reference = false, $use_rounding = false, $precise = false ) {
return $this->get_UK_grid_nums( $letters, $x, $y, $return_type, $deny_bad_reference, $use_rounding, $precise, true );
}
public function add_grid_units( $x, $y = false, $return_type = false, $precise = false ) {
list( $x, $y, $return_type, $precise ) = $this->expand_params( $x, 2, 2, array( $y, $return_type, $precise ) );
if( $return_type ) {
if( $precise ) {
$x = sprintf( '%F', $x );
$y = sprintf( '%F', $y );
} else {
//avoid stupid -0 by adding 0
$x = round($x) + 0;
$y = round($y) + 0;
}
return ( ( $x < 0 ) ? 'W ' : 'E ' ) . abs($x) . 'm, ' . ( ( $y < 0 ) ? 'S ' : 'N ' ) . abs($y) . 'm';
} else {
return array( abs($x), ( $x < 0 ) ? -1 : 1, abs($y), ( $y < 0 ) ? -1 : 1 );
}
}
public function parse_grid_nums( $coords, $return_type = false, $deny_bad_coords = false, $strict_nums = false, $precise = false ) {
if( is_array($coords) ) {
//passed an array, so extract the numbers from it
$matches = array('','',$coords[0]*$coords[1],'',$coords[2]*$coords[3]);
$precise = $deny_bad_coords;
} else {
//look for two floats either side of a comma (extra captures ensure array indexes remain the same as $flexi's)
$rigid = "/^(\s*)(-?[\d\.]+)\s*,(\s*)(-?[\d\.]+)\s*$/";
//look for two integers either side of a comma or space, with optional leading E/W or N/S, and trailing m
//[EW][-]<float>[m][, ][NS][-]<float>[m]
$flexi = "/^\s*([EW]?)\s*(-?[\d\.]+)(?:\s*M)?(?:\s+|\s*,\s*)([NS]?)\s*(-?[\d\.]+)(?:\s*M)?\s*$/";
if( !preg_match( $strict_nums ? $rigid : $flexi, strtoupper($coords), $matches ) ) {
//invalid format
if( $deny_bad_coords ) {
return false;
}
//assume 0,0
$matches = array( '', '', '0', '', '0' );
}
$matches[2] *= ( $matches[1] == 'W' ) ? -1 : 1;
$matches[4] *= ( $matches[3] == 'S' ) ? -1 : 1;
if( $deny_bad_coords && ( is_nan($matches[2]) || is_nan($matches[4]) ) ) {
return false;
}
}
return $this->generic_xy_coords( $return_type, $precise, $matches[2], $matches[4] );
}
//shifting grid coordinates between datums or by geoids using bilinear interpolation
private $shiftcache = array();
//shift parameter sets
//OSTN15 1 km resolution data set
private $shiftset_OSTN15 = array(
'name' => 'OSTN15',
'xspacing' => 1000,
'yspacing' => 1000,
'xstartsat' => 0,
'ystartsat' => 0,
'recordstart' => 1,
'recordcols' => 701
);
public function get_shift_set($name) {
$name = 'shiftset_'.$name;
if( isset( $this->$name ) ) {
return $this->$name;
}
return false;
}
public function create_shift_set( $name, $xspacing, $yspacing, $xstartsat = 0, $ystartsat = 0, $recordcols = 0, $recordstart = 0 ) {
if( !is_string($name) || !strlen($name) || !( $xspacing * 1 ) || !( $yspacing * 1 ) ) {
return false;
}
return array( 'name' => $name, 'xspacing' => $xspacing, 'yspacing' => $yspacing, 'xstartsat' => $xstartsat, 'ystartsat' => $ystartsat, 'recordstart' => $recordstart, 'recordcols' => $recordcols );
}
public function create_shift_record( $easting, $northing, $eastshift, $northshift, $geoidheight, $datum = null ) {
//this does not actually get used in this format - it exists only to force the right columns to exist
return array( 'name' => $easting.','.$northing, 'se' => $eastshift, 'sn' => $northshift, 'sg' => $geoidheight, 'datum' => $datum );
}
public function cache_shift_records( $name, $records ) {
if( !is_string($name) || !strlen($name) ) {
return false;
}
if( !isset($this->shiftcache[$name]) || !$this->shiftcache[$name] ) {
//nothing has ever been cached for this shift set, create the cache
$this->shiftcache[$name] = array();
}
if( !$records ) {
return true;
}
//allow it to accept either a single record, or an array of records
if( !is_array($records) || isset($records['name']) ) {
$records = array($records);
}
foreach( $records as $onerecord ) {
if( !is_array($onerecord) || !isset($onerecord['name']) || !isset($onerecord['se']) || !isset($onerecord['sn']) || !isset($onerecord['sg']) ) {
return false;
}
$this->shiftcache[$name][$onerecord['name']] = array(
'se' => $onerecord['se'],
'sn' => $onerecord['sn'],
'sg' => $onerecord['sg'],
'datum' => $onerecord['datum']
);
}
return true;
}
public function delete_shift_cache( $name = '' ) {
if( strlen($name) ) {
unset($this->shiftcache[$name]);
} else {
$this->shiftcache = array();
}
}
private function get_shift_records( $name, $records, $fetch_records ) {
if( !isset($this->shiftcache[$name]) || !$this->shiftcache[$name] ) {
//nothing has ever been cached for this shift set, create the cache
$this->shiftcache[$name] = array();
}
$unknown_records = array();
$recordnames = array();
foreach( $records as $onerecord ) {
$recordnames[] = $id = $onerecord['easting'] . ',' . $onerecord['northing'];
if( !isset($this->shiftcache[$name][$id]) ) {
$unknown_records[] = $onerecord;
$this->shiftcache[$name][$id] = false; //false gets detected as "set", null does not
}
}
//fetch records if needed
if( count($unknown_records) ) {
$this->cache_shift_records( $name, call_user_func( $fetch_records, $name, $unknown_records ) );
}
$shifts = array();
foreach( $recordnames as $id ) {
$shifts[] = $this->shiftcache[$name][$id];
}
return $shifts;
}
private function get_shifts( $x, $y, $shiftset, $fetch_records ) {
//transformation according to Transformations and OSGM15 User Guide
//https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/for-developers
if( !is_numeric($x) || !is_numeric($y) ) {
//JavaScript casts to NaN, PHP casts to 0, this forces the same response
return false;
}
//grid cell
$east_index = floor( ( $x - $shiftset['xstartsat'] ) / $shiftset['xspacing'] );
$north_index = floor( ( $y - $shiftset['ystartsat'] ) / $shiftset['yspacing'] );
//easting and northing of cell southwest corner
$x0 = $east_index * $shiftset['xspacing'] + $shiftset['xstartsat'];
$y0 = $north_index * $shiftset['yspacing'] + $shiftset['ystartsat'];
//offset within grid
$dx = $x - $x0;
$dy = $y - $y0;
//retrieve shifts for the four corners of the cell
$records = array();
$records[0] = array( 'easting' => $x0, 'northing' => $y0 );
//allow points to sit on the last line by setting them to the same values as each other (this is an improvement over the OS algorithm)
if( $dx ) {
$records[1] = array( 'easting' => $x0 + $shiftset['xspacing'], 'northing' => $y0 );
} else {
$records[1] = $records[0];
}
if( $dy ) {
$records[2] = array( 'easting' => $records[1]['easting'], 'northing' => $y0 + $shiftset['yspacing'] );
$records[3] = array( 'easting' => $x0, 'northing' => $records[2]['northing'] );
} else {
$records[2] = $records[1];
$records[3] = $records[0];
}
if( $shiftset['recordcols'] ) {
//unique, sequential record number starting from southwest eg. OSTN15
$records[0]['record'] = $east_index + ( $north_index * $shiftset['recordcols'] ) + $shiftset['recordstart'];
$records[1]['record'] = $dx ? ( $records[0]['record'] + 1 ) : $records[0]['record'];
$records[3]['record'] = $dy ? ( $records[0]['record'] + $shiftset['recordcols'] ) : $records[0]['record'];
$records[2]['record'] = $dx ? ( $records[3]['record'] + 1 ) : $records[3]['record'];
}
//fetch records if needed
$shifts = $this->get_shift_records( $shiftset['name'], $records, $fetch_records );
if( !$shifts[0] || !$shifts[1] || !$shifts[2] || !$shifts[3] ) {
return false;
}
//offset values
$t = $dx / $shiftset['xspacing'];
$u = $dy / $shiftset['yspacing'];
//bilinear interpolated shifts
$se = ( ( 1 - $t ) * ( 1 - $u ) * $shifts[0]['se'] ) + ( $t * ( 1 - $u ) * $shifts[1]['se'] ) + ( $t * $u * $shifts[2]['se'] ) + ( ( 1 - $t ) * $u * $shifts[3]['se'] );
$sn = ( ( 1 - $t ) * ( 1 - $u ) * $shifts[0]['sn'] ) + ( $t * ( 1 - $u ) * $shifts[1]['sn'] ) + ( $t * $u * $shifts[2]['sn'] ) + ( ( 1 - $t ) * $u * $shifts[3]['sn'] );
$sg = ( ( 1 - $t ) * ( 1 - $u ) * $shifts[0]['sg'] ) + ( $t * ( 1 - $u ) * $shifts[1]['sg'] ) + ( $t * $u * $shifts[2]['sg'] ) + ( ( 1 - $t ) * $u * $shifts[3]['sg'] );
$datums = null;
if( isset($shifts[0]['datum']) || isset($shifts[1]['datum']) || isset($shifts[2]['datum']) || isset($shifts[3]['datum']) ) {
$datums = array(
isset($shifts[0]['datum']) ? $shifts[0]['datum'] : null,
isset($shifts[1]['datum']) ? $shifts[1]['datum'] : null,
isset($shifts[2]['datum']) ? $shifts[2]['datum'] : null,
isset($shifts[3]['datum']) ? $shifts[3]['datum'] : null
);
}
return array( $se, $sn, $sg, $datums );
}
public function shift_grid( $easting, $northing, $height, $shiftset = false, $fetch_records = false, $return_type = false, $deny_bad_coords = false, $precise = false ) {
list( $easting, $northing, $height, $shiftset, $fetch_records, $return_type, $deny_bad_coords, $precise ) =
$this->expand_params( $easting, 2, 3, array( $northing, $height, $shiftset, $fetch_records, $return_type, $deny_bad_coords, $precise ) );
if( !$shiftset || !( $shiftset['xspacing'] * $shiftset['yspacing'] ) || !is_callable($fetch_records) ) {
return false;
}
$shifts = $this->get_shifts( $easting, $northing, $shiftset, $fetch_records );
if( !$shifts ) {
if( $deny_bad_coords ) {
return false;
} else {
return $this->generic_xy_coords( $return_type, $precise, NAN, NAN, is_numeric($height) ? NAN : null, array() );
}
}
//easting, northing, height
$easting += $shifts[0];
$northing += $shifts[1];
if( is_numeric($height) ) {
$height -= $shifts[2];
}
return $this->generic_xy_coords( $return_type, $precise, $easting, $northing, $height, $shifts[3] );
}
public function reverse_shift_grid( $easting, $northing, $height, $shiftset = false, $fetch_records = false, $return_type = false, $deny_bad_coords = false, $precise = false ) {
//transformation according to Transformations and OSGM15 User Guide
//https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/for-developers
list( $easting, $northing, $height, $shiftset, $fetch_records, $return_type, $deny_bad_coords, $precise ) =
$this->expand_params( $easting, 2, 3, array( $northing, $height, $shiftset, $fetch_records, $return_type, $deny_bad_coords, $precise ) );
if( !$shiftset || !( $shiftset['xspacing'] * $shiftset['yspacing'] ) || !is_callable($fetch_records) ) {
return false;
}
$shifts = array( 0, 0 );
$loopcount = 0;
do {
$oldshifts = $shifts;
//get the shifts for the current estimated point, starting with the initial easting and northing
//apply the new shift and try again with the shift from the new point
$shifts = $this->get_shifts( is_numeric($easting) ? ( $easting - $shifts[0] ) : $easting, is_numeric($northing) ? ( $northing - $shifts[1] ) : $northing, $shiftset, $fetch_records );
if( !$shifts ) {
if( $deny_bad_coords ) {
return false;
} else {
return $this->generic_xy_coords( $return_type, $precise, NAN, NAN, is_numeric($height) ? NAN : null, array() );
}
}
//see if the shift changes with each iteration
//unlikely to ever hit 100 iterations, even when taken to extremes
} while( ( abs( $oldshifts[0] - $shifts[0] ) > 0.0001 || abs( $oldshifts[1] - $shifts[1] ) > 0.0001 ) && ++$loopcount < 100 );
if( is_numeric($height) ) {
$height += $shifts[2];
}
return $this->generic_xy_coords( $return_type, $precise, $easting - $shifts[0], $northing - $shifts[1], $height, $shifts[3] );
}
//ellipsoid parameters used during grid->lat/long conversions and Helmert transformations
private $ellipsoid_Airy_1830 = array(
//Airy 1830 (OS)
'a' => 6377563.396,
'b' => 6356256.909
);
private $ellipsoid_Airy_1830_mod = array(
//Airy 1830 modified (OSI)
'a' => 6377340.189,
'b' => 6356034.447
);
private $ellipsoid_GRS80 = array(
//GRS80 (ETRS89, early versions of GPS)
'a' => 6378137,
'b' => 6356752.3141
);
private $ellipsoid_WGS84 = array(
//WGS84 (GPS)
'a' => 6378137,
'b' => 6356752.31425 // https://georepository.com/ellipsoid_7030/WGS-84.html 4 November 2020
);
private $datumset_OSGB36 = array(
//Airy 1830 (OS)
'a' => 6377563.396,
'b' => 6356256.909,
'F0' => 0.9996012717,
'E0' => 400000, 'N0' => -100000,
'Phi0' => 49,
'Lambda0' => -2
);
private $datumset_OSTN15 = array(
//OSGB36 parameters with GRS80 ellipsoid (OSTN15)
'a' => 6378137,
'b' => 6356752.3141,
'F0' => 0.9996012717,
'E0' => 400000, 'N0' => -100000,
'Phi0' => 49,
'Lambda0' => -2
);
private $datumset_Ireland_1965 = array(
//Airy 1830 modified (OSI)
'a' => 6377340.189,
'b' => 6356034.447,
'F0' => 1.000035,
'E0' => 200000, 'N0' => 250000,
'Phi0' => 53.5,
'Lambda0' => -8
);
private $datumset_IRENET95 = array(
//ITM (OSI) taken from http://en.wikipedia.org/wiki/Irish_Transverse_Mercator
//uses GRS80, which is so close to WGS84 that there is no point in transforming with Helmert before projecting
'a' => 6378137,
'b' => 6356752.3141,
'F0' => 0.999820,
'E0' => 600000, 'N0' => 750000,
'Phi0' => 53.5,
'Lambda0' => -8
);
//UPS (uses WGS84), taken from http://www.epsg.org/guides/ number 7 part 2 "Coordinate Conversions and Transformations including Formulas"
//officially defined in http://earth-info.nga.mil/GandG/publications/tm8358.2/TM8358_2.pdf
private $datumset_UPS_North = array(
'a' => 6378137,
'b' => 6356752.31425,
'F0' => 0.994,
'E0' => 2000000, 'N0' => 2000000,
'Phi0' => 90,
'Lambda0' => 0
);
private $datumset_UPS_South = array(
'a' => 6378137,
'b' => 6356752.31425,
'F0' => 0.994,
'E0' => 2000000, 'N0' => 2000000,
'Phi0' => -90,
'Lambda0' => 0
);
public function get_ellipsoid($name) {
$name = 'ellipsoid_'.$name;
if( isset( $this->$name ) ) {
return $this->$name;
}
return false;
}
public function create_ellipsoid( $a, $b ) {
return array('a'=>$a,'b'=>$b);
}
public function get_datum($name) {
$name = 'datumset_'.$name;
if( isset( $this->$name ) ) {
return $this->$name;
}
return false;
}
public function create_datum( $ellip, $F0, $E0, $N0, $Phi0, $Lambda0 ) {
if( !isset($ellip['a']) ) {
return false;
}
return array('a'=>$ellip['a'],'b'=>$ellip['b'],'F0'=>$F0,'E0'=>$E0,'N0'=>$N0,'Phi0'=>$Phi0,'Lambda0'=>$Lambda0);
}
//conversions between 12345,67890 grid references and latitude/longitude formats using transverse mercator projection
public $COORDS_OS_UK = 1;
public $COORDS_OSI = 2;
public $COORDS_GPS_UK = 3;
public $COORDS_GPS_IRISH = 4;
public $COORDS_GPS_ITM = 5;
public $COORDS_GPS_IRISH_HELMERT = 6;
public function grid_to_lat_long( $E, $N, $type = false, $return_type = false ) {
//horribly complex conversion according to "A guide to coordinate systems in Great Britain" Annexe C:
//http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/
//http://www.movable-type.co.uk/scripts/latlong-gridref.html shows an alternative script for JS, which also says what some OS variables represent
list( $E, $N, $type, $return_type ) = $this->expand_params( $E, 2, 2, array( $N, $type, $return_type ) );
//get appropriate ellipsoid semi-major axis 'a' (metres) and semi-minor axis 'b' (metres),
//grid scale factor on central meridian, and true origin (grid and lat-long) from Annexe A
//extracts variables called $a,$b,$F0,$E0,$N0,$Phi0,$Lambda0
if( is_array($type) ) {
extract( $type, EXTR_SKIP );
} elseif( $type == $this->COORDS_OS_UK || $type == $this->COORDS_GPS_UK ) {
extract( $this->datumset_OSGB36, EXTR_SKIP );
} elseif( $type == $this->COORDS_GPS_ITM ) {
extract( $this->datumset_IRENET95, EXTR_SKIP );
} elseif( $type == $this->COORDS_OSI || $type == $this->COORDS_GPS_IRISH || $type == $this->COORDS_GPS_IRISH_HELMERT ) {
extract( $this->datumset_Ireland_1965, EXTR_SKIP );
}
if( !isset($F0) ) {
//invalid type
return false;
}
//PHP will not allow expressions in the arrays as they are defined inline as class properties, so do the conversion to radians here
$Phi0 = deg2rad($Phi0);
//eccentricity-squared from Annexe B B1
//$e2 = ( ( $a * $a ) - ( $b * $b ) ) / ( $a * $a );
$e2 = 1 - ( $b * $b ) / ( $a * $a ); //optimised
//C1
$n = ( $a - $b ) / ( $a + $b );
//pre-compute values that will be re-used many times in the C3 formula
$n2 = $n * $n;
$n3 = pow( $n, 3 );
$n_parts1 = ( 1 + $n + 1.25 * $n2 + 1.25 * $n3 );
$n_parts2 = ( 3 * $n + 3 * $n2 + 2.625 * $n3 );
$n_parts3 = ( 1.875 * $n2 + 1.875 * $n3 );
$n_parts4 = ( 35 / 24 ) * $n3;
//iterate to find latitude (when $N - $N0 - $M < 0.01mm)
$Phi = $Phi0;
$M = 0;
$loopcount = 0;
do {
//C6 and C7
$Phi += ( ( $N - $N0 - $M ) / ( $a * $F0 ) );
//C3
$M = $b * $F0 * (
$n_parts1 * ( $Phi - $Phi0 ) -
$n_parts2 * sin( $Phi - $Phi0 ) * cos( $Phi + $Phi0 ) +
$n_parts3 * sin( 2 * ( $Phi - $Phi0 ) ) * cos( 2 * ( $Phi + $Phi0 ) ) -
$n_parts4 * sin( 3 * ( $Phi - $Phi0 ) ) * cos( 3 * ( $Phi + $Phi0 ) )
); //meridonal arc
//due to number precision, it is possible to get infinite loops here for extreme cases (especially for invalid ellipsoid numbers)
//in tests, upto 6 loops are needed for grid 25 times Earth circumference - if it reaches 100, assume it must be infinite, and break out
} while( abs( $N - $N0 - $M ) >= 0.00001 && ++$loopcount < 100 ); //0.00001 == 0.01 mm
//pre-compute values that will be re-used many times in the C2 and C8 formulae
$sin_Phi = sin( $Phi );
$sin2_Phi = $sin_Phi * $sin_Phi;
$tan_Phi = tan( $Phi );
$sec_Phi = 1 / cos( $Phi );
$tan2_Phi = $tan_Phi * $tan_Phi;
$tan4_Phi = $tan2_Phi * $tan2_Phi;
//C2
$Rho = $a * $F0 * ( 1 - $e2 ) * pow( 1 - $e2 * $sin2_Phi, -1.5 ); //meridional radius of curvature
$Nu = $a * $F0 / sqrt( 1 - $e2 * $sin2_Phi ); //transverse radius of curvature
$Eta2 = $Nu / $Rho - 1;
//pre-compute more values that will be re-used many times in the C8 formulae
$Nu3 = pow( $Nu, 3 );
$Nu5 = pow( $Nu, 5 );
//C8 parts
$VII = $tan_Phi / ( 2 * $Rho * $Nu );
$VIII = ( $tan_Phi / ( 24 * $Rho * $Nu3 ) ) * ( 5 + 3 * $tan2_Phi + $Eta2 - 9 * $tan2_Phi * $Eta2 );
$IX = ( $tan_Phi / ( 720 * $Rho * $Nu5 ) ) * ( 61 + 90 * $tan2_Phi + 45 * $tan4_Phi );
$X = $sec_Phi / $Nu;
$XI = ( $sec_Phi / ( 6 * $Nu3 ) ) * ( ( $Nu / $Rho ) + 2 * $tan2_Phi );
$XII = ( $sec_Phi / ( 120 * $Nu5 ) ) * ( 5 + 28 * $tan2_Phi + 24 * $tan4_Phi );
$XIIA = ( $sec_Phi / ( 5040 * pow( $Nu, 7 ) ) ) * ( 61 + 662 * $tan2_Phi + 1320 * $tan4_Phi + 720 * pow( $tan_Phi, 6 ) );
//C8, C9
$Edif = $E - $E0;
$latitude = rad2deg( $Phi - $VII * $Edif * $Edif + $VIII * pow( $Edif, 4 ) - $IX * pow( $Edif, 6 ) );
$longitude = $Lambda0 + rad2deg( $X * $Edif - $XI * pow( $Edif, 3 ) + $XII * pow( $Edif, 5 ) - $XIIA * pow( $Edif, 7 ) );
if( $type == $this->COORDS_GPS_UK ) {
list( $latitude, $longitude ) = $this->Helmert_transform( $latitude, $longitude, $this->ellipsoid_Airy_1830, $this->Helmert_OSGB36_to_WGS84, $this->ellipsoid_WGS84 );
} elseif( $type == $this->COORDS_GPS_IRISH ) {
list( $latitude, $longitude ) = $this->polynomial_transform( $latitude, $longitude, $this->polycoeff_OSiLPS );
} elseif( $type == $this->COORDS_GPS_IRISH_HELMERT ) {
list( $latitude, $longitude ) = $this->Helmert_transform( $latitude, $longitude, $this->ellipsoid_Airy_1830_mod, $this->Helmert_Ireland65_to_WGS84, $this->ellipsoid_WGS84 );
}
return $this->generic_lat_long_coords( $return_type, $latitude, $longitude );
}
public function lat_long_to_grid( $latitude, $longitude, $type = false, $return_type = false, $precise = false) {
//horribly complex conversion according to "A guide to coordinate systems in Great Britain" Annexe C:
//http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/
//http://www.movable-type.co.uk/scripts/latlong-gridref.html shows an alternative script for JS, which also says what some OS variables represent
list( $latitude, $longitude, $type, $return_type, $precise ) = $this->expand_params( $latitude, 2, 2, array( $longitude, $type, $return_type, $precise ) );
//convert back to local ellipsoid coordinates
if( $type == $this->COORDS_GPS_UK ) {
list( $latitude, $longitude ) = $this->Helmert_transform( $latitude, $longitude, $this->ellipsoid_WGS84, $this->Helmert_WGS84_to_OSGB36, $this->ellipsoid_Airy_1830 );
} elseif( $type == $this->COORDS_GPS_IRISH ) {
list( $latitude, $longitude ) = $this->reverse_polynomial_transform( $latitude, $longitude, $this->polycoeff_OSiLPS );
} elseif( $type == $this->COORDS_GPS_IRISH_HELMERT ) {
list( $latitude, $longitude ) = $this->Helmert_transform( $latitude, $longitude, $this->ellipsoid_WGS84, $this->Helmert_WGS84_to_Ireland65, $this->ellipsoid_Airy_1830_mod );
}
//get appropriate ellipsoid semi-major axis 'a' (metres) and semi-minor axis 'b' (metres),
//grid scale factor on central meridian, and true origin (grid and lat-long) from Annexe A
//extracts variables called $a,$b,$F0,$E0,$N0,$Phi0,$Lambda0
if( is_array($type) ) {
extract( $type, EXTR_SKIP );
} elseif( $type == $this->COORDS_OS_UK || $type == $this->COORDS_GPS_UK ) {
extract( $this->datumset_OSGB36, EXTR_SKIP );
} elseif( $type == $this->COORDS_GPS_ITM ) {
extract( $this->datumset_IRENET95, EXTR_SKIP );
} elseif( $type == $this->COORDS_OSI || $type == $this->COORDS_GPS_IRISH || $type == $this->COORDS_GPS_IRISH_HELMERT ) {
extract( $this->datumset_Ireland_1965, EXTR_SKIP );
}
if( !isset($F0) ) {
//invalid type
return false;
}
//PHP will not allow expressions in the arrays as they are defined inline as class properties, so do the conversion to radians here
$Phi0 = deg2rad($Phi0);
$Phi = deg2rad($latitude);
$Lambda = $longitude - $Lambda0;
//force Lambda between -180 and 180
$Lambda = $this->wrap_around_longitude($Lambda);
$Lambda = deg2rad($Lambda);
//eccentricity-squared from Annexe B B1
//$e2 = ( ( $a * $a ) - ( $b * $b ) ) / ( $a * $a );
$e2 = 1 - ( $b * $b ) / ( $a * $a ); //optimised
//C1
$n = ( $a - $b ) / ( $a + $b );
//pre-compute values that will be re-used many times in the C2, C3 and C4 formulae
$sin_Phi = sin( $Phi );
$sin2_Phi = $sin_Phi * $sin_Phi;
$cos_Phi = cos( $Phi );
$cos3_Phi = pow( $cos_Phi, 3 );
$cos5_Phi = pow( $cos_Phi, 5 );
$tan_Phi = tan( $Phi );
$tan2_Phi = $tan_Phi * $tan_Phi;
$tan4_Phi = $tan2_Phi * $tan2_Phi;
$n2 = $n * $n;
$n3 = pow( $n, 3 );
//C2
$Nu = $a * $F0 / sqrt( 1 - $e2 * $sin2_Phi ); //transverse radius of curvature
$Rho = $a * $F0 * ( 1 - $e2 ) * pow( 1 - $e2 * $sin2_Phi, -1.5 ); //meridional radius of curvature
$Eta2 = $Nu / $Rho - 1;
//C3, meridonal arc
$M = $b * $F0 * (
( 1 + $n + 1.25 * $n2 + 1.25 * $n3 ) * ( $Phi - $Phi0 ) -
( 3 * $n + 3 * $n2 + 2.625 * $n3 ) * sin( $Phi - $Phi0 ) * cos( $Phi + $Phi0 ) +
( 1.875 * $n2 + 1.875 * $n3 ) * sin( 2 * ( $Phi - $Phi0 ) ) * cos( 2 * ( $Phi + $Phi0 ) ) -
( 35 / 24 ) * $n3 * sin( 3 * ( $Phi - $Phi0 ) ) * cos( 3 * ( $Phi + $Phi0 ) )
);
//C4
$I = $M + $N0;
$II = ( $Nu / 2 ) * $sin_Phi * $cos_Phi;
$III = ( $Nu / 24 ) * $sin_Phi * $cos3_Phi * ( 5 - $tan2_Phi + 9 * $Eta2 );
$IIIA = ( $Nu / 720 ) * $sin_Phi * $cos5_Phi * ( 61 - 58 * $tan2_Phi + $tan4_Phi );
$IV = $Nu * $cos_Phi;
$V = ( $Nu / 6 ) * $cos3_Phi * ( ( $Nu / $Rho ) - $tan2_Phi );
$VI = ( $Nu / 120 ) * $cos5_Phi * ( 5 - 18 * $tan2_Phi + $tan4_Phi + 14 * $Eta2 - 58 * $tan2_Phi * $Eta2 );
$N = $I + $II * $Lambda * $Lambda + $III * pow( $Lambda, 4 ) + $IIIA * pow( $Lambda, 6 );
$E = $E0 + $IV * $Lambda + $V * pow( $Lambda, 3 ) + $VI * pow( $Lambda, 5 );
return $this->generic_xy_coords( $return_type, $precise, $E, $N );
}
//UTM - freaky format consisting of 60 transverse mercators
public function utm_to_lat_long( $zone, $north = false, $x = false, $y = false, $ellip = false, $return_type = false, $deny_bad_reference = false ) {
list( $zone, $north, $x, $y, $ellip, $return_type, $deny_bad_reference ) =
$this->expand_params( $zone, 4, 4, array( $north, $x, $y, $ellip, $return_type, $deny_bad_reference ) );
if( is_string( $zone ) ) {
//manually move params
if( !is_array($ellip) ) {
//if it is an array, expand_params will have already moved it over
$deny_bad_reference = $y;
$return_type = $x;
$ellip = $north;
}
$zone = strtoupper($zone);
if( preg_match( "/^\s*(0?[1-9]|[1-5][0-9]|60)\s*([C-HJ-NP-X]|NORTH|SOUTH)\s*(-?[\d\.]+)\s*[,\s]\s*(-?[\d\.]+)\s*$/", $zone, $parsed_zone ) ) {
//matched the shorthand 30U 1234 5678
//[01-60](<letter>|north|south)<float>[, ]<float>
//radix parameter is needed since numbers often start with a leading 0 and must not be treated as octal
$zone = $parsed_zone[1] * 1;
if( strlen( $parsed_zone[2] ) > 1 ) {
$north = ( $parsed_zone[2] == 'NORTH' ) ? 1 : -1;
} else {
$north = ( $parsed_zone[2] > 'M' ) ? 1 : -1;
}
$x = $parsed_zone[3] * 1;
$y = $parsed_zone[4] * 1;
} else if( preg_match( "/^\s*(-?[\d\.]+)\s*[A-Z]*\s*[,\s]\s*(-?[\d\.]+)\s*[A-Z]*\s*[\s,]\s*ZONE\s*(0?[1-9]|[1-5][0-9]|60)\s*,?\s*([NS])/", $zone, $parsed_zone ) ) {
//matched the longhand 630084 mE 4833438 mN, zone 17, Northern Hemisphere
//<float>[letters][, ]<float>[letters][, ]zone[01-60][,][NS]...
$zone = $parsed_zone[3] * 1;
$north = ( $parsed_zone[4] == 'N' ) ? 1 : -1;
$x = $parsed_zone[1] * 1;
$y = $parsed_zone[2] * 1;
} else {
//make it reject it
$zone = 0;
}
}
if( is_array($ellip) && !isset($ellip['a']) ) {
//invalid ellipsoid must be rejected early, before returning error values or constructing a datum
return false;
}
if( !is_numeric($zone) || $zone < 1 || $zone > 60 || !is_numeric($x) || !is_numeric($y) ) {
if( $deny_bad_reference ) {
//invalid format
return false;
}
//default coordinates put it at 90,0 lat/long - use dms_to_dd to get the right return_type
return $this->dms_to_dd(array(90,0,0,1,0,0,0,0),$return_type);
}
$ellipsoid = is_array($ellip) ? $ellip : $this->ellipsoid_WGS84;
$ellipsoid = array( 'a' => $ellipsoid['a'], 'b' => $ellipsoid['b'], 'F0' => 0.9996, 'E0' => 500000, 'N0' => ( $north < 0 ) ? 10000000 : 0, 'Phi0' => 0, 'Lambda0' => ( 6 * $zone ) - 183 );
return $this->grid_to_lat_long($x,$y,$ellipsoid,$return_type);
}
public function lat_long_to_utm( $latitude, $longitude = false, $ellip = false, $format = false, $return_type = false, $deny_bad_coords = false, $precise = false ) {
list( $latitude, $longitude, $ellip, $format, $return_type, $deny_bad_coords, $precise ) =
$this->expand_params( $latitude, 2, 2, array( $longitude, $ellip, $format, $return_type, $deny_bad_coords, $precise ) );
if( is_array($ellip) && !isset($ellip['a']) ) {
//invalid ellipsoid must be rejected early, before returning error values or constructing a datum
return false;
}
//force the longitude between -180 and 179.99...9
$longitude = $this->wrap_around_longitude( $longitude, true );
if( !is_numeric($longitude) || !is_numeric($latitude) || $latitude > 84 || $latitude < -80 ) {
if( $deny_bad_coords ) {
//invalid format
return false;
}
//default coordinates put it at ~0,0 lat/long
if( !is_numeric($longitude) ) {
$longitude = 0;
}
if( !is_numeric($latitude) ) {
$latitude = 0;
}
if( $latitude > 84 ) {
//out of range, return appropriate polar letter, and bail out
$zoneletter = $format ? 'North' : ( ( $longitude < 0 ) ? 'Y' : 'Z' );
$zone = $x = $y = 0;
}
if( $latitude < -80 ) {
//out of range, return appropriate polar letter, and bail out
$zoneletter = $format ? 'South' : ( ( $longitude < 0 ) ? 'A' : 'B' );
$zone = $x = $y = 0;
}
}
if( !isset($zoneletter) || !$zoneletter ) {
//add hacks to work out if it lies in the non-standard zones
if( $latitude >= 72 && $longitude >= 6 && $longitude < 36 ) {
//band X, these parts are moved around
if( $longitude < 9 ) {
$zone = 31;
} else if( $longitude < 21 ) {
$zone = 33;
} else if( $longitude < 33 ) {
$zone = 35;
} else {
$zone = 37;
}
} else if( $latitude >= 56 && $latitude < 64 && $longitude >= 3 && $longitude < 6 ) {
//band U, this part of zone 31 is moved into zone 32
$zone = 32;
} else {
//yay for standards!
$zone = floor( $longitude / 6 ) + 31;
}
//get the band letter
if( $format ) {
$zoneletter = ( $latitude < 0 ) ? 'South' : 'North';
} else {
$zoneletter = floor( $latitude / 8 ) + 77; //67 is ASCII C
if( $zoneletter > 72 ) {
//skip I
$zoneletter++;
}
if( $zoneletter > 78 ) {
//skip O
$zoneletter++;
}
if( $zoneletter > 88 ) {
//X is as high as it goes
$zoneletter = 88;
}
$zoneletter = chr($zoneletter);
}
//do actual transformation - happens inside "if(!$zoneletter)" because this is the non-error path
$ellipsoid = is_array($ellip) ? $ellip : $this->ellipsoid_WGS84;
$ellipsoid = array( 'a' => $ellipsoid['a'], 'b' => $ellipsoid['b'], 'F0' => 0.9996, 'E0' => 500000, 'N0' => ( $latitude < 0 ) ? 10000000 : 0, 'Phi0' => 0, 'Lambda0' => ( 6 * $zone ) - 183 );
$tmpcoords = $this->lat_long_to_grid($latitude,$longitude,$ellipsoid);
if( !$tmpcoords ) { return false; }
$x = $tmpcoords[0];
$y = $tmpcoords[1];
}
if( $return_type ) {
if( $precise ) {
$x = sprintf( '%F', $x );
$y = sprintf( '%F', $y );
} else {
//avoid stupid -0 by adding 0
$x = round($x) + 0;
$y = round($y) + 0;
}
if( $format ) {
return $x . 'mE, ' . $y . 'mN, Zone ' . $zone . ', ' . $zoneletter . 'ern Hemisphere';
}
return $zone . $zoneletter . ' ' . $x . ' ' . $y;
} else {
return array( $zone, ( $latitude < 0 ) ? -1 : 1, $x+0, $y+0, $zoneletter );
}
}
//basic polar stereographic pojections
//formulae according to http://www.epsg.org/guides/ number 7 part 2 "Coordinate Conversions and Transformations including Formulas"
public function polar_to_lat_long( $easting, $northing, $datum = false, $return_type = false ) {
list( $easting, $northing, $datum, $return_type ) = $this->expand_params( $easting, 2, 2, array( $northing, $datum, $return_type ) );
if( !$datum ) {
return false;
}
if( !isset($datum['F0']) || ( $datum['Phi0'] != 90 && $datum['Phi0'] != -90 ) ) {
//invalid type
return false;
}
$a = $datum['a'];
$b = $datum['b'];
$k0 = $datum['F0'];
$FE = $datum['E0'];
$FN = $datum['N0'];
$Phi0 = $datum['Phi0'];
$Lambda0 = $datum['Lambda0'];
//eccentricity-squared
$e2 = 1 - ( $b * $b ) / ( $a * $a ); //optimised
$e = sqrt($e2);
$Rho = sqrt( pow( $easting - $FE, 2 ) + pow( $northing - $FN, 2 ) );
$t = $Rho * sqrt( pow( 1 + $e, 1 + $e ) * pow( 1 - $e, 1 - $e ) ) / ( 2 * $a * $k0 );
if( $Phi0 < 0 ) {
//south
$x = 2 * atan($t) - M_PI / 2;
} else {
//north
$x = M_PI / 2 - 2 * atan($t);
}
//pre-compute values that will be re-used many times in the Phi formula
$e4 = $e2 * $e2;
$e6 = $e4 * $e2;
$e8 = $e4 * $e4;
$Phi = $x + ( $e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360 ) * sin( 2 * $x ) +
( 7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520 ) * sin( 4 * $x ) +
( 7 * $e6 / 120 + 81 * $e8 / 1120 ) * sin( 6 * $x ) +
( 4279 * $e8 / 161280 ) * sin( 8 * $x );
//longitude
//formulas here are wrong in the epsg guide; atan(foo/bar) should have been atan2(foo,bar) or it is wrong for half of the grid
if( $Phi0 < 0 ) {
//south
$Lambda = $Lambda0 + atan2( $easting - $FE, $northing - $FN );
} else {
//north
$Lambda = $Lambda0 + atan2( $easting - $FE, $FN - $northing );
}
$latitude = rad2deg($Phi);
$longitude = rad2deg($Lambda) + $Lambda0;
return $this->generic_lat_long_coords( $return_type, $latitude, $longitude );
}
public function lat_long_to_polar( $latitude, $longitude, $datum = false, $return_type = false, $precise = false ) {
list( $latitude, $longitude, $datum, $return_type, $precise ) = $this->expand_params( $latitude, 2, 2, array( $longitude, $datum, $return_type, $precise ) );
if( !$datum ) {
return false;
}
if( !isset($datum['F0']) || ( $datum['Phi0'] != 90 && $datum['Phi0'] != -90 ) ) {
//invalid type
return false;
}
$a = $datum['a'];
$b = $datum['b'];
$k0 = $datum['F0'];
$FE = $datum['E0'];
$FN = $datum['N0'];
$Lambda0 = $datum['Lambda0'];
$Phi = deg2rad($latitude);
$Lambda = deg2rad( $longitude - $Lambda0 );
//eccentricity-squared
$e2 = 1 - ( $b * $b ) / ( $a * $a ); //optimised
$e = sqrt($e2);
//t
$sinPhi = sin( $Phi );
if( $latitude < 0 ) {
//south
$t = tan( ( M_PI / 4 ) + ( $Phi / 2 ) ) / pow( ( 1 + $e * $sinPhi ) / ( 1 - $e * $sinPhi ), $e / 2 );
} else {
//north
$t = tan( ( M_PI / 4 ) - ( $Phi / 2 ) ) * pow( ( 1 + $e * $sinPhi ) / ( 1 - $e * $sinPhi ), $e / 2 );
}
//Rho
$Rho = 2 * $a * $k0 * $t / sqrt( pow( 1 + $e, 1 + $e ) * pow( 1 - $e, 1 - $e ) );
if( $latitude < 0 ) {
//south
$N = $FN + $Rho * cos( $Lambda );
} else {
//north - origin is *down*
$N = $FN - $Rho * cos( $Lambda );
}
$E = $FE + $Rho * sin( $Lambda );
return $this->generic_xy_coords( $return_type, $precise, $E, $N );
}
public function ups_to_lat_long( $hemisphere, $x = false, $y = false, $return_type = false, $deny_bad_reference = false, $min_length = false ) {
list( $hemisphere, $x, $y, $return_type, $deny_bad_reference, $min_length ) =
$this->expand_params( $hemisphere, 3, 3, array( $x, $y, $return_type, $deny_bad_reference, $min_length ) );
if( !is_numeric($x) || !is_numeric($y) ) {
//a single string 'X 12345 67890', split into parts
//manually move params
$min_length = $return_type;
$deny_bad_reference = $y;
$return_type = $x;
//(A|B|Y|Z|N|S|north|south)[,]<float>[, ]<float>
if( !preg_match( "/^\s*([ABNSYZ]|NORTH|SOUTH)\s*,?\s*(-?[\d\.]+)\s*[\s,]\s*(-?[\d\.]+)\s*$/i", $hemisphere, $hemisphere ) ) {
//invalid format, will get handled later
$x = $y = null;
} else {
$x = $hemisphere[2];
$y = $hemisphere[3];
$hemisphere = $hemisphere[1];
}
}
if( !is_numeric($x) || !is_numeric($y) || ( !is_numeric($hemisphere) && ( !is_string($hemisphere) || !preg_match( "/^([ABNSYZ]|NORTH|SOUTH)$/i", $hemisphere ) ) ) || ( $min_length && ( $x < 800000 || $y < 800000 ) ) ) {
if( $deny_bad_reference ) {
//invalid format
return false;
}
//default coordinates put it at 0,0 lat/long - use dms_to_dd to get the right return_type
return $this->dms_to_dd(array(0,0,0,0,0,0,0,0),$return_type);
}
if( is_string($hemisphere) ) {
$hemisphere = strtoupper($hemisphere);
} else {
$hemisphere = ( $hemisphere < 0 ) ? 'S' : 'N';
}
if( $hemisphere == 'N' || $hemisphere == 'NORTH' || $hemisphere == 'Y' || $hemisphere == 'Z' ) {
$hemisphere = $this->datumset_UPS_North;
} else {
$hemisphere = $this->datumset_UPS_South;
}
return $this->polar_to_lat_long($x*1,$y*1,$hemisphere,$return_type);
}
public function lat_long_to_ups( $latitude, $longitude = false, $format = false, $return_type = false, $deny_bad_coords = false, $precise = false ) {
list( $latitude, $longitude, $format, $return_type, $deny_bad_coords, $precise ) =
$this->expand_params( $latitude, 2, 2, array( $longitude, $format, $return_type, $deny_bad_coords, $precise ) );
//force the longitude between -179.99...9 and 180
$longitude = $this->wrap_around_longitude($longitude);
if( !is_numeric($longitude) || !is_numeric($latitude) || $latitude > 90 || $latitude < -90 || ( $latitude < 83.5 && $latitude > -79.5 ) ) {
if( $deny_bad_coords ) {
//invalid format
return false;
}
//default coordinates put it as 90,0 in OUT_OF_GRID zone
$tmp = array( 2000000, 2000000 );
$letter = 'OUT_OF_GRID';
} else {
$tmp = $this->lat_long_to_polar( $latitude, $longitude, ( $latitude < 0 ) ? $this->datumset_UPS_South : $this->datumset_UPS_North );
if( $latitude < 0 ) {
if( $format ) {
$letter = 'S';
} else if( $longitude < 0 ) {
$letter = 'A';
} else {
$letter = 'B';
}
} else {
if( $format ) {
$letter = 'N';
} else if( $longitude < 0 ) {
$letter = 'Y';
} else {
$letter = 'Z';
}
}
}
if( $return_type ) {
if( $precise ) {
$tmp[0] = sprintf( '%F', $tmp[0] );
$tmp[1] = sprintf( '%F', $tmp[1] );
} else {
//avoid stupid -0 by adding 0
$tmp[0] = round($tmp[0]) + 0;
$tmp[1] = round($tmp[1]) + 0;
}
return $letter.' '.$tmp[0].' '.$tmp[1];
} else {
return array( ( $latitude < 0 ) ? -1 : 1, $tmp[0]+0, $tmp[1]+0, $letter );
}
}
//conversions between different latitude/longitude formats
public function dd_to_dms( $N, $E = false, $only_dm = false, $return_type = false ) {
//decimal degrees (49.5,-123.5) to degrees-minutes-seconds (49°30'0"N, 123°30'0"W)
list( $N, $E, $only_dm, $return_type ) = $this->expand_params( $N, 2, 2, array( $E, $only_dm, $return_type ) );
$N_abs = abs($N);
$E_abs = abs($E);
$degN = floor($N_abs);
$degE = floor($E_abs);
if( $only_dm ) {
$minN = 60 * ( $N_abs - $degN );
$secN = 0;
$minE = 60 * ( $E_abs - $degE );
$secE = 0;
} else {
//the approach used here is careful to respond consistently to floating point errors for all of degrees/minutes/seconds
//errors should not cause one to be increased while another is decreased (which could cause eg. 5 minutes 60 seconds)
$minN = 60 * $N_abs;
$secN = ( $minN - floor( $minN ) ) * 60;
$minN %= 60;
$minE = 60 * $E_abs;
$secE = ( $minE - floor( $minE ) ) * 60;
$minE %= 60;
}
//avoid stupid -0 by adding 0
$degN += 0;
$minN += 0;
$secN += 0;
$degE += 0;
$minE += 0;
$secE += 0;
if( $return_type ) {
//sprintf to produce simple numbers instead of scientific notation (also reduces accuracy to 6 decimal places)
$deg = is_string($return_type) ? $return_type : '°';
$quot = is_string($return_type) ? '"' : '"';
if( $only_dm ) {
//careful not to round up to 60 minutes when displaying
return
$degN . $deg . sprintf( '%.8F', ( $minN >= 59.999999995 ) ? 59.99999999 : $minN ) . "'" . ( ( $N < 0 ) ? 'S' : 'N' ) . ', ' .
$degE . $deg . sprintf( '%.8F', ( $minE >= 59.999999995 ) ? 59.99999999 : $minE ) . "'" . ( ( $E < 0 ) ? 'W' : 'E' );
} else {
//careful not to round up to 60 seconds when displaying
return
$degN . $deg . $minN . "'" . sprintf( '%F', ( $secN >= 59.9999995 ) ? 59.999999 : $secN ) . $quot . ( ( $N < 0 ) ? 'S' : 'N' ) . ', ' .
$degE . $deg . $minE . "'" . sprintf( '%F', ( $secE >= 59.9999995 ) ? 59.999999 : $secE ) . $quot . ( ( $E < 0 ) ? 'W' : 'E' );
}
} else {
return array( $degN, $minN, $secN, ( $N < 0 ) ? -1 : 1, $degE, $minE, $secE, ( $E < 0 ) ? -1 : 1 );
}
}
public function dd_format( $N, $E = false, $no_units = false, $return_type = false ) {
//different formats of decimal degrees (49.5,-123.5)
list( $N, $E, $no_units, $return_type ) = $this->expand_params( $N, 2, 2, array( $E, $no_units, $return_type ) );
if( $no_units && $return_type ) {
return $this->generic_lat_long_coords( $return_type, $N, $E );
}
$E = $this->wrap_around_longitude($E);
//rounding needs to happen before testing for signs and letters, rounds *down* negative numbers just like sprintf
if( $return_type && $E <= -179.99999999995 ) {
$E = 180;
}
if( $no_units ) {
$lat_mul = $long_mul = 1;
} else {
$lat_mul = $return_type ? ( ( $N < 0 ) ? 'S' : 'N' ) : ( ( $N < 0 ) ? -1 : 1 );
$long_mul = $return_type ? ( ( $E < 0 ) ? 'W' : 'E' ) : ( ( $E < 0 ) ? -1 : 1 );
$N = abs($N);
$E = abs($E);
}
if( $return_type ) {
$deg = is_string($return_type) ? $return_type : '°';
return sprintf( '%.10F', $N ) . $deg . $lat_mul . ', ' . sprintf( '%.10F', $E ) . $deg . $long_mul;
}
return array( $N, 0, 0, $lat_mul, $E, 0, 0, $long_mul );
}
public function dms_to_dd( $dms, $return_type = false, $deny_bad_coords = false, $allow_unitless = false ) {
//degrees-minutes-seconds (49°30'0"N, 123°30'0"W) to decimal degrees (49.5,-123.5)
if( is_array( $dms ) ) {
//passed an array of values, which can be unshifted once to get the right positions
$latlong = $dms;
array_unshift($latlong,'','');
array_splice($latlong,6,0,'');
} else {
$latlong = null;
//matches [-]<float>,[-]<float> which could also be grid coordinates, so this is kept behind an option
if( $allow_unitless && preg_match("/^\s*(-?)([\d\.]+)\s*,\s*(-?)([\d\.]+)\s*$/", strtoupper($dms), $latlong ) ) {
$latlong = array('',$latlong[1],$latlong[2],'0','0','',$latlong[3],$latlong[4],'0','0','');
}
//matches [-]<float>[NS],[-]<float>[EW] with non-optional NSEW
if( !$latlong && preg_match("/^\s*(-?)([\d\.]+)\s*([NS])\s*,?\s*(-?)([\d\.]+)\s*([EW])\s*$/", strtoupper($dms), $latlong ) ) {
$latlong = array('',$latlong[1],$latlong[2],'0','0',$latlong[3],$latlong[4],$latlong[5],'0','0',$latlong[6]);
}
//simple regex ;) ... matches upto 3 sets of number[non-number] per northing/easting (allows for DMS, DM or D)
//uses \D+ to avoid encoding conflicts with ° character which can become multibyte (u flag just fails to match)
//smart-quote versions of ' and " are also multi-byte characters
//it must not treat NSEW as if they are part of the degree symbol with decimal degrees, or other separators
//note that this cannot accept HTML strings from dd_to_dms as it will not match "
//[-]<float><separator>[<float><separator>[<float><separator>]]([NS][,]|,)[-]<float><separator>[<float><separator>[<float><separator>]][EW]
//Captures -, <float>, <float>, <float>, [NS], -, <float>, <float>, <float>, [EW]
if( !$latlong && !preg_match( "/^\s*(-?)([\d\.]+)\D+?\s*(?:([\d\.]+)\D+?\s*(?:([\d\.]+)\D+?\s*)?)?(?:([NS])\s*,?|,)\s*(-?)([\d\.]+)\D+?\s*(?:([\d\.]+)\D+?\s*(?:([\d\.]+)\D+?\s*)?)?([EW]?)\s*$/", strtoupper($dms), $latlong ) ) {
//invalid format
if( $deny_bad_coords ) {
return false;
}
//assume 0,0
$latlong = array('','','0','0','0','N','','0','0','0','E');
}
if( $deny_bad_coords && ( is_nan($latlong[2]) || is_nan($latlong[7]) ) ) {
return false;
}
}
if( !$latlong[3] ) { $latlong[3] = 0; }
if( !$latlong[4] ) { $latlong[4] = 0; }
if( !$latlong[8] ) { $latlong[8] = 0; }
if( !$latlong[9] ) { $latlong[9] = 0; }
$lat = $latlong[2] + $latlong[3] / 60 + $latlong[4] / 3600;
if( $latlong[1] ) { $lat *= -1; }
if( $latlong[5] == 'S' || $latlong[5] == -1 ) { $lat *= -1; }
$long = $latlong[7] + $latlong[8] / 60 + $latlong[9] / 3600;
if( $latlong[6] ) { $long *= -1; }
if( $latlong[10] == 'W' || $latlong[10] == -1 ) { $long *= -1; }
return $this->generic_lat_long_coords( $return_type, $lat, $long );
}
//Helmert transform parameters used during Helmert transformations
//OSGB<->WGS84 parameters taken from "6.6 Approximate WGS84 to OSGB36/ODN transformation"
//http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/guide6.html
private $Helmert_WGS84_to_OSGB36 = array(
'tx' => -446.448,
'ty' => 125.157,
'tz' => -542.060,
's' => 20.4894,
'rx' => -0.1502,
'ry' => -0.2470,
'rz' => -0.8421
);
private $Helmert_OSGB36_to_WGS84 = array(
'tx' => 446.448,
'ty' => -125.157,
'tz' => 542.060,
's' => -20.4894,
'rx' => 0.1502,
'ry' => 0.2470,
'rz' => 0.8421
);
//Ireland65<->WGS84 parameters taken from http://en.wikipedia.org/wiki/Helmert_transformation
private $Helmert_WGS84_to_Ireland65 = array(
'tx' => -482.53,
'ty' => 130.596,
'tz' => -564.557,
's' => -8.15,
'rx' => 1.042,
'ry' => 0.214,
'rz' => 0.631
);
private $Helmert_Ireland65_to_WGS84 = array(
'tx' => 482.53,
'ty' => -130.596,
'tz' => 564.557,
's' => 8.15,
'rx' => -1.042,
'ry' => -0.214,
'rz' => -0.631
);
public function get_transformation($name) {
$name = 'Helmert_'.$name;
if( isset( $this->$name ) ) {
return $this->$name;
}
return false;
}
public function create_transformation( $tx, $ty, $tz, $s, $rx, $ry, $rz ) {
return array('tx'=>$tx,'ty'=>$ty,'tz'=>$tz,'s'=>$s,'rx'=>$rx,'ry'=>$ry,'rz'=>$rz);
}
public function Helmert_transform( $N, $E, $H, $from, $via = false, $to = false, $return_type = false ) {
//conversion according to formulae listed on http://www.movable-type.co.uk/scripts/latlong-convert-coords.html
//parts taken from http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/
list( $N, $E, $H, $from, $via, $to, $return_type ) = $this->expand_params( $N, 2, 3, array( $E, $H, $from, $via, $to, $return_type ) );
if( !$from || !isset($from['a']) || !$via || !isset($via['tx']) || !$to || !isset($to['a']) ) {
return false;
}
$has_height = is_numeric($H);
if( !$has_height ) {
$H = 0;
}
//work in radians
$N = deg2rad($N);
$E = deg2rad($E);
//convert polar coords to cartesian
//eccentricity-squared of source ellipsoid from Annexe B B1
$e2 = 1 - ( $from['b'] * $from['b'] ) / ( $from['a'] * $from['a'] );
$sin_Phi = sin( $N );
$cos_Phi = cos( $N );
//transverse radius of curvature
$Nu = $from['a'] / sqrt( 1 - $e2 * $sin_Phi * $sin_Phi );
$x = ( $Nu + $H ) * $cos_Phi * cos( $E );
$y = ( $Nu + $H ) * $cos_Phi * sin( $E );
$z = ( ( 1 - $e2 ) * $Nu + $H ) * $sin_Phi;
//extracts variables called $tx,$ty,$tz,$s,$rx,$ry,$rz
extract( $via, EXTR_SKIP );
//convert seconds to radians
$rx *= M_PI / 648000;
$ry *= M_PI / 648000;
$rz *= M_PI / 648000;
//convert ppm to pp_one, and add one to avoid recalculating
$s = $s / 1000000 + 1;
//apply Helmert transform (algorithm notes incorrectly show rx instead of rz in $x1 line)
$x1 = $tx + $s * $x - $rz * $y + $ry * $z;
$y1 = $ty + $rz * $x + $s * $y - $rx * $z;
$z1 = $tz - $ry * $x + $rx * $y + $s * $z;
//convert cartesian coords back to polar
//eccentricity-squared of destination ellipsoid from Annexe B B1
$e2 = 1 - ( $to['b'] * $to['b'] ) / ( $to['a'] * $to['a'] );
$p = sqrt( $x1 * $x1 + $y1 * $y1 );
$Phi = atan2( $z1, $p * ( 1 - $e2 ) );
$Phi1 = 2 * M_PI;
$accuracy = 0.000001 / $to['a']; //0.01 mm accuracy, though the OS transform itself only has 4-5 metres
$loopcount = 0;
//due to number precision, it is possible to get infinite loops here for extreme cases (especially for invalid parameters)
//in tests, upto 4 loops are needed - if it reaches 100, assume it must be infinite, and break out
while( abs( $Phi - $Phi1 ) > $accuracy && $loopcount++ < 100 ) {
$sin_Phi = sin( $Phi );
$Nu = $to['a'] / sqrt( 1 - $e2 * $sin_Phi * $sin_Phi );
$Phi1 = $Phi;
$Phi = atan2( $z1 + $e2 * $Nu * $sin_Phi, $p );
}
$Lambda = atan2( $y1, $x1 );
$H = ( $p / cos( $Phi ) ) - $Nu;
//done converting - return in degrees
$latitude = rad2deg($Phi);
$longitude = rad2deg($Lambda);
return $this->generic_lat_long_coords( $return_type, $latitude, $longitude, $has_height ? $H : null );
}
//polynomial transforming lat-long coordinates, such as OSi/LPS (Irish) Polynomial Transformation
private $polycoeff_OSiLPS = array(
'order' => 3,
'positive' => true,
'latm' => 53.5,
'longm' => -7.7,
'k0' => 0.1,
'A' => array(
array(0.763,0.123,0.183,-0.374),
array(-4.487,-0.515,0.414,13.110),
array(0.215,-0.570,5.703,113.743),
array(-0.265,2.852,-61.678,-265.898)
),
'B' => array(
array(-2.810,-4.680,0.170,2.163),
array(-0.341,-0.119,3.913,18.867),
array(1.196,4.877,-27.795,-284.294),
array(-0.887,-46.666,-95.377,-853.950)
)
);
public function get_polynomial_coefficients($name) {
$name = 'polycoeff_'.$name;
if( isset( $this->$name ) ) {
return $this->$name;
}
return false;
}
public function create_polynomial_coefficients( $order, $positive, $latm, $longm, $k0, $A, $B ) {
if( !isset($A) || !is_array($A) || !isset($A[0]) || !is_array($A[0]) || !isset($B) || !is_array($B) || !isset($B[0]) || !is_array($B[0]) ) {
return false;
}
return array( 'order' => $order, 'positive' => $positive, 'latm' => $latm, 'longm' => $longm, 'k0' => $k0, 'A' => $A, 'B' => $B );
}
private function polynomial3( $latitude, $longitude, $latm, $longm, $k0, $A, $B ) {
//transformation according to Transformations and OSGM15 User Guide
//https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/for-developers
$U = $k0 * ( $latitude - $latm );
$V = $k0 * ( $longitude - $longm );
//0 order
$dlat = ( $A[0][0] +
//1st order
$A[1][0] * $U + $A[0][1] * $V + $A[1][1] * $U * $V +
//2nd order
$A[2][0] * $U*$U + $A[0][2] * $V*$V + $A[2][1] * $U*$U * $V + $A[1][2] * $U * $V*$V + $A[2][2] * $U*$U * $V*$V +
//3rd order
$A[3][0] * $U*$U*$U + $A[0][3] * $V*$V*$V + $A[3][1] * $U*$U*$U * $V + $A[1][3] * $U * $V*$V*$V + $A[3][2] * $U*$U*$U * $V*$V + $A[2][3] * $U*$U * $V*$V*$V + $A[3][3] * $U*$U*$U * $V*$V*$V
) / 3600;
//0 order
$dlong = ( $B[0][0] +
//1st order
$B[1][0] * $U + $B[0][1] * $V + $B[1][1] * $U * $V +
//2nd order
$B[2][0] * $U*$U + $B[0][2] * $V*$V + $B[2][1] * $U*$U * $V + $B[1][2] * $U * $V*$V + $B[2][2] * $U*$U * $V*$V +
//3rd order
$B[3][0] * $U*$U*$U + $B[0][3] * $V*$V*$V + $B[3][1] * $U*$U*$U * $V + $B[1][3] * $U * $V*$V*$V + $B[3][2] * $U*$U*$U * $V*$V + $B[2][3] * $U*$U * $V*$V*$V + $B[3][3] * $U*$U*$U * $V*$V*$V
) / 3600;
return array( $dlat, $dlong );
}
public function polynomial_transform( $latitude, $longitude, $coefficients = false, $return_type = false ) {
list( $latitude, $longitude, $coefficients, $return_type ) = $this->expand_params( $latitude, 2, 2, array( $longitude, $coefficients, $return_type ) );
if( is_array($coefficients) && $coefficients['order'] == 3 ) {
//currently, third order polynomial is the only supported type
$shifts = $this->polynomial3( $latitude, $longitude, $coefficients['latm'], $coefficients['longm'], $coefficients['k0'], $coefficients['A'], $coefficients['B'] );
} else {
return false;
}
return $this->generic_lat_long_coords( $return_type, $coefficients['positive'] ? ( $latitude + $shifts[0] ) : ( $latitude - $shifts[0] ), $coefficients['positive'] ? ( $longitude + $shifts[1] ) : ( $longitude - $shifts[1] ) );
}
public function reverse_polynomial_transform( $latitude, $longitude, $coefficients = false, $return_type = false ) {
list( $latitude, $longitude, $coefficients, $return_type ) = $this->expand_params( $latitude, 2, 2, array( $longitude, $coefficients, $return_type ) );
$slat = $slong = 0;
$latshifted = $latitude;
$longshifted = $longitude;
$loopcount = 0;
do {
//get the shifts for the current estimated point, starting with the initial easting and northing
if( is_array($coefficients) && isset($coefficients['order']) && $coefficients['order'] == 3 ) {
//currently, third order polynomial is the only supported type
$shifts = $this->polynomial3( $latshifted, $longshifted, $coefficients['latm'], $coefficients['longm'], $coefficients['k0'], $coefficients['A'], $coefficients['B'] );
} else {
return false;
}
if( !$shifts ) {
return $this->generic_lat_long_coords( $return_type, NAN, NAN );
}
$difflat = $slat - $shifts[0];
$difflong = $slong - $shifts[1];
$slat = $shifts[0];
$slong = $shifts[1];
//apply the new shift and try again with the shift from the new point
$latshifted = $coefficients['positive'] ? ( $latitude - $slat ) : ( $latitude + $slat );
$longshifted = $coefficients['positive'] ? ( $longitude - $slong ) : ( $longitude + $slong );
//see if the shift changes with each iteration
//unlikely to ever hit 100 iterations, even when taken to extremes
} while( ( abs( $difflat ) > 1E-12 || abs( $difflong ) > 1E-12 ) && ++$loopcount < 100 );
return $this->generic_lat_long_coords( $return_type, $latshifted, $longshifted );
}
//geodesics
//Vincenty formulas with precision detection from https://www.movable-type.co.uk/scripts/latlong-vincenty.html and additional error handling
//direct
public function get_geodesic_destination( $latfrom, $longfrom, $geodeticdistance, $bearingfrom = false, $ellipsoid = false, $return_type = false ) {
list( $latfrom, $longfrom, $geodeticdistance, $bearingfrom, $ellipsoid, $return_type ) =
$this->expand_params( $latfrom, 2, 2, array( $longfrom, $geodeticdistance, $bearingfrom, $ellipsoid, $return_type ) );
list( $geodeticdistance, $bearingfrom, $ellipsoid, $return_type ) = $this->expand_params( $geodeticdistance, 2, 2, array( $bearingfrom, $ellipsoid, $return_type ) );
if( !is_array($ellipsoid) || !isset($ellipsoid['a']) || !isset($ellipsoid['b']) ) {
return false;
}
$f = ( $ellipsoid['a'] - $ellipsoid['b'] ) / $ellipsoid['a'];
$longfrom = $this->wrap_around_longitude($longfrom);
$lat1 = deg2rad($latfrom);
$long1 = deg2rad($longfrom);
$bearingfrom = deg2rad($bearingfrom);
//U is called "reduced latitude", latitude on the auxiliary sphere
$tanU1 = ( 1 - $f ) * tan($lat1);
//use trig identities instead of atan then cos or sin
$cosU1 = 1 / sqrt( 1 + pow( $tanU1, 2 ) );
$sinU1 = $tanU1 * $cosU1;
$s1 = atan2( $tanU1, cos($bearingfrom) );
//a = forward azimuth of the geodesic at the equator, if it were extended that far
$sinA = $cosU1 * sin($bearingfrom);
$sin2A = pow( $sinA, 2 );
$cos2A = 1 - $sin2A;
$u2 = $cos2A * ( pow( $ellipsoid['a'], 2 ) - pow( $ellipsoid['b'], 2 ) ) / pow( $ellipsoid['b'], 2 );
$A = 1 + $u2 * ( 4096 + $u2 * ( $u2 * ( 320 - 175 * $u2 ) - 768 ) ) / 16384;
$B = $u2 * ( 256 + $u2 * ( $u2 * ( 74 - 47 * $u2 ) - 128 ) ) / 1024;
//s = angular separation between points on auxiliary sphere
$s = $geodeticdistance / ( $ellipsoid['b'] * $A ); //initial approximation
$loopcount = 0;
do {
//SM = angular separation between the midpoint of the line and the equator
$cos2SM = cos( 2 * $s1 + $s );
$cos22SM = pow( $cos2SM, 2 );
$sin2s = sin($s);
//difference in angular separation between auxiliary sphere and ellipsoid
$deltaS = $B * $sin2s * ( $cos2SM + $B * ( cos($s) * ( 2 * pow( $cos2SM, 2 ) - 1 ) - $B * $cos2SM * ( 4 * pow( $sin2s, 2 ) - 3 ) * ( 4 * $cos22SM - 3 ) / 6 ) / 4 );
$sold = $s;
$s = $geodeticdistance / ( $ellipsoid['b'] * $A ) + $deltaS;
} while( abs( $s - $sold ) < 1E-12 && ++$loopcount < 500 );
$sins = sin($s);
$coss = cos($s);
$sinbearing = sin($bearingfrom);
$cosbearing = cos($bearingfrom);
$lat2 = atan2( $sinU1 * $coss + $cosU1 * $sins * $cosbearing, ( 1 - $f ) * sqrt( $sin2A + pow( $sinU1 * $sins - $cosU1 * $coss * $cosbearing, 2 ) ) );
//difference in longitude of the points on the auxiliary sphere
$lambda = atan2( $sins * $sinbearing, $cosU1 * $coss - $sinU1 * $sins * $cosbearing );
$C = $f * $cos2A * ( 4 + $f * ( 4 - 3 * $cos2A ) ) / 16;
//longitudinal difference of the points on the ellipsoid
$L = $lambda - ( 1 - $C) * $f * $sinA * ( $s + $C * $sins * ( $cos2SM + $C * $coss * ( 2 * $cos22SM - 1 ) ) );
$long2 = $long1 + $L;
$longto = $this->wrap_around_longitude(rad2deg($long2));
$latto = rad2deg($lat2);
$bearing2 = atan2( $sinA, -1 * ( $sinU1 * $sins - $cosU1 * $coss * $cosbearing ) );
$bearingto = $this->wrap_around_azimuth(rad2deg($bearing2));
if( $return_type ) {
//sprintf to produce simple numbers instead of scientific notation
$bearingto = sprintf( '%F', $bearingto );
//bearing can be 360, after rounding
if( $bearingto == '360.000000' ) {
$bearingto = '0.000000';
}
$deg = is_string($return_type) ? $return_type : '°';
return $this->generic_lat_long_coords( $return_type, $latto, $longto ) . ' @ ' . $bearingto . $deg ;
} else {
//avoid stupid -0 by adding 0
return array( $latto + 0, $longto + 0, $bearingto + 0 );
}
}
//inverse
public function get_geodesic( $latfrom, $longfrom, $latto, $longto = false, $ellipsoid = false, $return_type = false ) {
list( $latfrom, $longfrom, $latto, $longto, $ellipsoid, $return_type ) =
$this->expand_params( $latfrom, 2, 2, array( $longfrom, $latto, $longto, $ellipsoid, $return_type ) );
list( $latto, $longto, $ellipsoid, $return_type ) = $this->expand_params( $latto, 2, 2, array( $longto, $ellipsoid, $return_type ) );
if( !is_array($ellipsoid) || !isset($ellipsoid['a']) || !isset($ellipsoid['b']) ) {
return false;
}
$f = ( $ellipsoid['a'] - $ellipsoid['b'] ) / $ellipsoid['a'];
$longfrom = $this->wrap_around_longitude($longfrom);
$longto = $this->wrap_around_longitude($longto);
$lat1 = deg2rad($latfrom);
$long1 = deg2rad($longfrom);
$lat2 = deg2rad($latto);
$long2 = deg2rad($longto);
$L = $long2 - $long1; //longitudinal difference
//U is called "reduced latitude", latitude on the auxiliary sphere
$tanU1 = ( 1 - $f ) * tan($lat1);
$tanU2 = ( 1 - $f ) * tan($lat2);
//use trig identities instead of atan then cos or sin
$cosU1 = 1 / sqrt( 1 + pow( $tanU1, 2 ) );
$cosU2 = 1 / sqrt( 1 + pow( $tanU2, 2 ) );
$sinU1 = $tanU1 * $cosU1;
$sinU2 = $tanU2 * $cosU2;
$lambda = $L; //difference in longitude on an auxiliary sphere, initial approximation
//start detecting failures, not part of Vincenty formulas
$detectedproblem = false;
//~1 iteration for spheres, ~3-4 for spheroids
//extreme spheroids need more; a=4b can need 80+
//allow max 500 iterations, or bearing results flip around when resolution fails
$maxloops = 500;
//check if the points are somewhat far from each other, for use when handling precision errors
//either the latitudes are more than half a circle different, so that the points are near opposite poles
//or the latitudes are nearly opposite and longitudes are similarly different (points far apart near-ish the equator)
$antipodal = abs( $lat2 - $lat1 ) > M_PI / 2 || ( abs( $lat1 + $lat2 ) < M_PI / 2 && cos($L) <= 0 );
$prolate = $ellipsoid['b'] > $ellipsoid['a'];
do {
$sinS = sqrt( pow( $cosU2 * sin($lambda), 2 ) + pow( $cosU1 * $sinU2 - $sinU1 * $cosU2 * cos($lambda), 2 ) );
//prolate spheroids often manage to resolve in spite of having a tiny sinS
//their antipodal distances are much harder because they do not travel over the pole, and need the full calculation
if( abs( $sinS ) < PHP_FLOAT_EPSILON && ( !$prolate || !$antipodal ) ) {
//not part of Vincenty formulas
//when the sine is smaller than the number precision, the iteration can give wildly inaccurate results, happens when:
//* the points are almost opposite each other; measure the arc via the north pole
//* the points are very close to each other, length becomes 0
//* the points are the same, length is 0
$detectedproblem = true;
$cos2A = 1; //how much of the transit is via the minor axis (smaller numbers use the major axis dimensions)
$cos2SM = $sinS = $cosS = 0; //these will not be used
$s = $antipodal ? M_PI : 0; //distance around the spheroid
break;
}
$cosS = $sinU1 * $sinU2 + $cosU1 * $cosU2 * cos($lambda);
$s = atan2( $sinS, $cosS ); //angular separation between points
//a = forward azimuth of the geodesic at the equator, if it were extended that far
$sinA = $cosU1 * $cosU2 * sin($lambda) / $sinS;
$cos2A = 1 - pow( $sinA, 2 );
//SM is angular separation between the midpoint of the line and the equator
//at an equatorial line, cos2A becomes 0, and dividing by 0 is invalid, so set cos2SM to 0 instead
$cos2SM = $cos2A ? ( $cosS - 2 * $sinU1 * $sinU2 / $cos2A ) : 0;
$C = $f * $cos2A * ( 4 + $f * ( 4 - 3 * $cos2A ) ) / 16;
$lambdaold = $lambda;
$lambda = $L + ( 1 - $C ) * $f * $sinA * ( $s + $C * $sinS * ( $cos2SM + $C * $cosS * ( -1 + 2 * pow( $cos2SM, 2 ) ) ) );
//lambda > PI when it fails to resolve, but also sometimes when it works - rely only on loop counting
} while( abs( $lambda - $lambdaold ) > 1E-12 && --$maxloops );
$detectedproblem = $detectedproblem || !$maxloops;
$u2 = $cos2A * ( pow( $ellipsoid['a'], 2 ) - pow( $ellipsoid['b'], 2 ) ) / pow( $ellipsoid['b'], 2 );
$A = 1 + $u2 * ( 4096 + $u2 * ( $u2 * ( 320 - 175 * $u2 ) - 768 ) ) / 16384;
$B = $u2 * ( 256 + $u2 * ( $u2 * ( 74 - 47 * $u2 ) - 128 ) ) / 1024;
//difference in angular separation between auxiliary sphere and ellipsoid
$deltaS = $B * $sinS * ( $cos2SM + $B * ( $cosS * ( 2 * $cos2SM - 1 ) - $B * $cos2SM * ( 4 * pow( $sinS, 2 ) - 3 ) * ( 4 * $cos2SM - 3 ) / 6 ) / 4 );
$geodeticdistance = $ellipsoid['b'] * $A * ( $s - $deltaS );
$problemdata = null;
if( $detectedproblem ) {
//error handling to try to get useful bearings - not part of the Vincenty formulas
$problemdata = array( 'detectiontype' => $maxloops ? 'BELOW_RESOLUTION' : 'DID_NOT_CONVERGE' );
if( $antipodal ) {
if( abs( $latfrom - $latto ) == 180 || ( $latfrom + $latto == 0 && abs( $longto - $longfrom ) == 180 ) ) {
$problemdata['separation'] = 'ANTIPODAL';
} else {
$problemdata['separation'] = 'NEARLY_ANTIPODAL';
}
} else {
$problemdata['separation'] = 'NEARLY_IDENTICAL';
}
if( $latfrom == $latto && ( $longfrom == $longto || $latfrom == 90 || $latfrom == -90 ) ) {
//identical points, nothing useful can be returned
$bearingfrom = $bearingto = 0;
$problemdata['distanceconfidence'] = 1;
$problemdata['azimuthconfidence'] = 1;
$problemdata['separation'] = 'IDENTICAL';
} elseif( $latfrom == 90 ) {
//it seems to always converge when one of the points is on a pole, but these fixes are provided just in case they are needed
//with Vincenty formulae, azimuth at 90,* is measured with 180 lined up with the longitude, like it would at 0,<longitude>
$bearingfrom = M_PI + $long1 - $long2;
//it comes down the longto meridian, so it will always be seen as 180, even when latto is -90
$bearingto = M_PI;
$problemdata['distanceconfidence'] = ( $latto == -90 ) ? 1 : 0.99;
$problemdata['azimuthconfidence'] = 1;
} elseif( $latfrom == -90 ) {
//azimuth at -90,* is measured with 0 lined up with the longitude
$bearingfrom = $long2 - $long1;
$bearingto = 0;
$problemdata['distanceconfidence'] = ( $latto == 90 ) ? 1 : 0.99;
$problemdata['azimuthconfidence'] = 1;
} elseif( $latto == 90 ) {
//always north up their own meridian
$bearingfrom = 0;
$bearingto = $long2 - $long1;
$problemdata['distanceconfidence'] = ( $latfrom == -90 ) ? 1 : 0.99;
$problemdata['azimuthconfidence'] = 1;
} elseif( $latto == -90 ) {
//always south
$bearingfrom = M_PI;
$bearingto = M_PI + $long1 - $long2;
$problemdata['distanceconfidence'] = ( $latfrom == 90 ) ? 1 : 0.99;
$problemdata['azimuthconfidence'] = 1;
} elseif( $antipodal ) {
$problemdata['distanceconfidence'] = ( $problemdata['detectiontype'] == 'DID_NOT_CONVERGE' ) ? 0.6 : 0.98;
if( $prolate ) {
//prolate spheroid, the direction will be east or west
if( sin( $long2 - $long1 ) < 0 ) {
$bearingfrom = M_PI / -2;
$bearingto = M_PI / -2;
} else {
$bearingfrom = M_PI / 2;
$bearingto = M_PI / 2;
}
$problemdata['azimuthconfidence'] = 0.9;
} else {
//with antipodal points, always travel via the poles; choose the closest one if the points are not quite identical, or north otherwise
if( $lat1 + $lat2 < 0 ) {
$bearingfrom = M_PI;
$bearingto = 0;
} else {
$bearingfrom = 0;
$bearingto = M_PI;
}
$problemdata['azimuthconfidence'] = ( abs( $longto - $longfrom ) == 180 ) ? 1 : 0.9;
}
//in this case, the length is left with whatever it was
} else {
//points close together; below number resolution required for trig functions
//can be calculated using basic geometry near the equator, but this is nonsense near the poles, back to being the original problem
//this cannot be calculated in any reliable way, but this approximation is better than nothing over much of the ellipsoid
$north = $latto > $latfrom;
$south = $latfrom > $latto;
$east = $longto > $longfrom || $longfrom - $longto >= 180;
$west = $longfrom > $longto || $longto - $longfrom > 180;
if( $north && $east ) {
$bearingfrom = $bearingto = M_PI / 4;
} elseif( !$north && !$south && $east ) {
$bearingfrom = $bearingto = M_PI / 2;
} elseif( $south && $east ) {
$bearingfrom = $bearingto = 3 * M_PI / 4;
} elseif( $south && !$east && !$west ) {
$bearingfrom = $bearingto = M_PI;
} elseif( $south && $west ) {
$bearingfrom = $bearingto = 5 * M_PI / 4;
} elseif( !$north && !$south && $west ) {
$bearingfrom = $bearingto = 3 * M_PI / 2;
} elseif( $north && $west ) {
$bearingfrom = $bearingto = 7 * M_PI / 4;
} else {
$bearingfrom = $bearingto = 0;
}
$problemdata['distanceconfidence'] = 0.98;
$problemdata['azimuthconfidence'] = 0.5 * cos($lat1);
}
} else {
//actual measurement
$bearingfrom = atan2( $cosU2 * sin($lambda), ( $cosU1 * $sinU2 ) - ( $sinU1 * $cosU2 * cos($lambda) ) );
$bearingto = atan2( $cosU1 * sin($lambda), ( $cosU1 * $sinU2 * cos($lambda) ) - ( $sinU1 * $cosU2 ) );
}
//Vincenty formula returns bearings of -180 to 180, and error handling can make numbers greater than expected bounds
$bearingfrom = $this->wrap_around_azimuth(rad2deg($bearingfrom));
$bearingto = $this->wrap_around_azimuth(rad2deg($bearingto));
if( $return_type ) {
//sprintf to produce simple numbers instead of scientific notation
$bearingfrom = sprintf( '%F', $bearingfrom );
$bearingto = sprintf( '%F', $bearingto );
//bearing can be 360, after rounding
if( $bearingfrom == '360.000000' ) {
$bearingfrom = '0.000000';
}
if( $bearingto == '360.000000' ) {
$bearingto = '0.000000';
}
$deg = is_string($return_type) ? $return_type : '°';
$angle = is_string($return_type) ? '>' : '>';
return ( ( $detectedproblem || $lambda > M_PI || !$maxloops ) ? '~ ' : '' ) . sprintf( '%F', $geodeticdistance ) . ' m, ' . $bearingfrom . $deg . ' =' . $angle . ' ' . $bearingto . $deg ;
} else {
//avoid stupid -0 by adding 0
return array( $geodeticdistance, $bearingfrom + 0, $bearingto + 0, $problemdata );
}
}
//lat-long based geoid grid
private $geoidcache = array();
private $geoidgrd_OSGM15_Belfast = array(
'name' => 'OSGM15_Belfast',
'latmin' => 51.000000,
'latmax' => 56.000000,
'longmin' => -11.500000,
'longmax' => -5.000000,
'latspacing' => 0.013333,
'longspacing' => 0.020000,
'latstepdp' => 6,
'longstepdp' => 6,
);
private $geoidgrd_OSGM15_Malin = array(
'name' => 'OSGM15_Malin',
'latmin' => 51.000000,
'latmax' => 56.000000,
'longmin' => -11.500000,
'longmax' => -5.000000,
'latspacing' => 0.013333,
'longspacing' => 0.020000,
'latstepdp' => 6,
'longstepdp' => 6,
);
private $geoidgrd_EGM96_ww15mgh = array(
'name' => 'EGM96_ww15mgh',
'latmin' => -90.000000,
'latmax' => 90.000000,
'longmin' => .000000,
'longmax' => 360.000000,
'latspacing' => .250000,
'longspacing' => .250000,
'latstepdp' => 6,
'longstepdp' => 6,
);
public function get_geoid_grd($name) {
$name = 'geoidgrd_'.$name;
if( isset( $this->$name ) ) {
return $this->$name;
}
return false;
}
public function create_geoid_grd( $name, $latmin, $latmax, $longmin, $longmax, $latspacing, $longspacing, $latstepdp, $longstepdp ) {
if( !is_string($name) || !strlen($name) || !( $latspacing * 1 ) || !( $longspacing * 1 ) || $latspacing < 0 || $longspacing < 0 || $latmax < $latmin || $longmax < $longmin ) {
return false;
}
return array(
'name' => $name,
'latmin' => $latmin,
'latmax' => $latmax,
'longmin' => $longmin,
'longmax' => $longmax,
'latspacing' => $latspacing,
'longspacing' => $longspacing,
'latstepdp' => $latstepdp,
'longstepdp' => $longstepdp
);
}
public function create_geoid_record( $recordindex, $geoid_undulation ) {
//this does not actually get used in this format - it exists only to force the right columns to exist
return array( 'name' => "r_$recordindex", 'height' => $geoid_undulation );
}
public function cache_geoid_records( $name, $records ) {
if( !is_string($name) || !strlen($name) ) {
return false;
}
if( !isset($this->geoidcache[$name]) || !$this->geoidcache[$name] ) {
//nothing has ever been cached for this geoid set, create the cache
$this->geoidcache[$name] = array();
}
if( !$records ) {
return true;
}
//allow it to accept either a single record, or an array of records
if( !is_array($records) || isset($records['name']) ) {
$records = array($records);
}
foreach( $records as $onerecord ) {
if( !is_array($onerecord) || !isset($onerecord['name']) || !isset($onerecord['height']) ) {
return false;
}
$this->geoidcache[$name][$onerecord['name']] = $onerecord['height'];
}
return true;
}
public function delete_geoid_cache( $name = '' ) {
if( strlen($name) ) {
unset($this->geoidcache[$name]);
} else {
$this->geoidcache = array();
}
}
private function get_geoid_records( $geoidgrd, $records, $fetch_records ) {
if( !isset($this->geoidcache[$geoidgrd['name']]) || !$this->geoidcache[$geoidgrd['name']] ) {
//nothing has ever been cached for this geoid set, create the cache
$this->geoidcache[$geoidgrd['name']] = array();
}
$unknown_records = array();
$recordnames = array();
foreach( $records as $onerecord ) {
$recordnames[] = $id = 'r_'.$onerecord['recordindex'];
if( !isset($this->geoidcache[$geoidgrd['name']][$id]) ) {
$onerecord['latitude'] = round( $onerecord['latitude'], $geoidgrd['latstepdp'] );
$onerecord['longitude'] = round( $onerecord['longitude'], $geoidgrd['longstepdp'] );
$unknown_records[] = $onerecord;
$this->geoidcache[$geoidgrd['name']][$id] = false; //false gets detected as "set", null does not
}
}
if( count($unknown_records) ) {
$this->cache_geoid_records( $geoidgrd['name'], call_user_func( $fetch_records, $geoidgrd['name'], $unknown_records ) );
}
$shifts = array();
foreach( $recordnames as $id ) {
$shifts[] = $this->geoidcache[$geoidgrd['name']][$id];
}
return $shifts;
}
public function apply_geoid( $latitude, $longitude, $height, $direction, $geoidgrd, $fetch_records = false, $return_type = false, $deny_bad_coords = false ) {
list( $latitude, $longitude, $height, $direction, $geoidgrd, $fetch_records, $return_type, $deny_bad_coords ) =
$this->expand_params( $latitude, 2, 3, array( $longitude, $height, $direction, $geoidgrd, $fetch_records, $return_type, $deny_bad_coords ) );
if( !$geoidgrd || !( $geoidgrd['latspacing'] * $geoidgrd['longspacing'] ) || $geoidgrd['latspacing'] < 0 || $geoidgrd['longspacing'] < 0 || $geoidgrd['latmax'] < $geoidgrd['latmin'] || $geoidgrd['longmax'] < $geoidgrd['longmin'] || !is_callable($fetch_records) ) {
return false;
}
$longitude = $this->wrap_around_longitude($longitude);
//OSI data has a rounding precision problem that causes it to miss its end point; 0.013333 is not the same as 5 / 375
$latrange = $geoidgrd['latmax'] - $geoidgrd['latmin'];
$latspacing = $latrange ? ( $latrange / round( $latrange / $geoidgrd['latspacing'] ) ) : $geoidgrd['latspacing'];
$longrange = $geoidgrd['longmax'] - $geoidgrd['longmin'];
$longperrow = round( $longrange / $geoidgrd['longspacing'] );
$longspacing = $longrange ? ( $longrange / $longperrow ) : $geoidgrd['longspacing'];
$longperrow++; //fencepost problem
if( $longitude < 0 && $longitude < $geoidgrd['longmin'] && $geoidgrd['longmax'] > 180 ) {
//geoid data may use -180 <= x <= 180, or it may use 0 <= x <= 360
//this also allows them to use odd ranges like -10 to 200
$longitude += 360;
}
if( $latitude < $geoidgrd['latmin'] || $latitude > $geoidgrd['latmax'] || $longitude < $geoidgrd['longmin'] || $longitude > $geoidgrd['longmax'] ) {
//the point is not within range, so this is not valid
if( $deny_bad_coords ) {
return false;
} else {
return $return_type ? 'NaN' : NAN;
}
}
//grid cell, latitude indexes are from the top, not the bottom
$lat_index = ceil( ( $geoidgrd['latmax'] - $latitude ) / $latspacing );
$long_index = floor( ( $longitude - $geoidgrd['longmin'] ) / $longspacing );
//lat-long of southwest corner
$x0 = $geoidgrd['longmin'] + $long_index * $longspacing;
$y0 = $geoidgrd['latmax'] - $lat_index * $latspacing;
//offset within grid
$dx = $longitude - $x0;
$dy = $latitude - $y0;
$records = array();
$records[0] = array( 'latitude' => $y0, 'longitude' => $x0, 'latindex' => $lat_index, 'longindex' => $long_index, 'recordindex' => $longperrow * $lat_index + $long_index );
//allow points to sit on the last line by setting them to the same values as each other (otherwise it could never work on the north pole)
if( $dx ) {
$records[1] = array( 'latitude' => $y0, 'longitude' => ( $long_index + 1 ) * $longspacing + $geoidgrd['longmin'], 'latindex' => $lat_index, 'longindex' => $long_index + 1 );
$records[1]['recordindex'] = $longperrow * $records[1]['latindex'] + $records[1]['longindex'];
} else {
$records[1] = $records[0];
}
if( $dy ) {
$records[2] = array( 'latitude' => $geoidgrd['latmax'] - ( $lat_index - 1 ) * $latspacing, 'longitude' => $records[1]['longitude'], 'latindex' => $lat_index - 1, 'longindex' => $dx ? ( $long_index + 1 ) : $long_index );
$records[3] = array( 'latitude' => $records[2]['latitude'], 'longitude' => $records[0]['longitude'], 'latindex' => $lat_index - 1, 'longindex' => $long_index );
$records[2]['recordindex'] = $longperrow * $records[2]['latindex'] + $records[2]['longindex'];
$records[3]['recordindex'] = $longperrow * $records[3]['latindex'] + $records[3]['longindex'];
} else {
$records[2] = $records[1];
$records[3] = $records[0];
}
//fetch records if needed
$shifts = $this->get_geoid_records( $geoidgrd, $records, $fetch_records );
if( !is_numeric($shifts[0]) || !is_numeric($shifts[1]) || !is_numeric($shifts[2]) || !is_numeric($shifts[3]) ) {
//this is within the grid, but data was not available
if( $deny_bad_coords ) {
return false;
} else {
return $return_type ? 'NaN' : NAN;
}
}
//offset values
$t = $dx / $longspacing;
$u = $dy / $latspacing;
//bilinear interpolated shifts
$sg = ( ( 1 - $t ) * ( 1 - $u ) * $shifts[0] ) + ( $t * ( 1 - $u ) * $shifts[1] ) + ( $t * $u * $shifts[2] ) + ( ( 1 - $t ) * $u * $shifts[3] );
$height = $direction ? ( $height - $sg ) : ( $height + $sg );
if( $return_type ) {
//sprintf to produce simple numbers instead of scientific notation
return sprintf( '%F', $height );
} else {
return $height;
}
}
//conversions between latitude/longitude and cartesian coordinates
function xyz_to_lat_long( $X, $Y, $Z = false, $ellipsoid = false, $return_type = false) {
list( $X, $Y, $Z, $ellipsoid, $return_type ) = $this->expand_params( $X, 3, 3, array( $Y, $Z, $ellipsoid, $return_type ) );
//conversion according to "A guide to coordinate systems in Great Britain" Annexe B:
//https://www.ordnancesurvey.co.uk/documents/resources/guide-coordinate-systems-great-britain.pdf
if( !is_array($ellipsoid) || !isset($ellipsoid['a']) || !isset($ellipsoid['b']) ) {
return false;
}
$e2 = 1 - pow( $ellipsoid['b'], 2 ) / pow( $ellipsoid['a'], 2 );
$lon = atan2( $Y * 1, $X * 1 );
$p = sqrt( pow( $X, 2 ) + pow( $Y, 2 ) );
$lat = atan2( $Z, $p * ( 1 - $e2 ) );
$loopcount = 0;
do {
$v = $ellipsoid['a'] / sqrt( 1 - ( $e2 * pow( sin($lat), 2 ) ) );
//https://gssc.esa.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion alternative
//height misses one iteration but the result is the same when diff is 0
//$height = ( $p / cos($lat) ) - $v;
//$lat = atan2( $Z, $p * ( 1 - ( $e2 * $v / ( $v + $height ) ) ) );
$oldlat = $lat;
$lat = atan2( $Z + ( $e2 * $v * sin( $lat ) ), $p );
//due to number precision, it is possible to get infinite loops here for extreme cases, when the difference is enough to roll over a digit then roll it back on the next iteration
//in tests, upto 10 loops are needed for extremes - if it reaches 100, assume it must be infinite, and break out
} while( abs( $lat - $oldlat ) > 0 && ++$loopcount < 100 );
$height = ( $p / cos($lat) ) - $v;
return $this->generic_lat_long_coords( $return_type, rad2deg($lat), rad2deg($lon), $height );
}
public function lat_long_to_xyz( $latitude, $longitude, $height = false, $ellipsoid = false, $return_type = false ) {
list( $latitude, $longitude, $height, $ellipsoid, $return_type ) = $this->expand_params( $latitude, 2, 3, array( $longitude, $height, $ellipsoid, $return_type ) );
if( !is_array($ellipsoid) || !isset($ellipsoid['a']) || !isset($ellipsoid['b']) ) {
return false;
}
$height = is_numeric($height) ? $height : 0;
$latitude = deg2rad($latitude);
$longitude = deg2rad($longitude);
$e2 = 1 - pow( $ellipsoid['b'], 2 ) / pow( $ellipsoid['a'], 2 );
$v = $ellipsoid['a'] / sqrt( 1 - ( $e2 * pow( sin($latitude), 2 ) ) );
$X = ( $v + $height ) * cos($latitude) * cos($longitude);
$Y = ( $v + $height ) * cos($latitude) * sin($longitude);
$Z = ( ( ( 1 - $e2 ) * $v ) + $height ) * sin($latitude);
if( $return_type ) {
return sprintf( '%F', $X ) . ' ' . sprintf( '%F', $Y ) . ' ' . sprintf( '%F', $Z );
} else {
return array( $X, $Y, $Z );
}
}
}
?>