class Silva::Transform
Encapsulates the hairy maths required to perform the various transforms. Checked against documentation provided by Orndance Survey: www.ordnancesurvey.co.uk/oswebsite/gps/osnetfreeservices/furtherinfo/questdeveloper.html www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/guide6.html
Constants
- AIRY1830
The Airy ellipsoid used by the OSGB36 datum.
- DEGREE_ROUNDING_PLACES
Decimal places to round results in degrees.
- E0
The easting of true origin
- F0
The scale factor on central meridian
- GRS80
The
GRS80
ellipsoid used by the WGS84 datum.- HELMERT_PARAMS
Helmert transform parameters
- LAMBDA0
The longitude of true origin, in radians
- N0
The northing of true origin
- PHI0
The latitude of true origin, in radians
Public Class Methods
Convert an OSGB36 easting, northing system to :osgb36 co-ordinate system
@param [Silva::System::En] osgb36 A valid :en location system. @return [Silva::System::Osgb36] A valid :osgb36 location system.
# File lib/silva/transform.rb, line 73 def self.en_to_osgb36(en) # Algorithm from: # http://www.ordnancesurvey.co.uk/oswebsite/gps/osnetfreeservices/furtherinfo/questdeveloper.html a, b = AIRY1830[:a], AIRY1830[:b] e = en.easting n = en.northing phi_prime = PHI0 m = 0 while n - N0 - m >= 0.0001 do phi_prime = (n - N0 - m) / (a * F0) + phi_prime m = meridional_arc(phi_prime) end nu, rho, eta2 = transverse_and_meridional_radii(phi_prime) vii = Math.tan(phi_prime) / (2 * rho * nu) viii = Math.tan(phi_prime)/(24 * rho * nu**3) * \ (5 + 3 * Math.tan(phi_prime)**2 + eta2 - 9 * Math.tan(phi_prime)**2 * eta2) ix = Math.tan(phi_prime) / (720 * rho * nu**5) * (61 + 90 * Math.tan(phi_prime)**2 + 45 * Math.tan(phi_prime)**4) x = (1 / Math.cos(phi_prime)) / nu xi = (1 / Math.cos(phi_prime)) / (6 * nu**3) * (nu / rho + 2 * Math.tan(phi_prime)**2) xii = (1 / Math.cos(phi_prime)) / (120 * nu**5) * (5 + 28 * Math.tan(phi_prime)**2 + 24 * Math.tan(phi_prime)**4) xiia = (1 / Math.cos(phi_prime)) / (5040 * nu**7) * \ (61 + 662 * Math.tan(phi_prime)**2 + 1320 * Math.tan(phi_prime)**4 + 720 * Math.tan(phi_prime)**6) phi = phi_prime - vii * (e - E0)**2 + viii * (e - E0)**4 - ix * (e - E0)**6 lambda = LAMBDA0 + x * (e - E0) - xi * (e - E0)**3 + xii * (e - E0)**5 - xiia * (e - E0)**7 System.create(:osgb36, :lat => to_deg(phi).round(DEGREE_ROUNDING_PLACES), \ :long => to_deg(lambda).round(DEGREE_ROUNDING_PLACES), :alt => 0) end
Convert an :osgb36 co-ordinate system to :en (eastings and northings)
@param [Silva::System::Osgb36] osgb36 A valid :wgs84 location system. @return [Silva::System::En] A valid :en location system.
# File lib/silva/transform.rb, line 63 def self.osgb36_to_en(osgb36) self.latlong_to_en(osgb36, AIRY1830) end
Convert an :osgb36 co-ordinate system :wgs84
@param [Silva::System::Osgb36] osgb36 A valid :osgb36 location system. @return [Silva::System::Wgs84] A valid :wgs84 location system.
# File lib/silva/transform.rb, line 43 def self.osgb36_to_wgs84(osgb36) helmert_transform(osgb36, :wgs84, AIRY1830, HELMERT_PARAMS.inject({}) { |h, (k, v)| h[k] = v * -1; h }, GRS80) end
Convert a :wgs84 co-ordinate system to :osgb36
@param [Silva::System::Wgs84] wgs84 A valid :wgs84 location system. @return [Silva::System::Osgb36] A valid :osgb36 location system.
# File lib/silva/transform.rb, line 53 def self.wgs84_to_osgb36(wgs84) helmert_transform(wgs84, :osgb36, GRS80, HELMERT_PARAMS, AIRY1830) end
Private Class Methods
Calculate eccentricity**2 given relevant ellipsoid.
# File lib/silva/transform.rb, line 220 def self.eccentricity_squared(ellipsoid) (ellipsoid[:a]**2 - ellipsoid[:b]**2) / ellipsoid[:a]**2 end
Transform
to/from :osgb36/:wgs84 co-ordinate systems. Algorithm from: www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/guide6.html Portions of code from: www.harrywood.co.uk/blog/2010/06/29/ruby-code-for-converting-to-uk-ordnance-survey-coordinate-systems-from-wgs84/
# File lib/silva/transform.rb, line 138 def self.helmert_transform(source_system, target_system, ellipsoid_1, transform, ellipsoid_2) phi = to_rad(source_system.lat) lambda = to_rad(source_system.long) h = source_system.alt a1 = ellipsoid_1[:a] # convert co-ordinates to 3D Cartesian. See: # http://www.ordnancesurvey.co.uk/oswebsite/gps/docs/A_Guide_to_Coordinate_Systems_in_Great_Britain.pdf e_sq1 = eccentricity_squared(ellipsoid_1) nu = a1 / Math.sqrt(1 - e_sq1 * Math.sin(phi)**2) x1 = (nu + h) * Math.cos(phi) * Math.cos(lambda) y1 = (nu + h) * Math.cos(phi) * Math.sin(lambda) z1 = ((1 - e_sq1) * nu + h) * Math.sin(phi) # apply Helmert transformation tx = transform[:tx] ty = transform[:ty] tz = transform[:tz] rx = transform[:rx] / 3600 * Math::PI / 180 ry = transform[:ry] / 3600 * Math::PI / 180 rz = transform[:rz] / 3600 * Math::PI / 180 s1 = transform[:s] / 1e6 + 1 x2 = tx + x1 * s1 - y1 * rz + z1 * ry y2 = ty + x1 * rz + y1 * s1 - z1 * rx z2 = tz - x1 * ry + y1 * rx + z1 * s1 # convert 3D Cartesian co-ordinates back to lat, long, alt a2 = ellipsoid_2[:a] precision = 4 / a2 e_sq2 = eccentricity_squared(ellipsoid_2) p = Math.sqrt(x2**2 + y2**2) phi = Math.atan2(z2, p * (1 - e_sq2)) phi_prime = 2 * Math::PI while ((phi - phi_prime).abs > precision) do nu = a2 / Math.sqrt(1 - e_sq2 * Math.sin(phi)**2) phi_prime = phi phi = Math.atan2(z2 + e_sq2 * nu * Math.sin(phi), p) end lambda = Math.atan2(y2, x2) h = p / Math.cos(phi) - nu System.create(target_system, :lat => to_deg(phi).round(DEGREE_ROUNDING_PLACES),\ :long => to_deg(lambda).round(DEGREE_ROUNDING_PLACES), :alt => h) end
Convert an :osgb36 or :wgs84 co-ordinate system to OSGB36 easting and northing.
# File lib/silva/transform.rb, line 110 def self.latlong_to_en(system, ellipsoid) phi = to_rad(system.lat) lambda = to_rad(system.long) nu, rho, eta2 = transverse_and_meridional_radii(phi, ellipsoid) m = meridional_arc(phi, ellipsoid) i = m + N0 ii = (nu / 2) * Math.sin(phi) * Math.cos(phi) iii = (nu / 24) * Math.sin(phi) * Math.cos(phi)**3 * (5 - Math.tan(phi)**2 + 9 * eta2) iiia = (nu / 720) * Math.sin(phi) * Math.cos(phi)**5 * (61 - 58 * Math.tan(phi)**2 + Math.tan(phi)**4) iv = nu * Math.cos(phi) v = (nu / 6) * Math.cos(phi)**3 * (nu / rho - Math.tan(phi)**2) vi = (nu / 120) * Math.cos(phi)**5 * \ (5 - 18 * Math.tan(phi)**2 + Math.tan(phi)**4 + 14 * eta2 - 58 * Math.tan(phi)**4 * eta2) n = i + ii * (lambda - LAMBDA0)**2 + iii * (lambda - LAMBDA0)**4 + iiia * (lambda - LAMBDA0)**6 e = E0 + iv * (lambda - LAMBDA0) + v * (lambda - LAMBDA0)**3 + vi * (lambda - LAMBDA0)**5 System.create(:en, :easting => e.round, :northing => n.round) end
Calculate M (meridional arc) given latitude and relevant ellipsoid
# File lib/silva/transform.rb, line 195 def self.meridional_arc(phi, ellipsoid = AIRY1830) a, b = ellipsoid[:a], ellipsoid[:b] n = n(ellipsoid) ma = (1 + n + (5.0 / 4.0) * n**2 + (5.0 / 4.0) * n**3) * (phi - PHI0) mb = (3 * n + 3 * n**2 + (21.0 / 8.0) * n**3) * Math.sin(phi - PHI0) * Math.cos(phi + PHI0) mc = ((15.0 / 8.0) * n**2 + (15.0 / 8.0) * n**3) * Math.sin(2 * (phi - PHI0)) * Math.cos(2 * (phi + PHI0)) md = (35.0 / 24.0) * n**3 * Math.sin(3 * (phi - PHI0)) * Math.cos(3 * (phi + PHI0)) b * F0 * (ma - mb + mc - md) end
Calculate required n parameter given the relevant ellipsoid
# File lib/silva/transform.rb, line 190 def self.n(ellipsoid) (ellipsoid[:a] - ellipsoid[:b]) / (ellipsoid[:a] + ellipsoid[:b]) end
Radians to degrees.
# File lib/silva/transform.rb, line 230 def self.to_deg(rads) rads * 180 / Math::PI end
Degrees to radians.
# File lib/silva/transform.rb, line 225 def self.to_rad(degrees) degrees * Math::PI / 180 end
Calculate nu, rho, eta2 (transverse and meridional radii) given latitude and relevant ellipsoid.
# File lib/silva/transform.rb, line 208 def self.transverse_and_meridional_radii(phi, ellipsoid = AIRY1830) a, b = ellipsoid[:a], ellipsoid[:b] e2 = eccentricity_squared(ellipsoid) nu = a * F0 / Math.sqrt(1 - e2 * Math.sin(phi)**2) rho = a * F0 * (1 - e2) / ((1 - e2 * Math.sin(phi)**2)**1.5) eta2 = nu / rho - 1 [nu, rho, eta2] end