Created
June 23, 2012 00:52
-
-
Save andflett/2976046 to your computer and use it in GitHub Desktop.
Convert OSGB36 Eastings/Northings Coords to WGS84 Latitude/Longitude
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
module Cartography | |
extend self | |
# Takes OSGB36 Easting/Northing coords | |
# and returns WGS84 Latitude and Longitude | |
def en_to_ll(easting, northing) | |
@OSGB36_ll = to_OSGB36(easting, northing) | |
to_WGS84(@OSGB36_ll[:latitude], @OSGB36_ll[:longitude]) | |
end | |
# Takes OSGB36 Easting/Northing coords and returns | |
# OSGB36 Latitude and Longitude | |
def to_OSGB36(easting, northing) | |
@OSGB_F0 = 0.9996012717 | |
@N0 = -100000.0 | |
@E0 = 400000.0 | |
@phi0 = deg_to_rad(49.0) | |
@lambda0 = deg_to_rad(-2.0) | |
@a = 6377563.396 | |
@b = 6356256.909 | |
@eSquared = ((@a * @a) - (@b * @b)) / (@a * @a) | |
@phi = 0.0 | |
@lambda = 0.0 | |
@E = easting | |
@N = northing | |
@n = (@a - @b) / (@a + @b) | |
@M = 0.0 | |
@phiPrime = ((@N - @N0) / (@a * @OSGB_F0)) + @phi0 | |
begin | |
@M = | |
(@b * @OSGB_F0)\ | |
* (((1 + @n + ((5.0 / 4.0) * @n * @n) + ((5.0 / 4.0) * @n * @n * @n))\ | |
* (@phiPrime - @phi0))\ | |
- (((3 * @n) + (3 * @n * @n) + ((21.0 / 8.0) * @n * @n * @n))\ | |
* Math.sin(@phiPrime - @phi0)\ | |
* Math.cos(@phiPrime + @phi0))\ | |
+ ((((15.0 / 8.0) * @n * @n) + ((15.0 / 8.0) * @n * @n * @n))\ | |
* Math.sin(2.0 * (@phiPrime - @phi0))\ | |
* Math.cos(2.0 * (@phiPrime + @phi0)))\ | |
- (((35.0 / 24.0) * @n * @n * @n)\ | |
* Math.sin(3.0 * (@phiPrime - @phi0))\ | |
* Math.cos(3.0 * (@phiPrime + @phi0)))) | |
@phiPrime += (@N - @N0 - @M) / (@a * @OSGB_F0) | |
end while ((@N - @N0 - @M) >= 0.001) | |
@v = @a * @OSGB_F0 * ((1.0 - @eSquared * sin_pow_2(@phiPrime)) ** -0.5) | |
@rho = | |
@a\ | |
* @OSGB_F0\ | |
* (1.0 - @eSquared)\ | |
* ((1.0 - @eSquared * sin_pow_2(@phiPrime)) ** -1.5) | |
@etaSquared = (@v / @rho) - 1.0 | |
@VII = Math.tan(@phiPrime) / (2 * @rho * @v) | |
@VIII = | |
(Math.tan(@phiPrime) / (24.0 * @rho * (@v ** 3.0)))\ | |
* (5.0\ | |
+ (3.0 * tan_pow_2(@phiPrime))\ | |
+ @etaSquared\ | |
- (9.0 * tan_pow_2(@phiPrime) * @etaSquared)) | |
@IX = | |
(Math.tan(@phiPrime) / (720.0 * @rho * (@v ** 5.0)))\ | |
* (61.0\ | |
+ (90.0 * tan_pow_2(@phiPrime))\ | |
+ (45.0 * tan_pow_2(@phiPrime) * tan_pow_2(@phiPrime))) | |
@X = sec(@phiPrime) / @v | |
@XI = | |
(sec(@phiPrime) / (6.0 * @v * @v * @v))\ | |
* ((@v / @rho) + (2 * tan_pow_2(@phiPrime))) | |
@XII = | |
(sec(@phiPrime) / (120.0 * (@v ** 5.0)))\ | |
* (5.0\ | |
+ (28.0 * tan_pow_2(@phiPrime))\ | |
+ (24.0 * tan_pow_2(@phiPrime) * tan_pow_2(@phiPrime))) | |
@XIIA = | |
(sec(@phiPrime) / (5040.0 * (@v ** 7.0)))\ | |
* (61.0\ | |
+ (662.0 * tan_pow_2(@phiPrime))\ | |
+ (1320.0 * tan_pow_2(@phiPrime) * tan_pow_2(@phiPrime))\ | |
+ (720.0\ | |
* tan_pow_2(@phiPrime)\ | |
* tan_pow_2(@phiPrime)\ | |
* tan_pow_2(@phiPrime))) | |
@phi = | |
@phiPrime\ | |
- (@VII * ((@E - @E0) ** 2.0))\ | |
+ (@VIII * ((@E - @E0) ** 4.0))\ | |
- (@IX * ((@E - @E0) ** 6.0)) | |
@lambda = | |
@lambda0\ | |
+ (@X * (@E - @E0))\ | |
- (@XI * ((@E - @E0) ** 3.0))\ | |
+ (@XII * ((@E - @E0) ** 5.0))\ | |
- (@XIIA * ((@E - @E0) ** 7.0)) | |
{ :latitude => rad_to_deg(@phi), :longitude => rad_to_deg(@lambda) } | |
end | |
# Takes OSGB36 Latitude and Longitude coords | |
# and returns WGS84 Latitude and Longitude | |
def to_WGS84(latitude,longitude) | |
@a = 6377563.396 | |
@b = 6356256.909 | |
@eSquared = ((@a * @a) - (@b * @b)) / (@a * @a) | |
@phi = deg_to_rad(latitude) | |
@lambda = deg_to_rad(longitude) | |
@v = @a / (Math.sqrt(1 - @eSquared * sin_pow_2(@phi))) | |
@H = 0 | |
@x = (@v + @H) * Math.cos(@phi) * Math.cos(@lambda) | |
@y = (@v + @H) * Math.cos(@phi) * Math.sin(@lambda) | |
@z = ((1 - @eSquared) * @v + @H) * Math.sin(@phi) | |
@tx = 446.448 | |
@ty = -124.157 | |
@tz = 542.060 | |
@s = -0.0000204894 | |
@rx = deg_to_rad( 0.00004172222) | |
@ry = deg_to_rad( 0.00006861111) | |
@rz = deg_to_rad( 0.00023391666) | |
@xB = @tx + (@x * (1 + @s)) + (-@rx * @y) + (@ry * @z) | |
@yB = @ty + (@rz * @x) + (@y * (1 + @s)) + (-@rx * @z) | |
@zB = @tz + (-@ry * @x) + (@rx * @y) + (@z * (1 + @s)) | |
@a = 6378137.000 | |
@b = 6356752.3141 | |
@eSquared = ((@a * @a) - (@b * @b)) / (@a * @a) | |
@lambdaB = rad_to_deg(Math.atan(@yB / @xB)) | |
@p = Math.sqrt((@xB * @xB) + (@yB * @yB)) | |
@phiN = Math.atan(@zB / (@p * (1 - @eSquared))) | |
(1..10).each do |i| | |
@v = @a / (Math.sqrt(1 - @eSquared * sin_pow_2(@phiN))) | |
@phiN1 = Math.atan((@zB + (@eSquared * @v * Math.sin(@phiN))) / @p) | |
@phiN = @phiN1 | |
end | |
@phiB = rad_to_deg(@phiN) | |
{ :latitude => @phiB, :longitude => @lambdaB } | |
end | |
private # Some Math | |
def deg_to_rad(degrees) | |
degrees / 180.0 * Math::PI | |
end | |
def rad_to_deg(r) | |
(r/Math::PI)*180 | |
end | |
def sin_pow_2(x) | |
Math.sin(x) * Math.sin(x) | |
end | |
def cos_pow_2(x) | |
Math.cos(x) * Math.cos(x) | |
end | |
def tan_pow_2(x) | |
Math.tan(x) * Math.tan(x) | |
end | |
def sec(x) | |
1.0 / Math.cos(x) | |
end | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment