earth.inc 6.7 KB

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