|
1
|
|
|
<?php |
|
2
|
|
|
|
|
3
|
|
|
namespace kalanis\kw_coordinates\Codecs; |
|
4
|
|
|
|
|
5
|
|
|
|
|
6
|
|
|
use kalanis\kw_coordinates\Interfaces; |
|
7
|
|
|
use kalanis\kw_coordinates\Support; |
|
8
|
|
|
|
|
9
|
|
|
|
|
10
|
|
|
/** |
|
11
|
|
|
* Class Sjtsk1 |
|
12
|
|
|
* @package kalanis\kw_coordinates\Codecs |
|
13
|
|
|
* @link https://blog.drndos.sk/2013/03/how-to-convert-s-jtsk-coordinates-to-wsg84/ |
|
14
|
|
|
* @link https://www.pecina.cz/krovak.html |
|
15
|
|
|
* |
|
16
|
|
|
* SJTSKtoWSG84.java Filip Bednárik <[email protected]> |
|
17
|
|
|
* package sk.datalan.semweb.czechpediacrawler; |
|
18
|
|
|
* |
|
19
|
|
|
* toWgs is port from java class |
|
20
|
|
|
* fromWgs is port from JavaScript code |
|
21
|
|
|
* |
|
22
|
|
|
* Do not ask me how both works...; but used, tested and works |
|
23
|
|
|
*/ |
|
24
|
|
|
class Sjtsk1 implements Interfaces\ICodecs |
|
25
|
|
|
{ |
|
26
|
|
|
const FG_SHIFT = 17.66666666666666; |
|
27
|
|
|
const dX_SHIFT = 589; |
|
28
|
|
|
const dY_SHIFT = 76; |
|
29
|
|
|
const dZ_SHIFT = 480; |
|
30
|
|
|
|
|
31
|
|
|
protected Interfaces\INumbers $pos; |
|
32
|
|
|
protected Interfaces\IFormatted $out; |
|
33
|
|
|
|
|
34
|
4 |
|
public function __construct(?Interfaces\INumbers $position = null, ?Interfaces\IFormatted $output = null) |
|
35
|
|
|
{ |
|
36
|
4 |
|
$this->pos = $position ?: new Support\Position(); |
|
37
|
4 |
|
$this->out = $output ?: new Support\JtskObject(); |
|
38
|
4 |
|
} |
|
39
|
|
|
|
|
40
|
2 |
|
public function fromLonLat(Interfaces\INumbers $position, array $params): Interfaces\IFormatted |
|
41
|
|
|
{ |
|
42
|
2 |
|
$d2r = pi() / 180; |
|
43
|
2 |
|
$a = 6378137.0; |
|
44
|
2 |
|
$f1 = 298.257223563; |
|
45
|
2 |
|
$dx = -570.69; |
|
46
|
2 |
|
$dy = -85.69; |
|
47
|
2 |
|
$dz = -462.84; |
|
48
|
2 |
|
$wx = 4.99821 / 3600 * pi() / 180; |
|
49
|
2 |
|
$wy = 1.58676 / 3600 * pi() / 180; |
|
50
|
2 |
|
$wz = 5.2611 / 3600 * pi() / 180; |
|
51
|
2 |
|
$m = -3.543e-6; |
|
52
|
|
|
|
|
53
|
2 |
|
$B = $position->getLatitude() * $d2r; |
|
54
|
2 |
|
$L = $position->getLongitude() * $d2r; |
|
55
|
2 |
|
$H = $position->getAltitude(); |
|
56
|
|
|
|
|
57
|
2 |
|
$e2 = 1 - $this->sqr(1 - (1 / $f1)); |
|
58
|
2 |
|
$rho = $a / sqrt(1 - ($e2 * $this->sqr(sin($B)))); |
|
59
|
2 |
|
$x1 = ($rho + $H) * cos($B) * cos($L); |
|
60
|
2 |
|
$y1 = ($rho + $H) * cos($B) * sin($L); |
|
61
|
2 |
|
$z1 = (((1 - $e2) * $rho) + $H) * sin($B); |
|
62
|
|
|
|
|
63
|
2 |
|
$x2 = $dx + (1 + $m) * ($x1 + ($wz * $y1) - ($wy * $z1)); |
|
64
|
2 |
|
$y2 = $dy + (1 + $m) * (((-$wz) * $x1) + $y1 + ($wx * $z1)); |
|
65
|
2 |
|
$z2 = $dz + (1 + $m) * (($wy * $x1) - ($wx * $y1) + $z1); |
|
66
|
|
|
|
|
67
|
2 |
|
$a = 6377397.15508; |
|
68
|
2 |
|
$f1 = 299.152812853; |
|
69
|
2 |
|
$ab = $f1 / ($f1 - 1); |
|
70
|
2 |
|
$p = sqrt($this->sqr($x2) + $this->sqr($y2)); |
|
71
|
2 |
|
$e2 = 1 - $this->sqr(1 - (1 / $f1)); |
|
72
|
2 |
|
$th = atan(($z2 * $ab) / $p); |
|
73
|
2 |
|
$st = sin($th); |
|
74
|
2 |
|
$ct = cos($th); |
|
75
|
2 |
|
$t = ($z2 + ($e2 * $ab * $a * ($st * $st * $st))) / ($p - ($e2 * $a * ($ct * $ct * $ct))); |
|
76
|
|
|
|
|
77
|
2 |
|
$B = atan($t); |
|
78
|
2 |
|
$H = sqrt(1 + ($t * $t)) * ($p - ($a / sqrt(1 + ((1 - $e2) * $t * $t)))); |
|
79
|
2 |
|
$L = 2 * atan($y2 / ($p + $x2)); |
|
80
|
|
|
|
|
81
|
2 |
|
$a = 6377397.15508; |
|
|
|
|
|
|
82
|
2 |
|
$e = 0.081696831215303; |
|
83
|
2 |
|
$n = 0.97992470462083; |
|
84
|
2 |
|
$rho0 = 12310230.12797036; |
|
85
|
2 |
|
$sinUQ = 0.863499969506341; |
|
86
|
2 |
|
$cosUQ = 0.504348889819882; |
|
87
|
2 |
|
$sinVQ = 0.420215144586493; |
|
88
|
2 |
|
$cosVQ = 0.907424504992097; |
|
89
|
2 |
|
$alpha = 1.000597498371542; |
|
90
|
2 |
|
$k2 = 1.00685001861538; |
|
91
|
|
|
|
|
92
|
2 |
|
$sinB = sin($B); |
|
93
|
2 |
|
$t = (1 - ($e * $sinB)) / (1 + ($e * $sinB)); |
|
94
|
2 |
|
$t = ($this->sqr(1 + $sinB) / (1 - $this->sqr($sinB))) * exp($e * log($t)); |
|
95
|
2 |
|
$t = $k2 * exp($alpha * log($t)); |
|
96
|
|
|
|
|
97
|
2 |
|
$sinU = ($t - 1) / ($t + 1); |
|
98
|
2 |
|
$cosU = sqrt(1 - ($sinU * $sinU)); |
|
99
|
2 |
|
$V = $alpha * $L; |
|
100
|
2 |
|
$sinV = sin($V); |
|
101
|
2 |
|
$cosV = cos($V); |
|
102
|
2 |
|
$cosDV = ($cosVQ * $cosV) + ($sinVQ * $sinV); |
|
103
|
2 |
|
$sinDV = ($sinVQ * $cosV) - ($cosVQ * $sinV); |
|
104
|
2 |
|
$sinS = ($sinUQ * $sinU) + ($cosUQ * $cosU * $cosDV); |
|
105
|
2 |
|
$cosS = sqrt(1 - $sinS * $sinS); |
|
106
|
2 |
|
$sinD = ($sinDV * $cosU) / $cosS; |
|
107
|
2 |
|
$cosD = sqrt(1 - ($sinD * $sinD)); |
|
108
|
|
|
|
|
109
|
2 |
|
$eps = $n * atan($sinD / $cosD); |
|
110
|
2 |
|
$rho = $rho0 * exp(-$n * log((1 + $sinS) / $cosS)); |
|
111
|
|
|
|
|
112
|
2 |
|
$CX = $rho * sin($eps); |
|
113
|
2 |
|
$CY = $rho * cos($eps); |
|
114
|
|
|
|
|
115
|
2 |
|
$obj = clone $this->out; |
|
116
|
|
|
|
|
117
|
2 |
|
$mul = isset($params['czech']) && boolval($params['czech']) ? 1 : -1; |
|
118
|
|
|
|
|
119
|
2 |
|
return $obj->setData( |
|
120
|
2 |
|
$mul * $CX, |
|
121
|
2 |
|
$mul * $CY, |
|
122
|
|
|
$H |
|
123
|
|
|
); |
|
124
|
|
|
} |
|
125
|
|
|
|
|
126
|
2 |
|
protected function sqr(float $x): float |
|
127
|
|
|
{ |
|
128
|
2 |
|
return $x * $x; |
|
129
|
|
|
} |
|
130
|
|
|
|
|
131
|
2 |
|
public function toLonLat(Interfaces\IFormatted $source, array $params): Interfaces\INumbers |
|
132
|
|
|
{ |
|
133
|
2 |
|
$mul = isset($params['czech']) && boolval($params['czech']) ? 1 : -1; |
|
134
|
|
|
|
|
135
|
2 |
|
$sjtskLatLong = $this->toSJTSKLatLong($mul * floatval(strval($source->getLatitude())), $mul * floatval(strval($source->getLongitude()))); |
|
136
|
2 |
|
return $this->toWSG84LatLong($sjtskLatLong[0], $sjtskLatLong[1]); |
|
137
|
|
|
} |
|
138
|
|
|
|
|
139
|
|
|
/** |
|
140
|
|
|
* @param float $x |
|
141
|
|
|
* @param float $y |
|
142
|
|
|
* @return array<float> |
|
143
|
|
|
*/ |
|
144
|
2 |
|
protected function toSJTSKLatLong(float $x, float $y): array |
|
145
|
|
|
{ |
|
146
|
2 |
|
$finRad = pi() * 49.5 / 180; |
|
147
|
2 |
|
$la0deg = 42 + 31.00 / 60 + 31.41725 / 3600; |
|
148
|
2 |
|
$ro1 = sqrt($x * $x + $y * $y); |
|
149
|
2 |
|
$rGauss = 6380065.5402; |
|
150
|
2 |
|
$betakv = pi() * 11.5 / 180; |
|
151
|
2 |
|
$fi0 = 59 + (42.00 / 60) + (42.69689 / 3600); |
|
152
|
2 |
|
$a = 6377397.155; |
|
153
|
2 |
|
$b = 6356078.96325; |
|
154
|
|
|
|
|
155
|
2 |
|
$e2 = sqrt(($a * $a - $b * $b) / ($b * $b)); |
|
156
|
2 |
|
$m = cos($betakv); |
|
157
|
2 |
|
$vn = sqrt(1 + $e2 * $e2 * cos($finRad) * cos($finRad)); |
|
158
|
2 |
|
$gamma = atan2($y, $x); |
|
159
|
2 |
|
$c = $rGauss * sin($betakv) / $m / pow(tan($betakv / 2), $m); |
|
160
|
|
|
|
|
161
|
2 |
|
$fi0rad = $fi0 * pi() / 180; |
|
162
|
2 |
|
$fin = atan2(tan($finRad), $vn); |
|
163
|
2 |
|
$lambv = $gamma / $m; |
|
164
|
2 |
|
$betav = 2 * atan(pow($ro1 / $c, 1 / $m)); |
|
165
|
2 |
|
$fi = asin(cos($betav) * sin($fi0rad) - sin($betav) * cos($fi0rad) * cos($lambv)); |
|
166
|
|
|
|
|
167
|
2 |
|
$n = sin($finRad) / sin($fin); |
|
168
|
2 |
|
$la = asin(sin($betav) * sin($lambv) / cos($fi)); |
|
169
|
2 |
|
$la0 = $la0deg * pi() / 180; |
|
170
|
|
|
|
|
171
|
2 |
|
$laRad = ($la0 - $la) / $n; |
|
172
|
2 |
|
$laFerro = $laRad * 180 / pi(); |
|
173
|
2 |
|
$sjtskLong = $laFerro - static::FG_SHIFT; |
|
174
|
2 |
|
$e1 = sqrt(($a * $a - $b * $b) / ($a * $a)); |
|
175
|
|
|
|
|
176
|
2 |
|
$k = tan(pi() / 4 + $fin / 2) / (pow(tan(pi() / 4 + $finRad / 2), $n) * pow((1 - $e1 * sin($finRad)) / (1 + $e1 * sin($finRad)), $n * $e1 / 2)); |
|
177
|
2 |
|
$fik = 2 * atan(pow(1 / $k * tan(pi() / 4 + $fi / 2), 1 / $n)) - pi() / 2; |
|
178
|
2 |
|
$tga = pow((1 - $e1 * sin($fik)) / (1 + $e1 * sin($fik)), $n * $e1 / 2); |
|
179
|
2 |
|
$fi1Rad = 2 * atan(pow(1 / $k / $tga * tan(pi() / 4 + $fi / 2), 1 / $n)) - pi() / 2; |
|
180
|
2 |
|
$tga2 = pow((1 - $e1 * sin($fi1Rad)) / (1 + $e1 * sin($fi1Rad)), $n * $e1 / 2); |
|
181
|
2 |
|
$fi2Rad = 2 * atan(pow(1 / $k / $tga2 * tan(pi() / 4 + $fi / 2), 1 / $n)) - pi() / 2; |
|
182
|
2 |
|
$tga3 = pow((1 - $e1 * sin($fi2Rad)) / (1 + $e1 * sin($fi2Rad)), $n * $e1 / 2); |
|
183
|
2 |
|
$fi3Rad = 2 * atan(pow(1 / $k / $tga3 * tan(pi() / 4 + $fi / 2), 1 / $n)) - pi() / 2; |
|
184
|
2 |
|
$tga4 = pow((1 - $e1 * sin($fi3Rad)) / (1 + $e1 * sin($fi3Rad)), $n * $e1 / 2); |
|
185
|
2 |
|
$fi4Rad = 2 * atan(pow(1 / $k / $tga4 * tan(pi() / 4 + $fi / 2), 1 / $n)) - pi() / 2; |
|
186
|
2 |
|
$tga5 = pow((1 - $e1 * sin($fi4Rad)) / (1 + $e1 * sin($fi4Rad)), $n * $e1 / 2); |
|
187
|
2 |
|
$firad = 2 * atan(pow(1 / $k / $tga5 * tan(pi() / 4 + $fi / 2), 1 / $n)) - pi() / 2; |
|
188
|
2 |
|
$sjtskLat = $firad * 180 / pi(); |
|
189
|
|
|
|
|
190
|
2 |
|
return [$sjtskLat, $sjtskLong]; |
|
191
|
|
|
} |
|
192
|
|
|
|
|
193
|
2 |
|
protected function toWSG84LatLong(float $SJTSKlat, float $SJTSKlong): Interfaces\INumbers |
|
194
|
|
|
{ |
|
195
|
2 |
|
$filRad = $SJTSKlat * pi() / 180; |
|
196
|
2 |
|
$la1Rad = $SJTSKlong * pi() / 180; |
|
197
|
2 |
|
$a1 = 6377397.155; |
|
198
|
2 |
|
$f1 = 0.00334277318217481; |
|
199
|
2 |
|
$a2 = 6378137.00; |
|
200
|
2 |
|
$f2 = 0.00335281066474748; |
|
201
|
2 |
|
$e1 = sqrt(2 * $f1 - $f1 * $f1); |
|
202
|
2 |
|
$M = $a1 * (1 - $e1 * $e1) / pow((1 - $e1 * $e1 * sin($filRad) * sin($filRad)), 1.5); |
|
203
|
2 |
|
$N = $a1 / sqrt(1 - $e1 * $e1 * sin($filRad) * sin($filRad)); |
|
204
|
2 |
|
$dlasec = ((-1 * static::dX_SHIFT * sin($la1Rad)) + static::dY_SHIFT * cos($la1Rad)) / ($N * cos($filRad) * sin(pi() / 180 / 3600)); |
|
205
|
2 |
|
$dfisec = ((-1 * static::dX_SHIFT * sin($filRad) * cos($la1Rad)) - static::dY_SHIFT * sin($filRad) * sin($la1Rad) + static::dZ_SHIFT * cos($filRad) + ($a1 * ($f2 - $f1) + $f1 * ($a2 - $a1)) * sin(2 * $filRad)) / ($M * sin(pi() / 180 / 3600)); |
|
206
|
2 |
|
$wsgLong = $SJTSKlong + $dlasec / 3600; |
|
207
|
2 |
|
$wsgLat = $SJTSKlat + $dfisec / 3600; |
|
208
|
2 |
|
return (clone $this->pos)->setData($wsgLong, $wsgLat); |
|
209
|
|
|
} |
|
210
|
|
|
} |
|
211
|
|
|
|
alex-kalanis /
kw_coordinates
Loading history...