earth.inc 7.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243
  1. <?php
  2. /**
  3. * @file
  4. * Trigonometry for calculating geographical distances.
  5. * All function arguments and return values measure distances in metres
  6. * and angles in degrees. The ellipsoid model is from the WGS-84 datum.
  7. * Ka-Ping Yee, 2003-08-11
  8. */
  9. /**
  10. * License clarification:
  11. *
  12. * On Feb 13, 2005, in message <Pine.LNX.4.58.0502131827510.5072@server1.LFW.org>,
  13. * the creator of these routines, Ka-Ping Yee, authorized these routines to be
  14. * distributed under the GPL.
  15. */
  16. // @codingStandardsIgnoreStart
  17. // @todo cleanup.
  18. //$earth_radius_semimajor = 6378137.0;
  19. //$earth_flattening = 1/298.257223563;
  20. //$earth_radius_semiminor = $earth_radius_semimajor * (1 - $earth_flattening);
  21. //$earth_eccentricity_sq = 2*$earth_flattening - pow($earth_flattening, 2);
  22. // I don't know what's up: PHP is hating on my global variables (commented out above),
  23. // so I have to write functions that return them! (-Ankur)
  24. // Commenting out the global variables above and replacing them with functions that
  25. // return the same values is the only thing I changed since, for some reason, my
  26. // PHP wasn't acknowledging these global variables.
  27. // This library is an original implementation of UCB CS graduate student, Ka-Ping Yee (http://www.zesty.ca).
  28. // @codingStandardsIgnoreEnd
  29. /**
  30. * Earth redius semimajor.
  31. */
  32. function earth_radius_semimajor() {
  33. return 6378137.0;
  34. }
  35. /**
  36. * Earth flattening.
  37. */
  38. function earth_flattening() {
  39. return (1 / 298.257223563);
  40. }
  41. /**
  42. * Earth radius semiminor.
  43. */
  44. function earth_radius_semiminor() {
  45. return (earth_radius_semimajor() * (1 - earth_flattening()));
  46. }
  47. /**
  48. * Earth ccentricity sq.
  49. */
  50. function earth_eccentricity_sq() {
  51. return (2 * earth_flattening() - pow(earth_flattening(), 2));
  52. }
  53. /**
  54. * Earth redius.
  55. *
  56. * Latitudes in all of U. S.: from -7.2 (American Samoa) to 70.5 (Alaska).
  57. * Latitudes in continental U. S.: from 24.6 (Florida) to 49.0 (Washington).
  58. * Average latitude of all U. S. zipcodes: 37.9.
  59. */
  60. function earth_radius($latitude = 37.9) {
  61. // Estimate the Earth's radius at a given latitude.
  62. // Default to an approximate average radius for the United States.
  63. $lat = deg2rad($latitude);
  64. $x = cos($lat) / earth_radius_semimajor();
  65. $y = sin($lat) / earth_radius_semiminor();
  66. return 1 / (sqrt($x * $x + $y * $y));
  67. }
  68. /**
  69. * Earth xyz.
  70. */
  71. function earth_xyz($longitude, $latitude, $height = 0) {
  72. // Convert longitude and latitude to earth-centered earth-fixed coordinates.
  73. // X axis is 0 long, 0 lat; Y axis is 90 deg E; Z axis is north pole.
  74. $long = deg2rad($longitude);
  75. $lat = deg2rad($latitude);
  76. $coslong = cos($long);
  77. $coslat = cos($lat);
  78. $sinlong = sin($long);
  79. $sinlat = sin($lat);
  80. $radius = earth_radius_semimajor() /
  81. sqrt(1 - earth_eccentricity_sq() * $sinlat * $sinlat);
  82. $x = ($radius + $height) * $coslat * $coslong;
  83. $y = ($radius + $height) * $coslat * $sinlong;
  84. $z = ($radius * (1 - earth_eccentricity_sq()) + $height) * $sinlat;
  85. return array($x, $y, $z);
  86. }
  87. /**
  88. * Earth arclength.
  89. */
  90. function earth_arclength($angle, $latitude = 37.9) {
  91. // Convert a given angle to earth-surface distance.
  92. return deg2rad($angle) * earth_radius($latitude);
  93. }
  94. /**
  95. * Earth distance.
  96. */
  97. function earth_distance($longitude1, $latitude1, $longitude2, $latitude2) {
  98. // Estimate the earth-surface distance between two locations.
  99. $long1 = deg2rad($longitude1);
  100. $lat1 = deg2rad($latitude1);
  101. $long2 = deg2rad($longitude2);
  102. $lat2 = deg2rad($latitude2);
  103. $radius = earth_radius(($latitude1 + $latitude2) / 2);
  104. $cosangle = cos($lat1) * cos($lat2) *
  105. (cos($long1) * cos($long2) + sin($long1) * sin($long2)) +
  106. sin($lat1) * sin($lat2);
  107. return acos($cosangle) * $radius;
  108. }
  109. /**
  110. * Returns the SQL fragment needed to add a column called 'distance' to a query that includes the location table.
  111. *
  112. * @param float $longitude
  113. * The measurement point.
  114. *
  115. * @param float $latitude
  116. * The measurement point.
  117. *
  118. * @param string $tbl_alias
  119. * If necessary, the alias name of the location table to work from. Only required when working with named {location} tables.
  120. *
  121. * @return string
  122. * SQL fragment.
  123. */
  124. function earth_distance_sql($longitude, $latitude, $tbl_alias = '') {
  125. // Make a SQL expression that estimates the distance to the given location.
  126. $long = deg2rad($longitude);
  127. $lat = deg2rad($latitude);
  128. $radius = earth_radius($latitude);
  129. // If the table alias is specified, add on the separator.
  130. $tbl_alias = empty($tbl_alias) ? $tbl_alias : ($tbl_alias . '.');
  131. $coslong = cos($long);
  132. $coslat = cos($lat);
  133. $sinlong = sin($long);
  134. $sinlat = sin($lat);
  135. return "(COALESCE(ACOS($coslat*COS(RADIANS({$tbl_alias}latitude))*($coslong*COS(RADIANS({$tbl_alias}longitude)) + $sinlong*SIN(RADIANS({$tbl_alias}longitude))) + $sinlat*SIN(RADIANS({$tbl_alias}latitude))), 0.00000)*$radius)";
  136. }
  137. /**
  138. * Earth longitude range.
  139. *
  140. * @todo This function uses earth_asin_safe so is not accurate for all possible
  141. * parameter combinations. This means this function doesn't work properly
  142. * for high distance values. This function needs to be re-written to work properly for
  143. * larger distance values.
  144. *
  145. * @see http://drupal.org/node/821628
  146. */
  147. function earth_longitude_range($longitude, $latitude, $distance) {
  148. // Estimate the min and max longitudes within $distance of a given location.
  149. $long = deg2rad($longitude);
  150. $lat = deg2rad($latitude);
  151. $radius = earth_radius($latitude) * cos($lat);
  152. if ($radius > 0) {
  153. $angle = abs($distance / $radius);
  154. $angle = min($angle, pi());
  155. }
  156. else {
  157. $angle = pi();
  158. }
  159. $minlong = $long - $angle;
  160. $maxlong = $long + $angle;
  161. if ($minlong < -pi()) {
  162. $minlong = $minlong + pi() * 2;
  163. }
  164. if ($maxlong > pi()) {
  165. $maxlong = $maxlong - pi() * 2;
  166. }
  167. return array(rad2deg($minlong), rad2deg($maxlong));
  168. }
  169. /**
  170. * Earth latitude range.
  171. */
  172. function earth_latitude_range($longitude, $latitude, $distance) {
  173. // Estimate the min and max latitudes within $distance of a given location.
  174. $long = deg2rad($longitude);
  175. $lat = deg2rad($latitude);
  176. $radius = earth_radius($latitude);
  177. $angle = $distance / $radius;
  178. $minlat = $lat - $angle;
  179. $maxlat = $lat + $angle;
  180. $rightangle = pi() / 2;
  181. // Wrapped around the south pole.
  182. if ($minlat < -$rightangle) {
  183. $overshoot = -$minlat - $rightangle;
  184. $minlat = -$rightangle + $overshoot;
  185. if ($minlat > $maxlat) {
  186. $maxlat = $minlat;
  187. }
  188. $minlat = -$rightangle;
  189. }
  190. // Wrapped around the north pole.
  191. if ($maxlat > $rightangle) {
  192. $overshoot = $maxlat - $rightangle;
  193. $maxlat = $rightangle - $overshoot;
  194. if ($maxlat < $minlat) {
  195. $minlat = $maxlat;
  196. }
  197. $maxlat = $rightangle;
  198. }
  199. return array(rad2deg($minlat), rad2deg($maxlat));
  200. }
  201. /**
  202. * This is a helper function to avoid errors when using the asin() PHP function.
  203. *
  204. * asin is only real for values between -1 and 1.
  205. * If a value outside that range is given it returns NAN (not a number), which
  206. * we don't want to happen.
  207. * So this just rounds values outside this range to -1 or 1 first.
  208. *
  209. * This means that calculations done using this function with $x outside the range
  210. * will not be accurate. The alternative though is getting NAN, which is an error
  211. * and won't give accurate results anyway.
  212. */
  213. function earth_asin_safe($x) {
  214. return asin(max(-1, min($x, 1)));
  215. }