Module: GeoCalc

Included in:
GeoPoint
Defined in:
lib/geo_calc/calculations.rb

Overview

                                                                                            • -

Latitude/longitude spherical geodesy formulae & scripts (c) Chris Veness 2002-2010            
 - www.movable-type.co.uk/scripts/latlong.html

Class Method Summary collapse

Instance Method Summary collapse

Class Method Details

.intersection(p1, brng1, p2, brng2) ⇒ Object

Returns the point of intersection of two paths defined by point and bearing

see http:#williams.best.vwh.net/avform.htm#Intersection

Parameters:

  • p1: (LatLon) —

    First point

  • brng1: (Number) —

    Initial bearing from first point

  • p2: (LatLon) —

    Second point

  • brng2: (Number) —

    Initial bearing from second point



140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
 
# File 'lib/geo_calc/calculations.rb', line 140

def self.intersection p1, brng1, p2, brng2
  lat1 = p1.lat.to_rad
  lon1 = p1.lon.to_rad

  lat2 = p2.lat.to_rad
  lon2 = p2.lon.to_rad

  brng13 = brng1.to_rad
  brng23 = brng2.to_rad
      
  dlat = lat2-lat1
  dlon = lon2-lon1;

  dist12 = 2*Math.asin( Math.sqrt( Math.sin(dlat/2)*Math.sin(dlat/2) + Math.cos(lat1)*Math.cos(lat2)*Math.sin(dlon/2)*Math.sin(dlon/2) ) )
  return nil if dist12 == 0

  # initial/final bearings between points
  brng_a = begin
    Math.acos( ( Math.sin(lat2) - Math.sin(lat1)*Math.cos(dist12) ) / ( Math.sin(dist12)*Math.cos(lat1) ) )
  rescue # protect against rounding
    0
  end
  
  brng_b = Math.acos( ( Math.sin(lat1) - Math.sin(lat2)*Math.cos(dist12) ) / ( Math.sin(dist12)*Math.cos(lat2) ) )

  brng12, brng21 = if Math.sin(lon2-lon1) > 0
    [brng_a, 2*Math::PI - brng_b] 
  else
    [2*Math::PI - brng_a, brng_b]
  end

  alpha1 = (brng13 - brng12 + Math::PI) % (2*Math::PI) - Math::PI  # angle 2-1-3
  alpha2 = (brng21 - brng23 + Math::PI) % (2*Math::PI) - Math::PI  # angle 1-2-3

  return nil if (Math.sin(alpha1)==0 && Math.sin(alpha2)==0) # infinite intersections
  return nil if (Math.sin(alpha1)*Math.sin(alpha2) < 0)    # ambiguous intersection

  # alpha1 = Math.abs(alpha1);
  # alpha2 = Math.abs(alpha2);
  # ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?

  alpha3 = Math.acos( -Math.cos(alpha1)*Math.cos(alpha2) + Math.sin(alpha1)*Math.sin(alpha2)*Math.cos(dist12) )

  dist13 = Math.atan2( Math.sin(dist12)*Math.sin(alpha1)*Math.sin(alpha2), Math.cos(alpha2)+Math.cos(alpha1)*Math.cos(alpha3) )

  lat3 = Math.asin( Math.sin(lat1)*Math.cos(dist13) + Math.cos(lat1)*Math.sin(dist13)*Math.cos(brng13) )

  dlon13 = Math.atan2( Math.sin(brng13)*Math.sin(dist13)*Math.cos(lat1), Math.cos(dist13)-Math.sin(lat1)*Math.sin(lat3) )

  lon3 = lon1 + dlon13;
  lon3 = (lon3 + Math::PI) % (2*Math::PI) - Math::PI  # normalise to -180..180º

  GeoPoint.new lat3.to_deg, lon3.to_deg
end

Instance Method Details

#bearing_to(point) ⇒ Object

Returns the (initial) bearing from this point to the supplied point, in degrees

see http:#williams.best.vwh.net/avform.htm#Crs

  • Point point: Latitude/longitude of destination point

Returns - Numeric: Initial bearing in degrees from North



48
49
50
51
52
53
54
55
56
57
58
 
# File 'lib/geo_calc/calculations.rb', line 48

def bearing_to point
  lat1 = lat.to_rad
  lat2 = point.lat.to_rad
  dlon = (point.lon - lon).to_rad

  y = Math.sin(dlon) * Math.cos(lat2)
  x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(dlon)
  bearing = Math.atan2(y, x)

  (bearing.to_deg + 360) % 360
end

#destination_point(brng, dist) ⇒ Object

