TransverseMercator::getOriginLatitude()

Completed

Branch — 3.0-dev (e499db)

by Doug

created 2017-10-29 15:49 UTC

↳ Parent: TransverseMercator

Size

Total Lines

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
dl	0
loc	1
ccs	0
cts	0
cp	0
c	0
b	0
f	0
nc	1

<?php
/**
 * PHPCoord
 * @package PHPCoord
 * @author Jonathan Stott
 * @author Doug Wright
 */
declare(strict_types=1);
namespace PHPCoord;

/**
 * Abstract class representing a Tranverse Mercator Projection
 * @author Doug Wright
 * @package PHPCoord
 */
abstract class TransverseMercator
{

    /**
     * X
     * @var int
     */
    protected $x;

    /**
     * Y
     * @var int
     */
    protected $y;

    /**
     * h
     * @var int
     */
    protected $h;

    /**
     * Reference ellipsoid used in this datum
     * @var RefEll
     */
    protected $refEll;

    /**
     * Cartesian constructor.
     * @param int $x
     * @param int $y
     * @param int $h
     * @param RefEll $refEll
     */
    public function __construct(int $x, int $y, int $h, RefEll $refEll)
    {
        $this->x = $x;
        $this->y = $y;
        $this->h = $h;
        $this->refEll = $refEll;
    }

    /**
     * String version of coordinate.
     * @return string
     */
    public function __toString(): string
    {
        return "({$this->x}, {$this->y}, {$this->h})";
    }

    /**
     * @return int
     */
    public function getX(): int
    {
        return $this->x;
    }

    /**
     * @return int
     */
    public function getY(): int
    {
        return $this->y;
    }

    /**
     * @return int
     */
    public function getH(): int
    {
        return $this->h;
    }

    /**
     * Reference ellipsoid used by this projection
     * @return RefEll
     */
    abstract public function getReferenceEllipsoid(): RefEll;

    /**
     * Scale factor at central meridian
     * @return float
     */
    abstract public function getScaleFactor(): float;

    /**
     * Northing of true origin
     * @return int
     */
    abstract public function getOriginNorthing(): int;

    /**
     * Easting of true origin
     * @return int
     */
    abstract public function getOriginEasting(): int;

    /**
     * Latitude of true origin
     * @return float
     */
    abstract public function getOriginLatitude(): float;

    /**
     * Longitude of true origin
     * @return float
     */
    abstract public function getOriginLongitude(): float;

