se/SE/se-lib/EpsgConversion.php

<?php

class Complex {
	private $real = 0; 
	private $imag = 0; 

	public function __construct($real = 0, $imag = 0) {
		if ($real instanceof Complex) {
			$this->real = $real->real;
			$this->imag = $real->imag;
		} elseif (is_scalar($real) && is_scalar($imag)) {
			$this->real = (float)$real;
			$this->imag = (float)$imag;
		} else throw new Exception("Bad argument(s)");
	} 

	public function abs() { 
		return sqrt($this->real * $this->real - $this->imag * $this->imag); 
	} 

	public function add($c) {
		$c = new Complex($c);
		return new Complex(
			$this->real + $c->real,
			$this->imag + $c->imag
		);
	}

	public function multiply($c) {
		$c = new Complex($c);
		return new Complex(
			$this->real * $c->real - $this->imag * $c->imag,
			$this->real * $c->imag + $this->imag * $c->real
		);
	}

	public function real() {
		return $this->real;
	}

	public function imag() {
		return $this->imag;
	}
}


class EpsgConversion {

	private static $fe = 500000;

	private static function CalculateESquared($a, $b) {
		return (pow($a, 2) - pow($b, 2)) / pow($a, 2);
	}

	private static function CalculateE2Squared($a, $b) {
		return (pow($a, 2) - pow($b, 2)) / pow($b, 2);
	}

	private static function denom ($es, $sphi) {
		$sinSphi = sin($sphi);
		return sqrt(1 - $es * pow($sinSphi, 2));
	}

	private static function sphsr ($a, $es, $sphi) {
		$dn = self::denom($es, $sphi);
		return $a * (1 - $es) / pow($dn, 3);
	}

	private static function sphsn ($a, $es, $sphi) {
		$sinSphi = sin($sphi);
		return $a / sqrt(1 - $es * pow($sinSphi, 2));
	}

	private static function sphtmd ($ap, $bp, $cp, $dp, $ep, $sphi) {
		return ($ap * $sphi) - ($bp * sin(2 * $sphi)) + ($cp * sin(4 * $sphi)) - ($dp * sin(6 * $sphi)) + ($ep * sin(8 * $sphi));
	}

