The Museum of HP Calculators

HP Forum Archive 21

 UTM <-> Lat/Lon for WP-34SMessage #1 Posted by Tom Grydeland on 5 Sept 2013, 1:31 p.m. Hi all, Since I saw WGS-84 constants in the catalogue, I thought this handy little machine should be able to do the conversion between Lat/Lon and UTM. After a short detour in Snyder's series, I realised that The Wikipedia entry had a development that was both terser and more accurate (I verified with a small Python implementation that I can transform to UTM and back with an introduced inaccuracy in the micrometer range even when the transformation was forced to use UTM zone 5 away from its "true" zone, once I had gone to n^4 in the expansions.) I have a working implementation on my WP-34s, and I believe the listing below is the same, modulo local label numbers. In a fully loaded calc, I must reduce to 80 global registers to run the code, but I'm sure it can be shortened. The '34 version goes only to order n^3, which means that the inaccuracy is around 0.2 mm instead of micrometers. Good enough for my purposes. Issues: - Does not handle southern hemisphere (users must add/subtract 1e7 meters themselves) - If the first conversion is UTM->Lat/Lon, zone must be prestored [C], but the code does not check for it. - No prompting or instructions Usage: - XEQ 'UTM' // Runs intialization and stops. At this point, the labels A, B and C are active as follows - [A] transforms latitude (Y) and longitude (X), both in decimal degrees, to Northing (Y) and Easting (X) in metres. Rest of stack is preserved. If zone was set manually (with [C]), conversion takes place in the given zone, otherwise zone is given by coordinates. - [B] transforms Northing (Y) and Easting (X) in metres to latitude (Y) and longitude (X), both in decimal degrees. Rest of stack is preserved. The zone must have been set (either with [C] or with a transformation to UTM [A]) before this function can be used, but the program does not check. - [C] is used to recall UTM zone, or to store/force a zone (when ENTRY? is true). Store 0 to return to automatically determined zone. Any number of transformations can be performed while the program is active. Global registers 10-23 should not be touched between transformations, and the remainder of the registers 00-36 will be modified during transformations. Flags 00-03 are also used by the program. Suggestions are welcome, as always. What an amazing little machine to have in one's pocket! Thanks to its creators! --T Edit: had to fix the last bug before reposting the code ```/* * Allocation of registers * r00 - r02 scratch vector T1 * r03 - r05 scratch vector T2 * r06 - r08 scratch vector T3 * r10 - r12 alpha_i * r13 - r15 beta_i * r16 - r18 delta_i * r19 - UTM zone number * r20 - lambda_0, UTM zone center longitude, in radians * r21 - k0 (0.9996) * r22 - A = a/(1+n) * (1 + n**2/4 + n**4/64) * r23 - n (0.00167922039) */ /* * Labels: * 'UTM' : main entry point, initializes fixed vectors (alpha, beta, delta) * and zone if the X value is an integer 1 <= X <= 60 * * A : lat lon -> utm_n utm_e * B : utm_n utm_e -> lat lon * C : sto/rcl UTM zone * * dot2 : two-vector product-sum (dot product) VDy VDx -> X = Y [dot] Z * dot3 : three-vector product-sum VDz VDy VDx -> (X [dot] Y [dot] Z) * v_trig : create vector of trig(i*x) terms; xi Vi -> xi * where 'trig' is one of [cos sin cosh sinh] for flags 1/2 of [0/0 1/0 0/1 1/1] * X points at first register to modify before call, * and returns the angle in X (in case next vector uses same) */ /* * Flags: * 00: set when UTM zone is chosen manually, cleared (default) otherwise * 01: set when v_trig computes sin(h), cleared for cos(h) * 02: set when v_trig uses hyperbolic versions * 03: set when dot3 is used for two vectors only, cleared otherwise * * During init * 01: set when power series in n has order-1 term * 02: set when power series in n has order-2 term */ LBL'UTM' // Initialization: RAD LocR 10 STOS .00 // TODO: Initialize UTM zone // r21: k0 - scale factor at central latitude . 9 9 9 6 STO 21 // r23: n - third flattening # Sf[^-1] 2 * DEC X 1/x STO 23 // r22: A - length of elliptic arc from equator to pole // n is on stack already STO Z INC Z x^2 STO Y # 64 1/x * # 4 1/x + * INC X RCL/ Y # Sa * STO 22 CF 00 // no UTM zone set SF 01 // n-term present SF 02 // n**2-term present CF 03 // alpha_1 # 1/2 # 2 # 3 / +/- STOS 00 // will be used again in beta_1 and delta_1 # 5 # 16 / XEQ npow STO 10 // beta_1 RCLS 00 # 37 # 96 / XEQ npow STO 13 //delta_1 RCLS 00 # 2 STO Z +/- XEQ npow STO 16 CF 01 // n-term absent // alpha_2 # 13 # 48 / # 3 # 5 / +/- XEQ npow STO 11 // beta_2 # 48 1/x # 15 1/x XEQ npow STO 14 // delta_2 # 7 # 3 / # 8 # 5 / +/- XEQ npow STO 17 CF 02 // n**2-term absent // alpha_3 # 61 # 240 / XEQ npow STO 12 // beta_3 # 17 # 48 / # 10 / XEQ npow STO 15 // delta_3 # 56 # 15 / XEQ npow STO 18 RCLS .00 RTN // initialization done, wait for orders /* Convert lat/lon to UTM E/N and zone * * lat/lon to UTM register allocations * r30 - phi (latitude in radians) * r31 - dla (delta lambda in radians) * r32 - t (intermediate value) * r33 - xi' (intermediate angle) * r34 - eta' (intermediate angle) * r.00-.07 saved stack * r.00 - Easting (overwrites longitude) * r.01 - Northing (overwrites latitude) * * On entry: * latitude in degrees in Y * longitude in degrees in X * */ LBL A LocR 10 STOS .00 FC? 00 XEQ zone // At this point, zone is in R19 and lambda0 in r20 DEG[->] // X: longitude in radians RCL- 20 STO 31 x<> Y DEG[->] STO 30 // Compute t SIN ATANH RCL L RCL 23 [sqrt] RCL L INC X / # 2 * <> XYXZ * ATANH * - SINH STO 32 // Compute xi' RCL 31 COS / ATAN STO 33 // Compute cos(2*j*xi') -> T1 [0-2] # 2 * CF 01 CF 02 # 0 XEQ v_trig // Compute sin(2*j*xi') -> T2 [3-5] SF 01 // sin # 3 XEQ v_trig // Compute eta' RCL 31 SIN RCL 32 x^2 INC X [sqrt] / ATANH STO 34 // Compute sinh(2*j*eta') -> T3 [6-8] # 2 * SF 02 # 6 XEQ v_trig // Compute easting # 10 // alpha # 0 // T1 = cos(2*j*xi') # 6 // T3 = sinh(2*j*eta') XEQ dot3 // sum(...) RCL+ 34 RCL* 21 RCL* 22 // A*k0*(eta + sum(...)) 5 EEX 5 + STO .00 // Compute cosh(2*j*eta') -> T1 [0-2] # 2 RCL* 34 CF 01 // cosh # 0 // T1 XEQ v_trig // Compute northing # 10 // alpha # 3 // T2 = sin(2*j*xi') # 0 // T1 = cosh(2*j*eta') XEQ dot3 // sum(...) RCL+ 33 RCL* 21 RCL* 22 // A*k0*(xi + sum(...)) // TODO: Add false northing in southern hemisphere STO .01 // place northing in Y in returned stack // Restore stack and return! RCLS .00 RTN /* Convert UTM E/N and zone to lat/lon * * UTM to lat/lon register allocations * r30 - phi (latitude in radians) * r31 - dla (delta lambda in radians) * r32 - chi (intermediate value) * r33 - xi' (intermediate angle) * r34 - eta' (intermediate angle) * r35 - xi (rescaled latitude) * r36 - eta (rescaled longitude) * r.00-.07 saved stack * * On entry: * northing in metres in Y * easting in metres in X * */ LBL B LocR 10 STOS .00 // FC? 00 // prompt for zone 5 EEX 5 - RCL/ 21 RCL/ 22 STO 36 // eta = (E-E0)/(k0*A) x<> Y // TODO: subtract false northing in southern hemisphere RCL/ 21 RCL/ 22 STO 35 // xi = (N-N0)/(k0*A) // Compute cos(2*j*xi) -> T1 [0-2] # 2 * CF 01 CF 02 # 0 XEQ v_trig // Compute sin(2*j*xi) -> T2 [3-5] SF 01 # 3 XEQ v_trig // Compute sinh(2*j*eta) -> T3 [6-8] # 2 RCL* 36 SF 02 // sinh # 6 XEQ v_trig // compute eta' # 13 // beta # 0 // T1 = cos(2*j*xi) # 6 // T3 = sinh(2*j*eta) XEQ dot3 RCL- 36 +/- STO 34 // eta' = eta - sum(...) // Compute cosh(2*j*eta) -> T1 [0-2] # 2 RCL* 36 CF 01 // cosh # 0 // T1 XEQ v_trig // Compute xi' # 13 // beta # 3 // T2 = sin(2*j*xi) # 0 // T1 = cosh(2*j*eta) XEQ dot3 RCL- 35 +/- STO 33 // xi' = xi - sum(...) // Compute chi SIN RCL 34 COSH / ASIN STO 32 // Compute phi, which is latitude in radians # 2 * SF 01 CF 02 // sin # 0 XEQ v_trig # 16 // delta # 0 // T1 = sin(2*j*chi) XEQ dot2 RCL+ 32 // phi = chi + sum(...) rad[->][degree] // latitude STO .01 // delta lambda, then longitude RCL 34 SINH RCL 33 COS / ATAN // delta lambda RCL+ 20 // Translate to correct UTM zone rad[->][degree] // longitude STO .00 RCLS .00 RTN LBL C ENTRY? // Setting zone JMP set_zone RCL 19 RTN set_zone:: STO 19 x[!=]0? SKIP 3 STO 20 CF 00 RTN # 6 * # 183 - DEG[->] STO 20 DROP RCL 19 SF 00 RTN // Determine correct UTM zone // Y: lat (in degrees) // X: lon (in degrees) // Stack is left undisturbed, zone in r19, lambda_0 in r20 zone:: LocR 10 // STOS bug, could be 8 (or 4) STOS .00 // lon in .00, lat in .01 # 6 / FLOOR # 31 + STO 19 RCL L // look for SW Norway anomaly x[!=]? 19 JMP zone_next # 56 x>? .01 // lat < 56? JMP zone_next # 64 x 64? JMP zone_next INC 19 JMP zone_done zone_next:: # 72 // look for Svalbard anomaly x>? .01 // lat < 72? JMP zone_done # 6 x>? .00 // lon < 6? JMP zone_done # 36 x 36? JMP zone_done RCL .00 # 3 + # 12 / FLOOR # 2 * # 31 + STO 19 zone_done:: RCL 19 # 6 * # 183 - DEG[->] STO 20 RCLS .00 RTN // Z: label to execute for xi -> trig(xi) // Y: xi // X: first register (of 3) v_trig:: STO I // start of vector in I, theta in J x<> Y STO J XEQ trig STO[->]I INC I # 2 RCL* J XEQ trig STO[->]I INC I # 3 RCL* J XEQ trig STO[->]I RCL J RTN trig:: FS? 01 JMP trig_sin FS? 02 COSH FC? 02 COS RTN trig_sin:: FS? 02 SINH FC? 02 SIN RTN // takes three power series coefficients from the stack, n^3 in x, n^2 in y n in z // used to compute alpha_i, beta_i, delta_i npow:: RCL* 23 FS? 02 + RCL* 23 FS? 01 + RCL* 23 RTN // Takes vector start numbers in X, and Y, drops them both // off the stack and returns their dot product in X. // Uses only local registers. dot2:: SF 03 // Takes vector start numbers in X, Y, and Z, drops them all // off the stack and returns their triple-dot product in X. // Uses only local registers. dot3:: LocR 4 STO .00 DROP STO .01 DROP FS? 03 SKIP 2 STO .02 DROP # 3 STO .03 DROP # 0 dot3_top:: RCL[->].00 RCL*[->].01 FC? 03 RCL*[->].02 + INC .00 INC .01 FC? 03 INC .02 DSE .03 JMP dot3_top CF 03 END ``` Edited: 5 Sept 2013, 5:27 p.m.

 Re: UTM <-> Lat/Lon test vectorsMessage #2 Posted by Tom Grydeland on 5 Sept 2013, 5:41 p.m.,in response to message #1 by Tom Grydeland Edit: mispasted Surtsey results. Must have been tired yesterday. These values computed with the series extended to n^4 A point on the equator, at the centre of a zone: `(0, 3) <-> (5e5, 0, 31)` Surtsey, Iceland; native zone and forced to a different zone: ```(63.303, -20.6047) <-> (519815.0488289792, 7019410.3838700978, 27) (63.303, -20.6047) <-> (-1211867.9773856564, 7516864.8519201223, 33)``` The Duma in Moscow, computing at the coordinates given by WP 1) Lat/Lon given: `(55.758056, 37.615278) <-> (413102.2660859063, 6180020.5373971062, 37)` 2) UTM given: `(55.758051124985862, 37.615273932357319) <-> (413102, 6180020, 37)` Edited: 6 Sept 2013, 3:22 a.m.

Go back to the main exhibit hall