Merge branch 'master' of bn.git:plabudda/se

SE/se-lib/EpsgConversion.php

@@ -1,5 +1,49 @@
 <?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;
@@ -31,8 +75,8 @@ class EpsgConversion {
 		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, $lon, $lat, $proj) {
-		if ($lon < 13.5 || $lon >= 25.5) throw new Exception("Wrong longitude - {$lon}");
+	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;
@@ -41,14 +85,14 @@ class EpsgConversion {
 			$ok = 0.9993;
 		} else {
 			$nfn = 0;
-			$olam = floor(($lon + 1.5)/3) * M_PI / 60;
-			$strf = floor(($lon + 1.5)/3) * 1000000;
+			$olam = floor(($x + 1.5)/3) * M_PI / 60;
+			$strf = floor(($x + 1.5)/3) * 1000000;
 			$ok = 0.999923;
 		}
 		
 		$result = new StdClass();
-		$latRad = $lat * M_PI / 180;
-		$lonRad = $lon * M_PI / 180;
+		$latRad = $y * M_PI / 180;
+		$lonRad = $x * M_PI / 180;
 		$recf = 1 / $f;
 		$b = $a * ($recf - 1) / $recf;
 		$eSquared = self::CalculateESquared($a, $b);
@@ -80,7 +124,7 @@ class EpsgConversion {
 		return $result;
 	}
 
-	private static function PUWGToLonLat ($a, $f, $lon, $lat, $proj) {
+	private static function PuwgToLonLat ($a, $f, $x, $y, $proj) {
 		$result = new StdClass();
 		$ok =  0.999923;
 		if ($proj == 1) $ok = 0.9993;
@@ -101,10 +145,10 @@ class EpsgConversion {
           
 		} else {
 			$nfn = 0;
-			$strf = floor($lon / 1000000) * 1000000;
+			$strf = floor($x / 1000000) * 1000000;
 			$olam = $strf / 60000000 * M_PI;
 		}
-		$tmd = ($lat - $nfn) / $ok;
+		$tmd = ($y - $nfn) / $ok;
 		$sr = self::sphsr($a, $eSquared, 0);
 		$ftphi = $tmd / $sr;
 		for ($i = 0; $i < 5; $i++) {
@@ -118,7 +162,7 @@ class EpsgConversion {
 		$c = cos($ftphi);
 		$t = $s / $c;
 		$eta = $e2Squared * pow($c, 2);
-		$de = $lon - self::$fe - $strf;
+		$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));
         
@@ -134,21 +178,191 @@ class EpsgConversion {
 		return $result;
 	}
 
-	public static function getEpsg($lon, $lat = 0) {
-		if ($lon < 13.5 || $lon >= 25.5) throw new Exception("Wrong longitude - {$lon}");
-		return floor(($lon + 1.5) / 3) + 2171;
+	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 LonLatToPUWGWGS84($lon, $lat) {
-		$result = self::LonLatToPUWG(6378137, 1 / 298.257223563, $lon, $lat, 4);
-		$result->epsg = self::getEpsg($lon);
+	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 PUWGToLonLatWGS84($lon, $lat) {
-		$result = self::PUWGToLonLat(6378137, 1 / 298.257223563, $lon, $lat, 4);
-		$result->epsg = self::getEpsg($result->x);
+	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);
+	}
 }

SE/se-lib/Route/GeoreferencesManager.php

@@ -64,7 +64,7 @@ class Route_GeoreferencesManager extends RouteBase {
 			if ($point['x'] > 0 && $point['y'] > 0) {
 				if ($point['gx'] && $point['gy']) {	
 					try {
-						$test = epsgConversion::LonLatToPUWGWGS84($point['gx'], $point['gy']);
+						$test = epsgConversion::Wgs84ToPuwg2000($point['gx'], $point['gy']);
 						if (abs($point['x'] - $test->x) > 0.001) $errors["unfuckable"][$point["ID"]]["errors"][] = "x";
 						if (abs($point['y'] - $test->y) > 0.001) $errors["unfuckable"][$point["ID"]]["errors"][] = "y";
 						if ($point['EPSG'] != $test->epsg) $errors["unfuckable"][$point["ID"]]["errors"][] = "epsg";
@@ -84,7 +84,7 @@ class Route_GeoreferencesManager extends RouteBase {
 		if (isset($errorsData["unfuckable"])) {
 			$echo = '<b>Znalezione błędy, które można automatycznie naprawić:</b></br/><pre>';
 			foreach ($errorsData["unfuckable"] as $error) {
-				$ok = epsgConversion::PUWGToLonLatWGS84($error['data']['x'], $error['data']['y']);
+				$ok = epsgConversion::Puwg2000ToWgs84($error['data']['x'], $error['data']['y']);
 				$echo .= "[{$error['data']['ID']}] x = <span style='color:" . ((in_array("x", $error["errors"])) ? "red" : "green") .
 					";'>{$error['data']['gx']}</span> ({$ok->x}), y = <span style='color:" . ((in_array("y", $error["errors"])) ? "red" : "green") .
 					";'>{$error['data']['gy']}</span> ({$ok->y}), epsg = <span style='color:" . ((in_array("epsg", $error["errors"])) ? "red" : "green") .
@@ -114,7 +114,7 @@ class Route_GeoreferencesManager extends RouteBase {
 		$updateErrors = [];
 		$updateSuccess = 0;
 		foreach ($errorsData as $ID => $error) {
-			$ok = epsgConversion::PUWGToLonLatWGS84($error['data']['x'], $error['data']['y']);
+			$ok = epsgConversion::Puwg2000ToWgs84($error['data']['x'], $error['data']['y']);
 			$sqlArr = [
 				"ID" => $ID,
 				"the_geom" => "GeomFromText('POINT({$ok->x} {$ok->y})')",
@@ -175,7 +175,7 @@ class Route_GeoreferencesManager extends RouteBase {
 			DB::getPDO()->query("CREATE TEMPORARY TABLE `{$tempTbl}` (SELECT * FROM `" . self::TABLE . "`)");
 			$result = DB::getPDO()->fetchall("SELECT ID, x(the_geom) AS gx, y(the_geom) AS gy FROM `{$tempTbl}` WHERE ST_IsEmpty(the_geom) = 0 AND x = 0 AND y = 0");
 			foreach ($result as $row) {
-				$puwg = epsgConversion::LonLatToPUWGWGS84($row['gx'], $row['gy']);
+				$puwg = epsgConversion::Wgs84ToPuwg2000($row['gx'], $row['gy']);
 				DB::getPDO()->query("UPDATE `{$tempTbl}` SET x = '{$puwg->x}', y = '{$puwg->y}' WHERE ID='{$row['ID']}'");
 			}
 
@@ -189,7 +189,7 @@ class Route_GeoreferencesManager extends RouteBase {
 					$closePoints[$lp] = $result[0]['ID'];
 				}
 				else $closePointsDetail[$lp] = false;
-				$wgs84[$lp] = epsgConversion::PUWGToLonLatWGS84($point['x'], $point['y']);
+				$wgs84[$lp] = epsgConversion::Puwg2000ToWgs84($point['x'], $point['y']);
 			}
 ?>
     <style>
@@ -269,7 +269,7 @@ class Route_GeoreferencesManager extends RouteBase {
 				$y = $data['points'][$lp]['y'];
 				$z = $data['points'][$lp]['z'];
 				try {
-					$wgs84 = epsgConversion::PUWGToLonLatWGS84($x, $y);
+					$wgs84 = epsgConversion::Puwg2000ToWgs84($x, $y);
 					$gx = $wgs84->x;
 					$gy = $wgs84->y;
 					$epsg = $wgs84->epsg;

SE/se-lib/Route/ViewTableAjax.php

@@ -654,8 +654,13 @@ class Route_ViewTableAjax extends RouteBase {
 		$points = explode(',', $matches[1]);
 		$csv = implode("\n", array_map(function ($point, $i) {
 			list($x, $y) = explode(" ", $point, 2);
-			$wgs84 = EpsgConversion::LonLatToPUWGWGS84($x, $y);
-			return $i++ . ',' . round($wgs84->y, 3) . ',' . round($wgs84->x, 3) . ',0,Pikieta';
+			$puwg2000 = EpsgConversion::Wgs84ToPuwg2000($x, $y);
+			try {
+				$z = EpsgConversion::GetZByWgs84($x, $y);
+			} catch (Exception $e) {
+				$z = 0;
+			}
+			return $i++ . ',' . round($puwg2000->y, 3) . ',' . round($puwg2000->x, 3) . ',' . round($z, 3) . ',Punkt';
 		}, $points, range(1, count($points))));
 		Response::sendCsv($csv, "{$table}.{$id}");
 	}