Returns the destination point from this point having travelled the given distance (in km) on the given initial bearing (bearing may vary before destination is reached)

see http:#williams.best.vwh.net/avform.htm#LL

  • Numeric bearing: Initial bearing in degrees

  • Numeric dist: Distance in km

Returns GeoPoint: Destination point



115
116
117
118
119
120
121
122
123
124
125
126
127
 
# File 'lib/geo_calc/calculations.rb', line 115

def destination_point brng, dist
  dist = dist / radius  # convert dist to angular distance in radians
  brng = brng.to_rad 
  lat1 = lat.to_rad
  lon1 = lon.to_rad

  lat2 = Math.asin( Math.sin(lat1) * Math.cos(dist) + Math.cos(lat1) * Math.sin(dist) * Math.cos(brng) )
  lon2 = lon1 + Math.atan2(Math.sin(brng) * Math.sin(dist) * Math.cos(lat1), Math.cos(dist) - Math.sin(lat1) * Math.sin(lat2))

  lon2 = (lon2 + 3*Math::PI) % (2*Math::PI) - Math::PI  # normalise to -180...+180   

  GeoPoint.new lat2.to_deg, lon2.to_deg
end

#distance_to(point, precision = 4) ⇒ Object

Returns the distance from this point to the supplied point, in km (using Haversine formula)

from: Haversine formula - R. W. Sinnott, “Virtues of the Haversine”,

Sky and Telescope, vol 68, no 2, 1984

GeoPoint point: Latitude/longitude of destination point

  • Numeric precision=4: number of significant digits to use for returned value

Returns - Numeric distance in km between this point and destination point



21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
 
# File 'lib/geo_calc/calculations.rb', line 21

def distance_to point, precision = 4
  # default 4 sig figs reflects typical 0.3% accuracy of spherical model
  precision ||= 4

  lat1 = lat.to_rad
  lon1 = lon.to_rad

  lat2 = point.lat.to_rad
  lon2 = point.lon.to_rad

  dlat = lat2 - lat1
  dlon = lon2 - lon1

  a = Math.sin(dlat/2) * Math.sin(dlat/2) + Math.cos(lat1) * Math.cos(lat2) * Math.sin(dlon/2) * Math.sin(dlon/2)
  c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a))
  d = radius * c
  d.round(precision)
end

#final_bearing_to(point) ⇒ Object

Returns final bearing arriving at supplied destination point from this point; the final bearing will differ from the initial bearing by varying degrees according to distance and latitude

  • GeoPoint point: Latitude/longitude of destination point

Returns Numeric: Final bearing in degrees from North



68
69
70
71
72
73
74
75
76
77
78
79
80
 
# File 'lib/geo_calc/calculations.rb', line 68

def final_bearing_to point
  # get initial bearing from supplied point back to this point...
  lat1 = point.lat.to_rad
  lat2 = lat.to_rad
  dlon = (lon - point.lon).to_rad

  y = Math.sin(dlon) * Math.cos(lat2)
  x = Math.cos(lat1)*Math.sin(lat2) - Math.sin(lat1)*Math.cos(lat2)*Math.cos(dlon)
  bearing = Math.atan2(y, x)

  # ... & reverse it by adding 180°
  (bearing.to_deg+180) % 360
end

#midpoint_to(point) ⇒ Object

Returns the midpoint between this point and the supplied point.

see http:#mathforum.org/library/drmath/view/51822.html for derivation

  • GeoPoint point: Latitude/longitude of destination point

Returns GeoPoint: Midpoint between this point and the supplied point



89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
 
# File 'lib/geo_calc/calculations.rb', line 89

def midpoint_to point
  lat1 = lat.to_rad
  lon1 = lon.to_rad;
  lat2 = point.lat.to_rad
  dlon = (point.lon - lon).to_rad

  bx = Math.cos(lat2) * Math.cos(dlon)
  by = Math.cos(lat2) * Math.sin(dlon)

  lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2), Math.sqrt( (Math.cos(lat1)+bx)*(Math.cos(lat1)+bx) + by*by) )

  lon3 = lon1 + Math.atan2(by, Math.cos(lat1) + bx)

  GeoPoint.new lat3.to_deg, lon3.to_deg
end

#rhumb_bearing_to(point) ⇒ Object

Returns the bearing from this point to the supplied point along a rhumb line, in degrees

  • GeoPoint point: Latitude/longitude of destination point

Returns Numeric: Bearing in degrees from North



232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
 
# File 'lib/geo_calc/calculations.rb', line 232