    /**
     * Convert this grid reference into a latitude and longitude
     * Formula for transformation is taken from OS document
     * "A Guide to Coordinate Systems in Great Britain"
     *
     * @param float $N map coordinate (northing) of point to convert
     * @param float $E map coordinate (easting) of point to convert
     * @param float $N0 map coordinate (northing) of true origin
     * @param float $E0 map coordinate (easting) of true origin
     * @param float $phi0 map coordinate (latitude) of true origin
     * @param float $lambda0 map coordinate (longitude) of true origin and central meridian
     * @return LatLng
     */
    public function convertToLatitudeLongitude($N, $E, $N0, $E0, $phi0, $lambda0): LatLng
    {

        $phi0 = deg2rad($phi0);
        $lambda0 = deg2rad($lambda0);

        $refEll = $this->getReferenceEllipsoid();
        $F0 = $this->getScaleFactor();

        $a = $refEll->getMaj();
        $b = $refEll->getMin();
        $eSquared = $refEll->getEcc();
        $n = ($a - $b) / ($a + $b);
        $phiPrime = (($N - $N0) / ($a * $F0)) + $phi0;

        do {
            $M =
                ($b * $F0)
                * (((1 + $n + ((5 / 4) * $n * $n) + ((5 / 4) * $n * $n * $n))
                        * ($phiPrime - $phi0))
                    - (((3 * $n) + (3 * $n * $n) + ((21 / 8) * $n * $n * $n))
                        * sin($phiPrime - $phi0)
                        * cos($phiPrime + $phi0))
                    + ((((15 / 8) * $n * $n) + ((15 / 8) * $n * $n * $n))
                        * sin(2 * ($phiPrime - $phi0))
                        * cos(2 * ($phiPrime + $phi0)))
                    - (((35 / 24) * $n * $n * $n)
                        * sin(3 * ($phiPrime - $phi0))
                        * cos(3 * ($phiPrime + $phi0))));
            $phiPrime += ($N - $N0 - $M) / ($a * $F0);
        } while (($N - $N0 - $M) >= 0.00001);
        $v = $a * $F0 * pow(1 - $eSquared * pow(sin($phiPrime), 2), -0.5);
        $rho =
            $a
            * $F0
            * (1 - $eSquared)
            * pow(1 - $eSquared * pow(sin($phiPrime), 2), -1.5);
        $etaSquared = ($v / $rho) - 1;
        $VII = tan($phiPrime) / (2 * $rho * $v);
        $VIII =
            (tan($phiPrime) / (24 * $rho * pow($v, 3)))
            * (5
                + (3 * pow(tan($phiPrime), 2))
                + $etaSquared
                - (9 * pow(tan($phiPrime), 2) * $etaSquared));
        $IX =
            (tan($phiPrime) / (720 * $rho * pow($v, 5)))
            * (61
                + (90 * pow(tan($phiPrime), 2))
                + (45 * pow(tan($phiPrime), 2) * pow(tan($phiPrime), 2)));
        $X = (1 / cos($phiPrime)) / $v;
        $XI =
            ((1 / cos($phiPrime)) / (6 * $v * $v * $v))
            * (($v / $rho) + (2 * pow(tan($phiPrime), 2)));
        $XII =
            ((1 / cos($phiPrime)) / (120 * pow($v, 5)))
            * (5
                + (28 * pow(tan($phiPrime), 2))
                + (24 * pow(tan($phiPrime), 4)));
        $XIIA =
            ((1 / cos($phiPrime)) / (5040 * pow($v, 7)))
            * (61
                + (662 * pow(tan($phiPrime), 2))
                + (1320 * pow(tan($phiPrime), 4))
                + (720
                    * pow(tan($phiPrime), 6)));
        $phi =
            $phiPrime
            - ($VII * pow($E - $E0, 2))
            + ($VIII * pow($E - $E0, 4))
            - ($IX * pow($E - $E0, 6));
        $lambda =
            $lambda0
            + ($X * ($E - $E0))
            - ($XI * pow($E - $E0, 3))
            + ($XII * pow($E - $E0, 5))
            - ($XIIA * pow($E - $E0, 7));

        return new LatLng(rad2deg($phi), rad2deg($lambda), 0, $refEll);
    }

    /**
     * Calculate the surface distance between this object and the one
     * passed in as a parameter.
     *
     * @param self $to object to measure the distance to

     * @return int
     */
    public function distance(self $to): int {

        if ($this->refEll != $to->refEll) {
            throw new \RuntimeException('Source and destination co-ordinates are not using the same ellipsoid');
        }

        //Because this is a 2D grid, we can use simple Pythagoras
        $distanceX = $to->getX()-$this->getX();
        $distanceY = $to->getY()-$this->getY();

        return (int) round(pow(pow($distanceX, 2) + pow($distanceY, 2), 0.5));

    }
}


