-
Notifications
You must be signed in to change notification settings - Fork 0
/
ringArea.php
64 lines (60 loc) · 1.82 KB
/
ringArea.php
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
<?php
/**
* Calculate the approximate area of the polygon were it projected onto
* the earth. Note that this area will be positive if ring is oriented
* clockwise, otherwise it will be negative.
*
* Reference:
* Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for
* Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
* Laboratory, Pasadena, CA, June 2007 http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409
*
* Returns:
* {float} The approximate signed geodesic area of the polygon in square
* kilometers.
*/
function ringArea( $coords ) {
$RADIUS = 6378137;
$p1 = null;
$p2 = null;
$p3 = null;
$lowerIndex = null;
$middleIndex = null;
$upperIndex = null;
$i = null;
$total = 0;
$coordsLength = count( $coords );
if ( $coordsLength > 2 ) {
for ( $i = 0; $i < $coordsLength; $i++ ) {
if ( $i === $coordsLength - 2 ) {
// i = N-2
$lowerIndex = $coordsLength - 2;
$middleIndex = $coordsLength - 1;
$upperIndex = 0;
} else
if ( $i === $coordsLength - 1 ) {
// i = N-1
$lowerIndex = $coordsLength - 1;
$middleIndex = 0;
$upperIndex = 1;
} else {
// i = 0 to N-3
$lowerIndex = $i;
$middleIndex = $i + 1;
$upperIndex = $i + 2;
}
$p1 = $coords[ $lowerIndex ];
$p2 = $coords[ $middleIndex ];
$p3 = $coords[ $upperIndex ];
$total += ( rad( $p3[ 0 ] ) - rad( $p1[ 0 ] ) ) * sin( rad( $p2[ 1 ] ) );
}
$total = ( $total * $RADIUS * $RADIUS ) / 2;
}
//return the total area in SQ KMS and abs will remove the minus (from numbers)
return abs($total) / 1000000;
}
//helper for ringArea function
function rad( $num ) {
return ( $num * M_PI ) / 180;
}
?>