The remit here is to take Eastings/Northings (A UK based Ordnance Survey Gridreferencing system) and convert to LatLngs (the grid-referencing done by most GPS systems, and most importantly, Google Maps).
The algorithms for doing these conversions are here, and the accompanying guide I found really interesting. Basically, it comes down to Eastings/Northings are based on a 2-dimensional grid (assuming that England is roughly flat), but LatLngs are based on a 3-dimensional sphere (I'm looking at you, planet Earth...)
This call for me had the perfect storm of lots of complicated maths, knowledge of 2 programming languages, and the sexy prospect of overlaying the results on Google Maps.
Total trophy moment when I got it to work....
*&---------------------------------------------------------------------*
*& Report YPD_PL1_LATLNG
*&
*&---------------------------------------------------------------------*
*&
*&
*&---------------------------------------------------------------------*
REPORT YPD_PL1_LATLNG.
data: pi type f,
radPHI0 type f,
radLAM0 type f,
af0 type f,
bf0 type f,
e2 type f,
Et type f,
phi0 type f,
lam0 type f,
w_deciPhi0 type f,
w_DeciLam0 type f,
o_lat type f,
o_lng type f.
constants: c_a type f value '6377563.396',
c_b type f value '6356256.909',
c_e0 type f value '400000',
c_n0 type f value '100000-',
c_f0 type f value '0.9996012717',
c_deciPhi0 type f value '49.00000000',
c_deciLam0 type f value '2.00000000-',
c_Phi0 type f VALUE '0.855211333',
C_Lam0 type f value '0.034906585-'.
parameters: p_east type int4 default '359536',
p_north type int4 default '403646'.
start-of-selection.
write: 'East:', p_east, ', North', p_north.
perform initialise_vars.
perform E_N_to_LAT using p_east p_north c_n0 c_f0 changing o_lat.
perform E_N_to_LNG using p_east p_north c_n0 c_f0 changing o_lng.
form E_N_to_Lat using p_East type int4
p_North type int4
p_n0 type f
p_f0 type f
changing e_n_to_lat type f.
data: af0 type f,
bf0 type f,
e2 type f,
n type f,
Et type f,
VII type f,
VIII type f,
IX type f,
PHid type f,
nu type f,
rho type f,
eta2 type f.
*'Un-project Transverse Mercator eastings and northings back to latitude.
*'Input: - _
* eastings (East) and northings (North) in meters; _
* ellipsoid axis dimensions (a & b) in meters; _
* eastings (e0) and northings (n0) of false origin in meters; _
* central meridian scale factor (f0) and _
* latitude (PHI0) and longitude (LAM0) of false origin in decimal degrees.
"'REQUIRES THE "Marc" AND "InitialLat" FUNCTIONS
*'Compute af0, bf0, e squared (e2), n and Et
af0 = c_a * p_f0.
bf0 = c_b * p_f0.
e2 = ( ( af0 ** 2 ) - ( bf0 ** 2 ) ) / ( af0 ** 2 ) .
n = ( af0 - bf0 ) / ( af0 + bf0 ).
Et = p_East - c_e0.
*'Compute initial value for latitude (PHI) in radians
perform InitialLat using p_north c_n0 af0 RadPHI0 n bf0 changing PHId.
*'Compute nu, rho and eta2 using value for PHId
nu = af0 / ( Sqrt( 1 - ( e2 * ( ( Sin( PHId ) ) ** 2 ) ) ) ).
rho = ( nu * ( 1 - e2 ) ) / ( 1 - ( e2 * ( SIN( PHId ) ) ** 2 ) ).
eta2 = ( nu / rho ) - 1.
*'Compute Latitude
VII = ( Tan( PHId ) ) / ( 2 * rho * nu ).
VIII = ( ( Tan( PHId ) ) / ( 24 * rho * ( nu ** 3 ) ) ) * ( 5 + ( 3 * ( ( Tan( PHId ) ) ** 2 ) ) + eta2 - ( 9 * eta2 * ( ( Tan( PHId ) ) ** 2 ) ) ).
IX = ( ( Tan( PHId ) ) / ( 720 * rho * ( nu ** 5 ) ) ) * ( 61 + ( 90 * ( ( Tan( PHId ) ) ** 2 ) ) + ( 45 * ( ( Tan( PHId ) ) ** 4 ) ) ).
E_N_to_Lat = ( 180 / Pi ) * ( PHId - ( ( Et ** 2 ) * VII ) + ( ( Et ** 4 ) * VIII ) - ( ( Et ** 6 ) * IX ) ).
WRITE:/ 'Lat Calculation'.
write:/ 'Lat:', e_n_to_lat.
write:/ 'VII', VII.
write:/ 'IX', IX.
endform.
form initialise_vars.
Pi = '3.14159265358979' .
"Convert angle measures to radians
RadPHI0 = c_DeciPhi0 * ( Pi / 180 ).
RadLAM0 = c_DeciLam0 * ( Pi / 180 ).
endform.
form marc using p_bf0 type f
p_n type f
p_phi0 type f
p_phi type f
changing p_marc type f.
data: phi_minus_phi0 type f,
phi_plus_phi0 type f,
phi_minus_phi0_times_2 type f,
phi_plus_phi0_times_2 type f,
phi_minus_phi0_times_3 type f,
phi_plus_phi0_times_3 type f.
data: l_test1 type f,
l_line1 type f,
l_line2 type f,
l_line3 type f,
l_line4 type f.
DATA: L_MARC2 TYPE F.
"Copied from VB Macro Calculation; They can nest function calculations, ABAP can't... stupid ABAP.
phi_minus_phi0 = p_phi - p_phi0.
phi_plus_phi0 = p_phi + p_phi0.
phi_minus_phi0_times_2 = phi_minus_phi0 * 2.
phi_plus_phi0_times_2 = phi_plus_phi0 * 2.
phi_minus_phi0_times_3 = phi_minus_phi0 * 3.
phi_plus_phi0_times_3 = phi_plus_phi0 * 3.
l_test1 = Sin( Phi_minus_phi0 ).
l_line1 = p_bf0.
l_line1 = ( ( ( 1 + p_n + ( ( 5 / 4 ) * ( p_n ** 2 ) ) + ( ( 5 / 4 ) * ( p_n ** 3 ) ) ) * ( p_PHI - p_PHI0 ) ) ).
l_line2 = ( ( ( 3 * p_n ) + ( 3 * ( p_n ** 2 ) ) + ( ( 21 / 8 ) * ( p_n ** 3 ) ) ) * ( Sin( phi_minus_phi0 ) ) * ( COS( phi_plus_phi0 ) ) ).
l_line3 = ( ( ( ( 15 / 8 ) * ( p_n ** 2 ) ) + ( ( 15 / 8 ) * ( p_n ** 3 ) ) ) * ( Sin( phi_minus_phi0_times_2 ) ) * ( Cos( phi_plus_phi0_times_2 ) ) ).
l_line4 = ( ( ( 35 / 24 ) * ( p_n ** 3 ) ) * ( Sin( phi_minus_phi0_times_3 ) ) * ( Cos( phi_plus_phi0_times_3 ) ) ).
* Marc = bf0 * (((1 + n + ((5 / 4) * (n ^ 2)) + ((5 / 4) * (n ^ 3))) * (PHI - PHI0)) _
* - (((3 * n) + (3 * (n ^ 2)) + ((21 / 8) * (n ^ 3))) * (Sin(PHI - PHI0)) * (Cos(PHI + PHI0))) _
* + ((((15 / 8) * (n ^ 2)) + ((15 / 8) * (n ^ 3))) * (Sin(2 * (PHI - PHI0))) * (Cos(2 * (PHI + PHI0)))) _
* - (((35 / 24) * (n ^ 3)) * (Sin(3 * (PHI - PHI0))) * (Cos(3 * (PHI + PHI0)))))
L_MARC2 = p_bf0 * ( l_line1 - l_line2 + l_line3 - l_line4 ).
p_marc = l_marc2.
* p_Marc = p_bf0 * ( ( ( 1 + p_n + ( ( 5 / 4 ) * ( p_n ** 2 ) ) + ( ( 5 / 4 ) * ( p_n ** 3 ) ) ) * ( p_PHI - PHI0 ) )
* - ( ( ( 3 * p_n ) + ( 3 * ( p_n ** 2 ) ) + ( ( 21 / 8 ) * ( p_n ** 3 ) ) ) * ( Sin( phi_minus_phi0 ) ) * ( COS( phi_plus_phi0 ) ) )
* + ( ( ( ( 15 / 8 ) * ( p_n ** 2 ) ) + ( ( 15 / 8 ) * ( p_n ** 3 ) ) ) * ( Sin( phi_minus_phi0_times_2 ) ) * ( Cos( phi_plus_phi0_times_2 ) ) )
* - ( ( ( 35 / 24 ) * ( p_n ** 3 ) ) * ( Sin( phi_minus_phi0_times_3 ) ) * ( Cos( phi_plus_phi0_times_3 ) ) ) ) .
endform.
form InitialLat using p_North type int4
p_n0 type f
P_afo type f
p_PHI0 type f
p_n type f
p_bfo type f
changing p_InitialLat type f.
*'Compute initial value for Latitude (PHI) IN RADIANS.
*'Input: - _
* northing of point (North) and northing of false origin (n0) in meters; _
* semi major axis multiplied by central meridian scale factor (af0) in meters; _
* latitude of false origin (PHI0) IN RADIANS; _
* n (computed from a, b and f0) and _
* ellipsoid semi major axis multiplied by central meridian scale factor (bf0) in meters.
*'REQUIRES THE "Marc" FUNCTION
*'THIS FUNCTION IS CALLED BY THE "E_N_to_Lat", "E_N_to_Long" and "E_N_to_C" FUNCTIONS
*'THIS FUNCTION IS ALSO USED ON IT'S OWN IN THE "Projection and Transformation Calculations.xls" SPREADSHEET
data: l_phi1 type f,
l_phi2 type f,
M type f,
delta type f,
l_myCounter type f.
"First PHI value (PHI1)
l_PHI1 = ( ( p_North - p_n0 ) / p_afo ) + P_PHI0.
*'Calculate M
perform marc using p_bfo
p_n
p_phi0
l_phi1
changing M.
*'Calculate new PHI value (PHI2)
l_PHI2 = ( ( p_North - p_n0 - M ) / p_afo ) + l_PHI1.
*'Iterate to get final value for InitialLat
Do 10000 times. "Prevent Infinity.
delta = p_North - p_n0 - M.
if ABS( delta ) < '0.00001'.
exit.
endif.
l_PHI2 = ( DELTA / P_afo ) + l_PHI1.
perform marc using p_bfo
p_n
p_phi0
l_phi2
changing M.
l_PHI1 = l_PHI2.
l_myCounter = l_myCounter + 1.
enddo.
p_InitialLat = l_PHI2.
endform.
form E_N_to_LNG using p_East type int4
p_North type int4
p_n0 type f
p_f0 type f
changing e_n_to_lng type f.
"Un-project Transverse Mercator eastings and northings back to longitude.
data: af0 type f,
bf0 type f,
e2 type f,
n type f,
Et type f,
VII type f,
VIII type f,
IX type f,
PHid type f,
nu type f,
rho type f,
eta2 type f,
X type f,
XI type f,
XII type f,
XIIA type f.
"Compute af0, bf0, e squared (e2), n and Et
af0 = c_a * p_f0.
bf0 = c_b * p_f0.
e2 = ( ( af0 ** 2 ) - ( bf0 ** 2 ) ) / ( af0 ** 2 ).
n = ( af0 - bf0 ) / ( af0 + bf0 ).
Et = p_East - c_e0.
*'Compute initial value for latitude (PHI) in radians
perform InitialLat using p_north c_n0 af0 RadPHI0 n bf0 changing PHId.
*'Compute nu, rho and eta2 using value for PHId
nu = af0 / ( Sqrt( 1 - ( e2 * ( ( Sin( PHId ) ) ** 2 ) ) ) ).
rho = ( nu * ( 1 - e2 ) ) / ( 1 - ( e2 * ( SIN( PHId ) ) ** 2 ) ).
eta2 = ( nu / rho ) - 1.
" Compute Longitude
X = ( ( Cos( PHId ) ) ** -1 ) / nu.
XI = ( ( ( Cos( PHId ) ) ** -1 ) / ( 6 * ( nu ** 3 ) ) ) * ( ( nu / rho ) + ( 2 * ( ( Tan( PHId ) ) ** 2 ) ) ).
XII = ( ( ( Cos( PHId ) ) ** -1 ) / ( 120 * ( nu ** 5 ) ) ) * ( 5 + ( 28 * ( ( Tan( PHId ) ) ** 2 ) ) + ( 24 * ( ( Tan( PHId ) ) ** 4 ) ) ).
XIIA = ( ( ( Cos( PHId ) ) ** -1 ) / ( 5040 * ( nu ** 7 ) ) ) * ( 61 + ( 662 * ( ( Tan( PHId ) ) ** 2 ) ) + ( 1320 * ( ( Tan( PHId ) ) ** 4 ) ) + ( 720 * ( ( Tan( PHId ) ) ** 6 ) ) ).
E_N_to_Lng = ( 180 / Pi ) * ( RadLAM0 + ( Et * X ) - ( ( Et ** 3 ) * XI ) + ( ( Et ** 5 ) * XII ) - ( ( Et ** 7 ) * XIIA ) ).
WRITE:/ 'Lng Calculation'.
write:/ 'Lat:', e_n_to_lng.
Endform.
No comments:
Post a Comment