1		<?php
2		/**
3		* PHPCoord
4		* @package PHPCoord
5		* @author Jonathan Stott
6		* @author Doug Wright
7		*/
8		declare(strict_types=1);
9		namespace PHPCoord;
10
11		/**
12		* Abstract class representing a Tranverse Mercator Projection
13		* @author Doug Wright
14		* @package PHPCoord
15		*/
16		abstract class TransverseMercator
17		{
18
19		/**
20		* X
21		* @var int
22		*/
23		protected $x;
24
25		/**
26		* Y
27		* @var int
28		*/
29		protected $y;
30
31		/**
32		* h
33		* @var int
34		*/
35		protected $h;
36
37		/**
38		* Reference ellipsoid used in this datum
39		* @var RefEll
40		*/
41		protected $refEll;
42
43		/**
44		* Cartesian constructor.
45		* @param int $x
46		* @param int $y
47		* @param int $h
48		* @param RefEll $refEll
49		*/
50	31	public function __construct(int $x, int $y, int $h, RefEll $refEll)
51		{
52	31	$this->x = $x;
53	31	$this->y = $y;
54	31	$this->h = $h;
55	31	$this->refEll = $refEll;
56		}
57
58		/**
59		* String version of coordinate.
60		* @return string
61		*/
62		public function __toString(): string
63		{
64		return "({$this->x}, {$this->y}, {$this->h})";
65		}
66
67		/**
68		* @return int
69		*/
70	2	public function getX(): int
71		{
72	2	return $this->x;
73		}
74
75		/**
76		* @return int
77		*/
78	2	public function getY(): int
79		{
80	2	return $this->y;
81		}
82
83		/**
84		* @return int
85		*/
86	1	public function getH(): int
87		{
88	1	return $this->h;
89		}
90
91		/**
92		* Reference ellipsoid used by this projection
93		* @return RefEll
94		*/
95		abstract public function getReferenceEllipsoid(): RefEll;
96
97		/**
98		* Scale factor at central meridian
99		* @return float
100		*/
101		abstract public function getScaleFactor(): float;
102
103		/**
104		* Northing of true origin
105		* @return int
106		*/
107		abstract public function getOriginNorthing(): int;
108
109		/**
110		* Easting of true origin
111		* @return int
112		*/
113		abstract public function getOriginEasting(): int;
114
115		/**
116		* Latitude of true origin
117		* @return float
118		*/
119		abstract public function getOriginLatitude(): float;
120
121		/**
122		* Longitude of true origin
123		* @return float
124		*/
125		abstract public function getOriginLongitude(): float;
126
127		/**
128		* Convert this grid reference into a latitude and longitude
129		* Formula for transformation is taken from OS document
130		* "A Guide to Coordinate Systems in Great Britain"
131		*
132		* @param float $N map coordinate (northing) of point to convert
133		* @param float $E map coordinate (easting) of point to convert
134		* @param float $N0 map coordinate (northing) of true origin
135		* @param float $E0 map coordinate (easting) of true origin
136		* @param float $phi0 map coordinate (latitude) of true origin
137		* @param float $lambda0 map coordinate (longitude) of true origin and central meridian
138		* @return LatLng
139		*/
140	7	public function convertToLatitudeLongitude($N, $E, $N0, $E0, $phi0, $lambda0): LatLng
141		{
142
143	7	$phi0 = deg2rad($phi0);
144	7	$lambda0 = deg2rad($lambda0);
145
146	7	$refEll = $this->getReferenceEllipsoid();
147	7	$F0 = $this->getScaleFactor();
148
149	7	$a = $refEll->getMaj();
150	7	$b = $refEll->getMin();
151	7	$eSquared = $refEll->getEcc();
152	7	$n = ($a - $b) / ($a + $b);
153	7	$phiPrime = (($N - $N0) / ($a * $F0)) + $phi0;
154
155		do {
156		$M =
157	7	($b * $F0)
158	7	* (((1 + $n + ((5 / 4) * $n * $n) + ((5 / 4) * $n * $n * $n))
159	7	* ($phiPrime - $phi0))
160	7	- (((3 * $n) + (3 * $n * $n) + ((21 / 8) * $n * $n * $n))
161	7	* sin($phiPrime - $phi0)
162	7	* cos($phiPrime + $phi0))
163	7	+ ((((15 / 8) * $n * $n) + ((15 / 8) * $n * $n * $n))
164	7	* sin(2 * ($phiPrime - $phi0))
165	7	* cos(2 * ($phiPrime + $phi0)))
166	7	- (((35 / 24) * $n * $n * $n)
167	7	* sin(3 * ($phiPrime - $phi0))
168	7	* cos(3 * ($phiPrime + $phi0))));
169	7	$phiPrime += ($N - $N0 - $M) / ($a * $F0);
170	7	} while (($N - $N0 - $M) >= 0.00001);
171	7	$v = $a * $F0 * pow(1 - $eSquared * pow(sin($phiPrime), 2), -0.5);
172		$rho =
173		$a
174	7	* $F0
175	7	* (1 - $eSquared)
176	7	* pow(1 - $eSquared * pow(sin($phiPrime), 2), -1.5);
177	7	$etaSquared = ($v / $rho) - 1;
178	7	$VII = tan($phiPrime) / (2 * $rho * $v);
179		$VIII =
180	7	(tan($phiPrime) / (24 * $rho * pow($v, 3)))
181		* (5
182	7	+ (3 * pow(tan($phiPrime), 2))
183	7	+ $etaSquared
184	7	- (9 * pow(tan($phiPrime), 2) * $etaSquared));
185		$IX =
186	7	(tan($phiPrime) / (720 * $rho * pow($v, 5)))
187		* (61
188	7	+ (90 * pow(tan($phiPrime), 2))
189	7	+ (45 * pow(tan($phiPrime), 2) * pow(tan($phiPrime), 2)));
190	7	$X = (1 / cos($phiPrime)) / $v;
191		$XI =
192	7	((1 / cos($phiPrime)) / (6 * $v * $v * $v))
193	7	* (($v / $rho) + (2 * pow(tan($phiPrime), 2)));
194		$XII =
195	7	((1 / cos($phiPrime)) / (120 * pow($v, 5)))
196		* (5
197	7	+ (28 * pow(tan($phiPrime), 2))
198	7	+ (24 * pow(tan($phiPrime), 4)));
199		$XIIA =
200	7	((1 / cos($phiPrime)) / (5040 * pow($v, 7)))
201		* (61
202	7	+ (662 * pow(tan($phiPrime), 2))
203	7	+ (1320 * pow(tan($phiPrime), 4))
204		+ (720
205	7	* pow(tan($phiPrime), 6)));
206		$phi =
207		$phiPrime
208	7	- ($VII * pow($E - $E0, 2))
209	7	+ ($VIII * pow($E - $E0, 4))
210	7	- ($IX * pow($E - $E0, 6));
211		$lambda =
212		$lambda0
213	7	+ ($X * ($E - $E0))
214	7	- ($XI * pow($E - $E0, 3))
215	7	+ ($XII * pow($E - $E0, 5))
216	7	- ($XIIA * pow($E - $E0, 7));
217
218	7	return new LatLng(rad2deg($phi), rad2deg($lambda), 0, $refEll);
219		}
220
221		/**
222		* Calculate the surface distance between this object and the one
223		* passed in as a parameter.
224		*
225		* @param self $to object to measure the distance to
		0 ignored issues – show Documentation introduced 2017-01-22 19:24 UTC by Report Bug Copy Issue Report Should the type for parameter `$to` not be `\self`? This check looks for `@param` annotations where the type inferred by our type inference engine differs from the declared type. It makes a suggestion as to what type it considers more descriptive. Most often this is a case of a parameter that can be null in addition to its declared types. Loading history...
226		* @return int
227		*/
228	1	public function distance(self $to): int {
229
230	1	if ($this->refEll != $to->refEll) {
231		throw new \RuntimeException('Source and destination co-ordinates are not using the same ellipsoid');
232		}
233
234		//Because this is a 2D grid, we can use simple Pythagoras
235	1	$distanceX = $to->getX()-$this->getX();
236	1	$distanceY = $to->getY()-$this->getY();
237
238	1	return (int) round(pow(pow($distanceX, 2) + pow($distanceY, 2), 0.5));
239
240		}
241		}
242

dvdoug / PHPCoord

Branch — 3.0-dev (e499db)

TransverseMercator::getOriginLatitude()

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like