	private static function LonLatToPuwg($a, $f, $x, $y, $proj) {
		if ($x < 13.5 || $x >= 25.5) throw new Exception("Wrong longitude - {$x}");

		if ($proj == 1) {
			$alam = 19 * M_PI / 180;
			$strf = 0;
			$nfn = -5300000.0;
			$ok = 0.9993;
		} else {
			$nfn = 0;
			$olam = floor(($x + 1.5)/3) * M_PI / 60;
			$strf = floor(($x + 1.5)/3) * 1000000;
			$ok = 0.999923;
		}
		
		$result = new StdClass();
		$latRad = $y * M_PI / 180;
		$lonRad = $x * M_PI / 180;
		$recf = 1 / $f;
		$b = $a * ($recf - 1) / $recf;
		$eSquared = self::CalculateESquared($a, $b);
		$e2Squared = self::CalculateE2Squared ($a, $b);
		$tn = ($a - $b) / ($a + $b);
		$ap = $a * (1 - $tn + 5 * (pow($tn, 2) - pow($tn, 3)) / 4 + 81 * (pow($tn, 4) - pow($tn, 5)) / 64);
		$bp = 3 * $a * ($tn - pow($tn, 2) + 7 * (pow($tn, 3) - pow($tn, 4)) / 8 + 55 * pow($tn, 5) / 64) / 2;
		$cp = 15 * $a * (pow($tn, 2) - pow($tn, 3) + 3 * (pow($tn, 4) - pow($tn, 5)) / 4) / 16;
		$dp = 35 * $a * (pow($tn, 3) - pow($tn, 4) + 11 * pow($tn, 5) / 16) / 48;
		$ep = 315 * $a * (pow($tn, 4) - pow($tn, 5)) / 512;
		$dlam = $lonRad - $olam;
		$s = sin($latRad);
		$c = cos($latRad);
		$t = $s / $c;
		$eta = $e2Squared * pow($c, 2);
		$sn = self::sphsn($a, $eSquared, $latRad);
		$tmd = self::sphtmd($ap, $bp, $cp, $dp, $ep, $latRad);
		$t1 = $tmd * $ok;
		$t2 = $sn * $s * $c * $ok / 2;
		$t3 = $sn * $s * pow($c, 3) * $ok * (5 - pow($t, 2) + 9 * $eta + 4 * pow($eta, 2)) / 24;
		$t4 = $sn * $s * pow($c, 5) * $ok * (61 - 58 * pow($t, 2) + pow($t, 4) + 270 * $eta - 330 * pow($t, 2) * $eta + 445 * pow($eta, 2) + 324 * pow($eta, 3) - 680 * pow($t, 2) * pow($eta, 2) + 88 * pow($eta, 4) - 600 * pow($t, 2) * pow($eta, 3) - 192 * pow($t, 2) * pow($eta, 4)) / 720;
		$t5 = $sn * $s * pow($c, 7) * $ok * (1385 - 3111 * pow($t, 2) + 543 * pow($t, 4) - pow($t, 6)) / 40320;
		$t6 = $sn * $c * $ok;
		$t7 = $sn * pow($c, 3) * $ok * (1 - pow($t, 2) + $eta) / 6;
		$t8 = $sn * pow($c, 5) * $ok * (5 - 18 * pow($t, 2) + pow($t, 4) + 14 * $eta - 58 * pow($t, 2) * $eta + 13 * pow($eta, 2) + 4 * pow($eta, 3) - 64 * pow($t, 2) * pow($eta, 2) - 24 * pow($t, 2) * pow($eta, 3)) / 120;
		$t9 = $sn * pow($c, 7) * $ok * (61 - 479 * pow($t, 2) + 179 * pow($t, 4) - pow($t, 6)) / 5040;
		$result->x = self::$fe + $strf + $dlam * $t6 + pow($dlam, 3) * $t7  + pow($dlam, 5) * $t8 + pow($dlam, 7) * $t9;
		$result->y = $nfn + $t1 + pow($dlam, 2) * $t2 + pow($dlam, 4) * $t3 + pow($dlam, 6) * $t4 + pow($dlam, 8) * $t5;
		return $result;
	}

	private static function PuwgToLonLat ($a, $f, $x, $y, $proj) {
		$result = new StdClass();
		$ok =  0.999923;
		if ($proj == 1) $ok = 0.9993;
		$recf = 1 / $f;
		$b = $a * ($recf - 1) / $recf;
		$eSquared = self::CalculateESquared($a, $b);
		$e2Squared = self::CalculateE2Squared($a, $b);
		$tn = ($a - $b) / ($a + $b);
		$ap = $a * (1 - $tn + 5 * (pow($tn, 2) - pow($tn, 3)) / 4 + 81 * (pow($tn, 4) - pow($tn, 5)) / 64);
		$bp = 3 * $a * ($tn - pow($tn, 2) + 7 * (pow($tn, 3) - pow($tn, 4)) / 8 + 55 * pow($tn, 5) / 64) / 2;
		$cp = 15 * $a * (pow($tn, 2) - pow($tn, 3) + 3 * (pow($tn, 4) - pow($tn, 5)) / 4) / 16;
		$dp = 35 * $a * (pow($tn, 3) - pow($tn, 4) + 11 * pow($tn, 5) / 16) / 48;
		$ep = 315 * $a * (pow($tn, 4) - pow($tn, 5)) / 512;
		if ($proj == 1) {
			$olam = 19 * M_PI / 180;
			$strf = 0;
			$nfn = -5300000;
          
		} else {
			$nfn = 0;
			$strf = floor($x / 1000000) * 1000000;
			$olam = $strf / 60000000 * M_PI;
		}
		$tmd = ($y - $nfn) / $ok;
		$sr = self::sphsr($a, $eSquared, 0);
		$ftphi = $tmd / $sr;
		for ($i = 0; $i < 5; $i++) {
			$t10 = self::sphtmd($ap, $bp, $cp, $dp, $ep, $ftphi);
			$sr = self::sphsr($a, $eSquared, $ftphi);
			$ftphi = $ftphi + ($tmd - $t10) / $sr;
		}
		$sr = self::sphsr($a, $eSquared, $ftphi);
		$sn = self::sphsn($a, $eSquared, $ftphi);
		$s = sin($ftphi);
		$c = cos($ftphi);
		$t = $s / $c;
		$eta = $e2Squared * pow($c, 2);
		$de = $x - self::$fe - $strf;
		$t10 = $t / (2 * $sr * $sn * pow($ok, 2));
		$t11 = $t * (5 + 3 * pow($t, 2) + $eta - 4 * pow($eta, 2) - 9 * pow($t, 2) * $eta) / (24 * $sr * pow($sn, 3) * pow($ok, 4));
        
		$t12 = $t * (61 + 90 * pow($t, 2) + 46 * $eta + 45 * pow($t, 4) - 252 * pow($t, 2) * $eta - 3 * pow($eta, 2) + 100 * pow($eta, 3) - 66 * pow($t, 2) * pow($eta, 2) - 90 * pow($t, 4) * $eta + 88 * pow($eta, 4) + 225 * pow($t, 4) * pow($eta, 2) + 84 * pow($t, 2) * pow($eta, 3) - 192 * pow($t, 2) * pow($eta, 4)) / (720 * $sr * pow($sn, 5) * pow($ok, 6));
		$t13 = $t * (1385 + 3633 * pow($t, 2) + 4095 * pow($t, 4) + 1575  * pow($t, 6)) / (40320 * $sr * pow($sn, 7) * pow($ok, 8));
		$t14 = 1 / ($sn * $c * $ok);
		$t15 = (1 + 2 * pow($t, 2) + $eta) / (6 * pow($sn, 3) * $c * pow($ok, 3));
		$t16 = 1 * (5 + 6 * $eta + 28 * pow($t, 2) - 3 * pow($eta, 2) + 8 * pow($t, 2) * $eta + 24 * pow($t, 4) - 4 * pow($eta, 3) + 4 * pow($t, 2) * pow($eta, 2) + 24 * pow($t, 2) * pow($eta, 3)) / (120 * pow($sn, 5) * $c * pow($ok, 5));
		$t17 = 1 * (61 + 662 * pow($t, 2) + 1320 * pow($t, 4) + 720 * pow($t, 6)) / (5040 * pow($sn, 7) * $c * pow($ok, 7));
		$dlam = $de * $t14 - pow($de, 3) * $t15 + pow($de, 5) * $t16 - pow($de, 7) * $t17;
		$result->x = ($olam + $dlam) * 180 / M_PI;
		$result->y = ($ftphi - pow($de, 2) * $t10 + pow($de, 4) * $t11 - pow($de, 6) * $t12 + pow($de, 8) * $t13) * 180 / M_PI;
		return $result;
	}

	public static function getEpsgPuwg2000($x, $y = 0) {
		if ($x < 13.5 || $x >= 25.5) throw new Exception("Wrong longitude - {$x}");
		return floor(($x + 1.5) / 3) + 2171;
	}

	public static function Wgs84ToPuwg2000($x, $y) {
		$result = self::LonLatToPuwg(6378137, 1 / 298.257223563, $x, $y, 4);
		$result->epsg = self::getEpsgPuwg2000($x);
		return $result;
	}

	public static function Puwg2000ToWgs84($x, $y) {
		$result = self::PuwgToLonLat(6378137, 1 / 298.257223563, $x, $y, 4);
		$result->epsg = self::getEpsgPuwg2000($result->x);
		return $result;
	}