def rhumb_bearing_to point
  lat1 = lat.to_rad
  lat2 = point.lat.to_rad

  dlon = (point.lon - lon).to_rad

  dphi = Math.log(Math.tan(lat2/2+Math::PI/4) / Math.tan(lat1/2+Math::PI/4))
  if dlon.abs > Math::PI
    dlon = dlon>0 ? -(2*Math::PI-dlon) : (2*Math::PI+dlon);
  end

  brng = Math.atan2(dlon, dphi);

  (brng.to_deg+360) % 360
end

#rhumb_destination_point(brng, dist) ⇒ Object

Returns the destination point from this point having travelled the given distance (in km) on the given bearing along a rhumb line

Parameters:

  • brng: (Number) —

    Bearing in degrees from North

  • dist: (Number) —

    Distance in km



256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
 
# File 'lib/geo_calc/calculations.rb', line 256

def rhumb_destination_point brng, dist
  d = dist / radius  # d = angular distance covered on earth's surface
  lat1 = lat.to_rad
  lon1 = lon.to_rad
  brng = brng.to_rad

  lat2 = lat1 + d*Math.cos(brng);
  dlat = lat2-lat1;
  dphi = Math.log(Math.tan(lat2/2+Math::PI/4)/Math.tan(lat1/2+Math::PI/4))

  q = begin
    dlat/dphi
  rescue 
    Math.cos(lat1) # E-W line gives dPhi=0
  end

  dlon = d*Math.sin(brng)/q
  # check for some daft bugger going past the pole

  if lat2.abs > Math::PI/2
    lat2 = lat2>0 ? Math::PI-lat2 : -(Math::PI-lat2)
  end
  lon2 = (lon1+dlon+3*Math::PI) % (2*Math::PI) - Math::PI

  GeoPoint.new lat2.to_deg, lon2.to_deg
end

#rhumb_distance_to(point) ⇒ Object

Returns the distance from this point to the supplied point, in km, travelling along a rhumb line

see http:#williams.best.vwh.net/avform.htm#Rhumb

  • GeoPoint point: Latitude/longitude of destination point

Returns Numeric: Distance in km between this point and destination point



203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
 
# File 'lib/geo_calc/calculations.rb', line 203

def rhumb_distance_to point
  lat1 = lat.to_rad
  lat2 = point.lat.to_rad

  dlat = (point.lat-lat).to_rad
  dlon = (point.lon-lon).abs.to_rad

  dphi = Math.log(Math.tan(lat2/2 + Math::PI/4) / Math.tan(lat1/2 + Math::PI/4))

  q = begin
    dlat / dphi
  rescue 
    Math.cos(lat1) # E-W line gives dPhi=0
  end
      
  # if dlon over 180° take shorter rhumb across 180° meridian:
  dlon = 2*Math::PI - dlon if (dlon > Math::PI)

  dist = Math.sqrt(dlat*dlat + q*q*dlon*dlon) * radius; 

  dist.round(4)  # 4 sig figures reflects typical 0.3% accuracy of spherical model
end

#to_lat(format = :dms, dp = 0) ⇒ Object

Returns the latitude of this point; signed numeric degrees if no format, otherwise format & dp as per Geo.to_lat

  • String [format]: Return value as ‘d’, ‘dm’, ‘dms’

  • Numeric [dp=0|2|4]: No of decimal places to display

Returns Numeric|String: Numeric degrees if no format specified, otherwise deg/min/sec



293
294
295
296
 
# File 'lib/geo_calc/calculations.rb', line 293

def to_lat format = :dms, dp = 0
  return lat if !format
  Geo.to_lat lat, format, dp
end

#to_lon(format, dp) ⇒ Object

Returns the longitude of this point; signed numeric degrees if no format, otherwise format & dp as per Geo.toLon()

Parameters:

  • [format]: (String) —

    Return value as ‘d’, ‘dm’, ‘dms’

  • [dp=0|2|4]: (Number) —

    No of decimal places to display



308
309
310
311
 
# File 'lib/geo_calc/calculations.rb', line 308

def to_lon format, dp
  return lon if !format
  Geo.to_lon lon, format, dp
end

#to_s(format = :dms, dp = 0) ⇒ Object 



322
323
324
325
326
327
328
329
330
331
 
# File 'lib/geo_calc/calculations.rb', line 322

def to_s format = :dms, dp = 0 
  format ||= :dms

  return '-,-' if !lat || !lon

  _lat = Geo.to_lat lat, format, dp
  _lon = Geo.to_lon lon, format, dp

  "#{_lat}, #{_lon}"
end