Zarządzanie punktami georeferencyjnymi (w tym klasa do przeliczania WGS84->EPSG:2177 i na odwrót)

SE/se-lib/EpsgConversion.php

@@ -0,0 +1,155 @@
+<?php
+
+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, $lon, $lat, $proj) {
+		if ($lon < 13.5 || $lon >= 25.5) throw new Exception("Wrong longitude - {$lon}");
+
+		if ($proj == 1) {
+			$alam = 19 * M_PI / 180;
+			$strf = 0;
+			$nfn = -5300000.0;
+			$ok = 0.9993;
+		} else {
+			$nfn = 0;
+			$olam = floor(($lon + 1.5)/3) * M_PI / 60;
+			$strf = floor(($lon + 1.5)/3) * 1000000;
+			$ok = 0.999923;
+		}
+		
+		$result = new StdClass();
+		$latRad = $lat * M_PI / 180;
+		$lonRad = $lon * 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, $lon, $lat, $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($lon / 1000000) * 1000000;
+			$olam = $strf / 1000000 * M_PI / 180;
+		}
+		$tmd = ($lat - $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 = $lon - 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));
+		$result->y = $ftphi - pow($de, 2) * $t10 + pow($de, 4) * $t11 - pow($de, 6) * $t12 + pow($de, 8) * $t13;
+		$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;
+		$this->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 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 LonLatToPUWGWGS84($lon, $lat) {
+		$result = self::LonLatToPUWG(6378137, 1 / 298.257223563, $lon, $lat, 4);
+		$result->epsg = self::getEpsg($lon);
+		return $result;
+	}
+
+	public static function PUWGToLonLatWGS84($lon, $lat) {
+		$result = self::PUWGToLonLat(6378137, 1 / 298.257223563, $lon, $lat, 4);
+		$result->epsg = self::getEpsg($result->x);
+		return $result;
+	}
+
+}

SE/se-lib/Route/GeoreferencesManager.php

@@ -0,0 +1,89 @@
+<?php
+
+Lib::loadClass('RouteBase');
+Lib::loadClass('EpsgConversion');
+
+class Route_GeoreferencesManager extends RouteBase {
+
+	protected $ACTION, $REFERER;
+	const TABLE = "WMS_MAP_GEOREFERENCES";
+
+	public function handleAuth() {
+		if (!User::logged()) {
+			User::authByRequest();
+		}
+		$this->REFERER = $_SERVER['HTTP_REFERER'];
+		$this->ACTION = V::get('action', '', $_GET, '');
+	}
+
+	public function defaultAction() {
+		SE_Layout::gora();
+		SE_Layout::menu();
+?>
+<div class="container" style="margin-top:20px">
+  <legend>Zarządzanie punktami georeferencyjnymi</legend>
+  <div class="form-group">
+    <div class="col-sm-12">
+      <a href="?_route=GeoreferencesManager&action=verifyPoints" class="btn-sm btn-primary">Zweryfikuj punkty</a>
+    </div>
+  </div>
+<?php
+		switch ($this->ACTION) {
+			case "verifyPoints":
+				$this->verifyPoints();
+				break;
+		}
+?>
+</div>
+<?php
+		SE_Layout::dol();
+	}
+
+	public function reinstallAction() {
+		$this->reinstall();
+		die('OK');
+	}
+
+	public function reinstall() {
+	}
+
+	private function verifyPoints() {
+		$points = DB::getPDO()->fetchall("select `ID`, `EPSG`, x(`the_geom`) as `gx`, y(`the_geom`) as `gy`, `x`, `y` from `" . self::TABLE . "`");
+		$errors = [];
+		foreach ($points as $point) {
+			if ($point['x'] && $point['y']) {
+				if ($point['gx'] && $point['gy']) {
+					try {
+						$test = epsgConversion::LatLonToPUWGWGS84($point['gx'], $point['gy']);
+					} catch (Exception $e) {
+						SE_Layout::alert('danger', $e->getMessage());
+					}
+					if ((abs($point['x'] - $test->x) > 1 ) || (abs($point['y'] - $test->y) > 1) || ($point['EPSG'] != $test->epsg)) {
+						$error = "[{$point['ID']}] x = <span style='color:" . ((abs($point['x'] - $test->x) > 1) ? "red" : "green") .
+							";'>{$point['x']}</span> ({$test->x}), y = <span style='color:" . ((abs($point['y'] - $test->y) > 1) ? "red" : "green") .
+							";'>{$point['y']}</span> ({$test->y}), epsg = <span style='color:" . (($point['EPSG'] != $test->epsg) ? "red" : "green") .
+							";'>{$point['EPSG']}</span>" . (($point['EPSG'] != $test->epsg) ? " ({$test->epsg})" : "");
+						$errors[] = $error;
+					}
+				}
+			}
+		}
+?>
+  <div class="form-group">
+    <div class="col-sm-12"><br/>
+<?php
+//$errors = [];
+//		if ($errors) echo "<pre>" . implode("<br/>", $errors) . "</pre>";
+		if ($errors) {
+			SE_Layout::alert('danger', "<pre>" . implode("<br/>", $errors) . "</pre>");
+?>
+    <div style="text-align:center;"><a href="?_route=GeoreferencesManager&action=repairPoints" class="btn btn-primary">Napraw wszystkie</a></div>
+<?php
+		} else SE_Layout::alert('success', "Wszystko OK"); //echo "<span style='color:green;'>Wszystko OK</span>";
+?>
+    </div>
+  </div>
+<?php
+	}
+
+}

SE/se-lib/TableAjax.php

@@ -436,6 +436,9 @@ class TableAjax extends ViewAjax {
 				<?php elseif ($this->allowTreeView()) : ?>
 					<a class="pull-right" style="padding:0 20px 0 0;" href="<?php echo "index.php?MENU_INIT=VIEWTREE_AJAX&ZASOB_ID={$this->_zasobID}"; ?>"><i class="glyphicon glyphicon-eye-open"></i> Drzewo</a>
 				<?php endif; ?>
+				<?php if ($this->_tbl == 'WMS_MAP_GEOREFERENCES') : ?>
+					<a class="pull-right" style="padding:0 20px 0 0;" href="<?php echo "index.php?_route=GeoreferencesManager"; ?>"><i class="glyphicon glyphicon glyphicon-wrench"></i> Zarządzaj punktami</a>
+				<?php endif; ?>
 			</ul>
 			<div id="<?php echo $this->_htmlID; ?>"></div>
 		</div>