	public static function Wgs84ToEpsg2180($x, $y) {
		// Parametry elipsoidy GRS-80i
		$e = 0.0818191910428;		//pierwszy mimośród elipsoidy
		$R0 = 6367449.14577;		//promieñ sfery Lagrange'a
		$Snorm = 0.000002;		//parametr normujący
		$xo = 5760000;			//parametr centrujący

		//Współczynniki wielomianu
		$a0 = 5765181.11148097;
		$a1 = 499800.81713800;
		$a2 = -63.81145283;
		$a3 = 0.83537915;
		$a4 = 0.13046891;
		$a5 = -0.00111138;
		$a6 = -0.00010504;

		// Parametry odwzorowania Gaussa-Kruegera dla układu PUWG92
		$L0_stopnie = 19;		//Początek układu wsp. PUWG92 (długość)
		$m0 = 0.9993;
		$x0 = -5300000;
		$y0 = 500000;

		// Zakres stosowalności metody
		$Bmin = 48 * M_PI / 180;
		$Bmax = 56 * M_PI / 180;
		$dLmin = -6 * M_PI / 180;
		$dLmax = 6 * M_PI / 180;

		// Weryfikacja danych wejściowych
		$B = $y * M_PI / 180;
		$dL_stopnie = $x - $L0_stopnie;
		$dL = $dL_stopnie * M_PI / 180;

		if ($B < $Bmin || $B > $Bmax) throw new Exception("Wrong lattitude - {$y}");
		if ($dL < $dLmin || $dL > $dLmax) throw new Exception("Wrong longitude - {$x}");

		//etap I - elipsoida na kulę
		$U = 1 - $e * sin($B);
		$V = 1 + $e * sin($B);
		$K = pow(($U /$V), ($e / 2));
		$C = $K * tan($B / 2 + M_PI / 4);
		$fi = 2 * atan($C) - M_PI / 2;
		$d_lambda = $dL;

		// etap II - kula na walec
		$p = sin($fi);
		$q = cos($fi) * cos($d_lambda);
		$r = 1 + cos($fi) * sin($d_lambda);
		$s = 1 - cos($fi) * sin($d_lambda);
		$XMERC = $R0 * atan($p / $q);
		$YMERC = 0.5 * $R0 * log($r / $s);

		//etap III - walec na płaszczyznę
		$Z = new Complex(($XMERC - $xo) * $Snorm, $YMERC * $Snorm);
		$Zgk = (new Complex())->add($a6)->multiply($Z)->add($a5)->multiply($Z)->add($a4)->multiply($Z)->add($a3)->multiply($Z)->add($a2)->multiply($Z)->add($a1)->multiply($Z)->add($a0);
		$Xgk = $Zgk->real();
		$Ygk = $Zgk->imag();

		//Przejście do układu aplikacyjnego
		$result = new StdClass();
		$result->x = $m0 * $Xgk + $x0;
		$result->y = $m0 * $Ygk + $y0;
		$result->epsg = 2180;

		return $result;
	}

	public static function Wgs84ToPuwg1992($x, $y) {
		return self::Wgs84ToEpsg2180($x, $y);
	}

