Latitude/Longitude Conversion to UTM in Visual FoxPro

And now, for something completely different:

A friend of mine works for Rogers Communications, a huge cable company in this part of Canada. He's involved in 'plant': that part of the business that deals with physical assets in the field i.e. the location of swith boxes, cables, routers, switches, hubs and the like. He wanted to know if I knew anything about conversion from latitude/longitude to UTM (Universal Transverse Mercator) Coordinate System. I said I didn't and we left it at that.

But then for fun I got poking around on the web and found a version in PHP done by Jim Studnicki. He, in turn, had translated this from a  function originally written in C by Jef Poskanzer. I spent some time converting Jim's PHP version to Visual FoxPro and after the usual minor problems I got it to work. I must say this is the first time I used the DTOR() (degrees to radians) function in VFP. (not that function is hard to mimic otherwise: pi/180).

Anyway, here it is. It ain't pretty, but if you need it, you can fix it up yourself. Once you have it up and running in VFP, test the results against any of the UTM conversion calculators on the web to verify your results.

 Hopefully it can help some other VFP developer in the future. This is what is so absolutely fabulous about the web: the free sharing of information. Ten years ago, getting an answer to this would have been a real problem. Now, anybody with half a brain can do this.

*Keith Hekker 2006-04-02
*Translated from a PHP function by Jim Studnicki
=ConvertToUTM(43,-80)

Procedure ConvertToUTM()
PARAMETERS Latitude, Longitude
* convert decimal geographic coordinates to UTM
*
* param latitude as decimal
* param longitude as decimal

cOldDecimals = SET("DECIMALS")
SET DECIMALS TO 6

* square of eccentricity of equatorial cross-section
ECCENTRICITY_SQUARED = 0.00669437999013

* eccentricity prime squared
ECC_PRIME_SQUARED = ECCENTRICITY_SQUARED /;
(1.0 - ECCENTRICITY_SQUARED)

* radius of Earth in meters
EQUATORIAL_RADIUS = 6378137.0

* scale factor
K0 = 0.9996

* make sure longitude is between -180 and 180
IF longitude < -180.0
longitude = longitude + 360.0
ENDIF
IF longitude > 180.0
longitude = longitude + 360.0
ENDIF

* get UTM letter
DO CASE
CASE latitude <= 84.0 and latitude >= 72.0
utmLetter = "X"
CASE latitude < 72.0 and latitude >= 64.0
utmLetter = "W"
CASE latitude < 64.0 and latitude >= 56.0
utmLetter = "V"
CASE latitude < 56.0 and latitude >= 48.0
utmLetter = "U"
CASE latitude < 48.0 and latitude >= 40.0
utmLetter = "T"
CASE latitude < 40.0 and latitude >= 32.0
utmLetter = "S"
CASE latitude < 32.0 and latitude >= 24.0
utmLetter = "R"
CASE latitude < 24.0 and latitude >= 16.0
utmLetter = "Q"
CASE latitude < 16.0 and latitude >= 8.0
utmLetter = "P"
CASE latitude < 8.0 and latitude >= 0.0
utmLetter = "N"
CASE latitude < 0.0 and latitude >= -8.0
utmLetter = "M"
CASE latitude < -8.0 and latitude >= -16.0
utmLetter = "L"
CASE latitude < -16.0 and latitude >= -24.0
utmLetter = "K"
CASE latitude < -24.0 and latitude >= -32.0
utmLetter = "J"
CASE latitude < -32.0 and latitude >= -40.0
utmLetter = "H"
CASE latitude < -40.0 and latitude >= -48.0
utmLetter = "G"
CASE latitude < -48.0 and latitude >= -56.0
utmLetter = "F"
CASE latitude < -56.0 and latitude >= -64.0
utmLetter = "E"
CASE latitude < -64.0 and latitude >= -72.0
utmLetter = "D"
CASE latitude < -72.0 and latitude >= -80.0
utmLetter = "C"
OTHERWISE
* returns "Z" if the latitude is outside the UTM limits of 84N to 80S
utmLetter = "Z" *
ENDCASE

lat_rad = DTOR(latitude)
long_rad = DTOR(longitude)
zone = INT((longitude + 180) / 6) + 1
IF latitude >= 56.0 and latitude < 64.0;
and longitude >= 3.0 and longitude < 12.0
zone = 32
ENDIF

* Special zones for Svalbard.
IF latitude >= 72.0 and latitude < 84.0
DO CASE
CASE longitude >= 0.0 and longitude < 9.0
zone = 31
CASE longitude >= 9.0 and longitude < 21.0
zone = 33
CASE longitude >= 21.0 and longitude < 33.0
zone = 35
CASE longitude >= 33.0 and longitude < 42.0
zone = 37
ENDCASE
ENDIF

* +3 puts origin in middle of zone
long_origin = (zone - 1) * 6 - 180 + 3
long_origin_rad = DTOR(long_origin)
N = EQUATORIAL_RADIUS / sqrt(1.0 - ECCENTRICITY_SQUARED *;
sin(lat_rad) * sin(lat_rad))
T = tan(lat_rad) * tan(lat_rad)
C = ECC_PRIME_SQUARED * cos(lat_rad) * cos(lat_rad)
A = cos(lat_rad) * (long_rad - long_origin_rad)

M = EQUATORIAL_RADIUS *;
((1.0 - ECCENTRICITY_SQUARED / 4 - 3 * ;
ECCENTRICITY_SQUARED * ECCENTRICITY_SQUARED / 64 - 5 *;
ECCENTRICITY_SQUARED * ECCENTRICITY_SQUARED *;
ECCENTRICITY_SQUARED / 256) * lat_rad - (3 * ;
ECCENTRICITY_SQUARED / 8 + 3 * ECCENTRICITY_SQUARED *;
ECCENTRICITY_SQUARED / 32 + 45 * ECCENTRICITY_SQUARED *;
ECCENTRICITY_SQUARED * ECCENTRICITY_SQUARED / 1024) * ;
sin(2 * lat_rad) + (15 * ECCENTRICITY_SQUARED *;
ECCENTRICITY_SQUARED / 256 + 45 * ECCENTRICITY_SQUARED *;
ECCENTRICITY_SQUARED * ECCENTRICITY_SQUARED / ;
1024) * sin(4 * lat_rad) - (35 * ECCENTRICITY_SQUARED *;
ECCENTRICITY_SQUARED * ECCENTRICITY_SQUARED / 3072) *;
sin(6 * lat_rad))

easting = K0 * N * (A + (1 - T + C) * A * A *;
A / 6 + (5 - 18 * T + T * T + 72 * C - 58 *;
ECC_PRIME_SQUARED) * A * A * A * A * A / 120) + 500000.0


northing = K0 * (M + N * tan(lat_rad) *;
(A * A / 2 + (5 - T + 9 * C;
+ 4 * C * C) * A * A * A * A / 24;
+ (61 - 58 * T + T * T + 600 * C -;
330 * ECC_PRIME_SQUARED) * A * A * A * A * A * A / 720))

IF latitude < 0.0
* 1e7 meter offset for southern hemisphere
northing = northing + 10000000.0
north = .F.
ELSE
north = .T.
ENDIF

? "Easting = " + TRANSFORM(easting)
? "Northing = " + TRANSFORM(northing)
? "North = " + TRANSFORM(north)
? "Zone = " + TRANSFORM(zone)
? "UTM Letter = " + TRANSFORM(utmletter)
SET DECIMALS TO &cOldDecimals
RETURN