/* --------------------------------------------------------------- */ /* Lat-Lon to UTM converter */ /* --------------------------------------------------------------- */ /* Copyright (c) Mike Cowlishaw, 2007-2012. All rights reserved. */ /* Parts Copyright (c) IBM, 2007. */ /* */ /* Permission to use, copy, modify, and distribute this software */ /* for any non-commercial purpose without fee is hereby granted, */ /* provided that the above copyright notice and this permission */ /* notice appear in all copies, and that notice and the date of */ /* any modifications be added to the software. */ /* */ /* This software is provided "as is". No warranties, whether */ /* express, implied, or statutory, including, but not limited to, */ /* implied warranties of merchantability and fitness for a */ /* particular purpose apply to this software. The author shall */ /* not, in any circumstances, be liable for special, incidental, */ /* or consequential damages, for any reason whatsoever. */ /* */ /* --------------------------------------------------------------- */ /* */ /* Two words are required, by default in decimal degrees: */ /* */ /* lat -- latitude, north positive [-80 -> +84] */ /* lon -- longitude, east positive [-180 -> +180] */ /* */ /* or these may be specified in MemoWiki format (Lat starts with N */ /* or S, Lon starts with E or W, for example: N43:18:12 and */ /* W4:34:13); if these forms are used they are assumed to be valid */ /* (pre-checked) and have at least degrees, with optional minutes */ /* and seconds, separated by ':'s. */ /* */ /* A third word may be specified giving the datum (and hence the */ /* ellipsoid). This may be one of: EUR50, EUR79, or WGS84. If */ /* not given, WGS84 is assumed. */ /* */ /* If a parameter is invalid or out of range, a string starting */ /* with '*' followed by text describing the error is returned. */ /* */ /* Otherwise returns UTM zone, X, and Y as three words; X and Y */ /* are in kilometers with three decimal places (to metres). */ /* e.g.: '30T 378.153 4798.667' */ /* */ /* When called as a command the return string is displayed. */ /* */ /* The formulae are derived from Redfearn's formulae, and give */ /* results which match maps to +/- 10m or better (i.e., as */ /* accurately as a 1:25,000 map can be read: +/- 5m in x or y). */ /* The results also exactly match those given by GeoTrans. */ /* --------------------------------------------------------------- */ -- 2007.08.27 initial version -- 2012.05.24 cosmetic updates gpsdata='EUR50 EUR79 WGS84' -- known datum arg lat lon datum . -- say lat ';' lon ';' datum ';' if left(lat, 1)='N' | left(lat, 1)='S' then do parse var lat ns +1 d':'m':'s if m='' then m=0 if s='' then s=0 lat=d+m/60+s/3600 if ns='S' then lat=-lat end if left(lon, 1)='E' | left(lon, 1)='W' then do parse var lon ew +1 d':'m':'s if m='' then m=0 if s='' then s=0 lon=d+m/60+s/3600 if ew='W' then lon=-lon end if \datatype(lat, 'n') | lat>90 | lat<-90 then call done '* Bad latitude:' lat if \datatype(lon, 'n') | lon>180 | lon<-180 then call done '* Bad longitude:' lon if lon=180 then lon=-180 -- algorithms may need this if lat>84 then call done '* Latitude too northerly for UTM' if lat<-80 then call done '* Latitude too southerly for UTM' if datum='' then datum='WGS84' if wordpos(datum, gpsdata)=0 then call done '* Unknown datum:' datum numeric digits 16 -- max possible for RxMath functions -- constants select when datum='WGS84' then do erad=6378137 -- ellipsoid radius ees=0.00669438 -- ellipsoid eccentricity squared end when datum='EUR50' | datum='EUR79' then do erad=6378388 -- Hayford/International ees=0.00672267 end end k0=0.9996 -- bias pi=3.1415926535897932 -- math -- load function package... if RxFuncQuery("MathLoadFuncs") then do call RxFuncAdd "MathLoadFuncs","rxmath","MathLoadFuncs" call MathLoadFuncs end -- determine the zone and zone origin longitude zone=trunc((lon+180)/6)+1 -- several ugly special cases in northern latitudes if lat>=56 & lat<64 & lon>=3 & lon<12 then zone=32 -- Norway if lat>=72 & lat<84 then do -- Arctics if lon<0 then nop else if lon<9 then zone=31 else if lon<21 then zone=33 else if lon<33 then zone=35 else if lon<42 then zone=37 end origin=(zone-1)*6-180 +3 -- +3 is origin offset in zone -- start calculating, using Redfearn's formulae ... ees2=ees**2 ees3=ees**3 M = erad*((1-ees/4 - 3*ees2/64 - 5*ees3/256) *d2r(lat), - (3*ees/8 + 3*ees2/32 + 45*ees3/1024)*rxcalcsin(2*lat), + (15*ees2/256 + 45*ees3/1024)*rxcalcsin(4*lat), - (35*ees3/3072)*rxcalcsin(6*lat)) eeps=ees/(1-ees) -- ee prime squared N = erad/rxcalcsqrt(1-ees*rxcalcsin(lat)**2); T = rxcalctan(lat)**2; C = eeps * rxcalccos(lat)**2; A = rxcalccos(lat)*d2r(lon-origin); utmx = (k0*N*(A+(1-T+C)*A**3/6, + (5-18*T+T*T+72*C-58*eeps)*A**5) + 500000) utmy = (k0*(M+N*rxcalctan(lat)*(A*A/2+(5-T+9*C+4*C*C)*A**4/24, + (61-58*T+T*T+600*C-330*eeps)*A**6/720))) -- add offset if southern hemisphere if lat<0 then utmy=utmy+10000000 utmx=format(utmx/1000,,3) -- to km, 3 decimals utmy=format(utmy/1000,,3) -- .. -- say 'UTM:' utmx utmy zone=zone||zonechar(lat) -- complete zone designator call done zone utmx utmy done: parse source . how . if how='FUNCTION' then exit arg(1) say arg(1) exit /* --------------------------------------------------------------- */ /* convert degrees to radians */ /* --------------------------------------------------------------- */ d2r: procedure expose pi return arg(1)*pi/180 /* --------------------------------------------------------------- */ /* Convert latitude to UTM zone character */ /* --------------------------------------------------------------- */ zonechar: procedure arg lat -- in degrees lat=lat+80 -- now 0 -> 164 lat=trunc(lat/8) -- now 0 -> 20, fraction removed if lat=20 then lat=19 -- N72 to N84 is all one band letters='CDEFGHJKLMNPQRSTUVWX' return substr(letters, lat+1, 1)