	public static function Epsg2180ToWgs84($x, $y) {
		$L0_stopnie = 19;		//Początek układu wsp. PUWG92 (długość)
		$m0 = 0.9993;
		$x0 = -5300000;
		$y0 = 500000;

		$R0 = 6367449.14577;		//promieñ sfery Lagrange'a
		$Snorm = 0.000002;		//parametr normujący
		$xo_prim=5765181.11148097;	//parametr centrujący

		// Współczynniki wielomianu
		$b0 = 5760000;
		$b1 = 500199.26224125;
		$b2 = 63.88777449;
		$b3 = -0.82039170;
		$b4 = -0.13125817;
		$b5 = 0.00101782;
		$b6 = 0.00010778;

		// Współczynniki szeregu tryg.
		$c2 = 0.0033565514856;
		$c4 = 0.0000065718731;
		$c6 = 0.0000000176466;
		$c8 = 0.0000000000540;

		//Przejście z układu aplikacyjnego
		$Xgk = ($x - $x0) / $m0;
		$Ygk = ($y - $y0) / $m0;

		//etap I - (Xgk, Ygk) -> (Xmerc, Ymerc)
		$Z = new Complex(($Xgk - $xo_prim) * $Snorm, $Ygk * $Snorm);
		$Zmerc = (new Complex())->add($b6)->multiply($Z)->add($b5)->multiply($Z)->add($b4)->multiply($Z)->add($b3)->multiply($Z)->add($b2)->multiply($Z)->add($b1)->multiply($Z)->add($b0);
		$Xmerc = $Zmerc->real();
		$Ymerc = $Zmerc->imag();

		//etap II - Xmerc,Ymerc -> fi, delta_lambda
		$alfa = $Xmerc / $R0;
		$beta = $Ymerc / $R0;
		$w = 2 * atan(exp($beta)) - M_PI / 2;
		$fi = asin(cos($w) * sin($alfa));
		$d_lambda = atan(tan($w) / cos($alfa));

		//etap III
		$B = $fi + $c2 * sin(2 * $fi) + $c4 * sin(4 *$fi) + $c6 * sin(6 * $fi) + $c8 * sin(8 * $fi);
		$dL = $d_lambda;

		//Obliczenia koncowe
		$result = new StdClass();
		$result->x = $dL / M_PI * 180 + $L0_stopnie;
		$result->y = $B / M_PI * 180;
		$result->epsg = 2180;

		return $result;
	}

	public static function Puwg1992ToWgs84($x, $y) {
		return self::Epsg2180ToWgs84($x, $y);
	}

	public static function GetZByEpsg2180($x, $y) {
		if (!(is_numeric($x) && is_numeric($y))) throw new Exception('Bad argument(s)');
		$table = 'CODGIK_NMT_100_EPSG2180';
		$query = "select x(`the_geom`) as `x`, y(`the_geom`) as `y`, `z` from `{$table}` where abs(x(`the_geom`) - {$x}) < 100 and abs(y(`the_geom`) - {$y}) < 100";
		try {
			$result = DB::getPDO()->fetchAll($query);
		} catch (Exception $e) {
			throw new Exception('Error while connecting to database');
		}
		$count = count($result);
		if (!$count) throw new Exception('Out of the area');
		$z = [];
		foreach ($result as $row) {
			$z[$row['x']][$row['y']] = $row['z'];
			$_x[$row['x']] = true;
			$_y[$row['y']] = true;
		}
		$xmin = min(array_keys($_x));
		$xmax = max(array_keys($_x));
		$ymin = min(array_keys($_y));
		$ymax = max(array_keys($_y));
		if ($count == 1) {
			return $z[$xmin][$ymin];
		} elseif ($count == 2) {
			if ($xmin == $xmax && $ymax - $ymin == 100) {
				return (($y - $ymin) * $z[$xmin][$ymax] + ($ymax - $y) * $z[$xmin][$ymin]) / 100;
			} elseif ($ymin == $ymax && $xmax - $xmin == 100) {
				return (($x - $xmin) * $z[$xmax][$ymin] + ($xmax - $x) * $z[$xmin][$ymin]) / 100;
			} else throw new Exception('Unknown error #2');
		} elseif ($count == 4) {
			if (!($ymax - $ymin == 100 && $xmax - $xmin == 100)) throw new Exception('Unknown error #3');
			return (($y - $ymin) * (($x - $xmin) * $z[$xmax][$ymax] + ($xmax - $x) * $z[$xmin][$ymax])
				+ ($ymax - $y) * (($x - $xmin) * $z[$xmax][$ymin] + ($xmax - $x) * $z[$xmin][$ymin])) / 10000;
		} else throw new Exception('Unknown error #1');
	}

	public static function GetZByWgs84($x, $y) {
		$epsg2180 = self::Wgs84ToEpsg2180($x, $y);
		return self::GetZByEpsg2180($epsg2180->x, $epsg2180->y);
	}
}