Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
map.on('click', function(event){
var x1_wgs;
var y1_wgs;
async function ajax(x,y){
const data = $.ajax({
url:'/geocoding/egsa_2_wgs.php?x='+x+'&y='+y,
async:true,
type:'GET',
dataType: 'JSON',
contentType: "application/json; charset=utf-8",
success: function (data) {
var x1_wgs = data['0'];
var y1_wgs = data['1'];
}
});
console.log(x1_wgs,y1_wgs);
}
if(startPoint.getGeometry() == null){
startPoint.setGeometry(new ol.geom.Point(event.coordinate));
var x = event.coordinate[0];
var y = event.coordinate[1];
ajax(x,y);
}
<?php
class EGSA_2_WGS{
public function egsa87xy_to_wgs84xy($x,$y,$z)
{
return(array ($x-199.723, $y+74.03, $z+246.018));
}
public function flh2xyz($ellipsoid,$phi,$lambda,$h)
//convert phi, lambda, h to x,y,z - h height above ellipsoid surface
//EGSA87 and HTRS07 ellipsoids are the same and differ only by a fraction of a millimeter to WGS84 ellipsoid in small semi-axis
{
$phi=deg2rad($phi);
$lambda=deg2rad($lambda);
$a=6378137; //equatorial radius
if($ellipsoid=="WGS84") $f = 1/298.257223563;
else if ($ellipsoid=="EGSA87" || $ellipsoid=="HTRS07") $f = 1/298.257222100882711243;
$e2=2*$f-$f*$f;
//$epsilon=$e2/(1-$e2);
$greekni=$a/sqrt(1-$e2*sin($phi)*sin($phi)); //prime vertical radius of curvature at latitude φ
$X=($greekni+$h)*cos($phi)*cos($lambda);
$Y=($greekni+$h)*cos($phi)*sin($lambda);
$Z=((1-$e2)*$greekni+$h)*sin($phi);
//echo "\nIn function:".__FUNCTION__."\n";
//print_r( get_defined_vars());
return(array($X,$Y,$Z));
}
public function to_fl_from_proj_xy($type, $east,$north)
//given x,y (east,north in greek UTM projection this function calulates the f,l on greek ellipsoid (GRS80)
//the $type can be either EGSA87 or HTRS07 which differ only in the False Northing value
{
$a=6378137; //equatorial radius
$f=1/298.257222100882711;
//$a=6377563.396; //airy
//$f=1/299.32496; //Airy
$e=sqrt(2*$f-$f*$f);
$l_orig=24*M_PI/180; //Longitude of natural origin
$f_orig=0*M_PI/180; //Latitude of natural origin
$k_orig=0.9996; //Scale factor at natural origin
$FE=500000; //False easting
if ($type=="EGSA87") $FN=0; //False northing
else if($type=="HTRS07") $FN=-2000000;
//$b=6356752.314140347; //grs80: polar radius
//$f=($a-$b)/$a;
$n=$f/(2-$f);
$B=$a/(1+$n)*(1+$n*$n/4+$n*$n*$n*$n/64);
$h1=$n/2-2/3*$n*$n+5/16*$n*$n*$n+41/180*$n*$n*$n*$n;
$h2=13/48*$n*$n-3/5*$n*$n*$n+557/1440*$n*$n*$n*$n;
$h3=61/240*$n*$n*$n-103/140*$n*$n*$n*$n;
$h4=49561/161280*$n*$n*$n*$n;
//meridional arc distance from equator to projection origin (Mo)
$Qo=asinh(tan($f_orig))-$e*atanh($e * sin($f_orig));
$bo=atan(sinh($Qo));
$ksio0=asin(sin($bo));
$ksio1=$h1*sin(2*$ksio0);
$ksio2=$h2*sin(4*$ksio0);
$ksio3=$h3*sin(6*$ksio0);
$ksio4=$h4*sin(8*$ksio0);
$ksio=$ksio0+$ksio1+$ksio2+$ksio3+$ksio4;
$Mo=$B*$ksio;
//up to here the same as toxy
$h1t=$n/2-2/3*$n*$n+37/96*$n*$n*$n-1/360*$n*$n*$n*$n;
$h2t=1/48*$n*$n+1/15*$n*$n*$n-437/1440*$n*$n*$n*$n;
$h3t=17/480*$n*$n*$n-37/840*$n*$n*$n*$n;
$h4t=4397/161280*$n*$n*$n*$n;
$itat=($east-$FE)/($B*$k_orig);
$ksit=($north-$FN+$k_orig*$Mo)/($B*$k_orig);
$ksi1t=$h1t*sin(2*$ksit)*cosh(2*$itat);
$ita1t=$h1t*cos(2*$ksit)*sinh(2*$itat);
$ksi2t=$h2t*sin(4*$ksit)*cosh(4*$itat);
$ita2t=$h2t*cos(4*$ksit)*sinh(4*$itat);
$ksi3t=$h3t*sin(6*$ksit)*cosh(6*$itat);
$ita3t=$h3t*cos(6*$ksit)*sinh(6*$itat);
$ksi4t=$h4t*sin(8*$ksit)*cosh(8*$itat);
$ita4t=$h4t*cos(8*$ksit)*sinh(8*$itat);
$ksi0t=$ksit-($ksi1t+$ksi2t+$ksi3t+$ksi4t);
$ita0t=$itat-($ita1t+$ita2t+$ita3t+$ita4t);
$vitat=asin(sin($ksi0t)/cosh($ita0t));
$Qt=asinh(tan($vitat));
$Qtt=$Qt+$e*atanh($e*tanh($Qt));
do {
$Qtt_previous=$Qtt;
$Qtt=$Qt+$e*atanh($e*tanh($Qtt));
} while(abs($Qtt-$Qtt_previous)>1e-12);
$phi=atan(sinh($Qtt));
$lambda=$l_orig+asin(tanh($ita0t)/cos($vitat));
return(array(rad2deg($phi),rad2deg($lambda)));
}
public function xyz2flh($ellipsoid,$xo,$yo,$zo)
//convert x,y,z to phi, lambda, h
//EGSA87 and HTRS07 ellipsoids are the same and differ only by a fraction of a millimeter to WGS84 ellipsoid in small semi-axis
{
$a=6378137; //equatorial radius
if($ellipsoid=="WGS84") $f = 1/298.257223563;
else if ($ellipsoid=="EGSA87" || $ellipsoid=="HTRS07") $f = 1/298.257222100882711243;
$e2=2*$f-$f*$f;
$epsilon=$e2/(1-$e2);
$bi=$a*(1-$f);
$p=sqrt($xo*$xo+$yo*$yo);
$q=atan($zo*$a/($p*$bi));
$phi=atan(($zo+$epsilon*$bi*sin($q)*sin($q)*sin($q))/($p-$e2*$a*cos($q)*cos($q)*cos($q)));
$greekni=$a/sqrt(1-$e2*sin($phi)*sin($phi)); //prime vertical radius of curvature at latitude φ
$lambda=atan($yo/$xo);
$h=$p/cos($phi)-$greekni;
//echo "\nIn function:".__FUNCTION__."\n";
//print_r( get_defined_vars());
return(array(rad2deg($phi),rad2deg($lambda),$h));
}
public function EGSA87_to_WGS84($xi,$psi) {
//converts x,y (easting, northing) of EGSA87 projection to φ,λ of WGS84
$temparr=$this->to_fl_from_proj_xy("EGSA87", $xi,$psi);
$temparr=$this->flh2xyz("EGSA87",$temparr[0],$temparr[1],0);
$temparr=$this->egsa87xy_to_wgs84xy($temparr[0],$temparr[1],$temparr[2]);
$temparr=$this->xyz2flh("WGS84",$temparr[0],$temparr[1],$temparr[2]);
// print_r($temparr);
return array($temparr[0],$temparr[1]);
//return array($temparr[0]);
}
}
$x = $_GET['x'];
$y = $_GET['y'];
$trans = new EGSA_2_WGS();
$returnValue1 = json_encode($trans->EGSA87_to_WGS84($x,$y));
header('Content-Type: application/json');
echo $returnValue1;
?>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment