Download or view CoordinateConversions.frink in plain text format
// Conversions between different geodetic coordinate systems.
// References:
// http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.HTM
// Requires use of a datum data structure from another file to
// represent the shape of the earth in that datum.
use Datum.frink
// Convert a set of coordinates from UTM to Lat/Long.
// easting and northing are given as dimensionless numbers
// (consider changing this to meters?), zone indicates the UTM zone as a
// string, e.g. "13E". (The letter is ignored,) and the datum is a Datum
// object (see Datum.frink)
// Returns: [lat, long]
UTMToLatLong[easting, northing, zone, datum is Datum] :=
{
y = northing
x = easting - 500000 // Relative to central meridian
k0 = 0.9996 // Scale along long0
// Meridional arc (that looks spelled funny)
M = y/k0 m
// println["M = $M"]
// Footprint latitude
// Add more terms to this?
mu = M/(datum.a (1 - datum.e^2/4 - 3 datum.e^4/64 - 5 datum.e^6/256))
// println["mu = $mu"] // format[mu, "degrees", 5]]
e1 = (1 - (1 - datum.e^2)^(1/2)) / (1 + (1 - datum.e^2)^(1/2))
// println["e1 = $e1"]
// Add more terms to these?
J1 = (3 e1/2 - 27 e1^3/32)
J2 = (21 e1^2/16 - 55 e1^4/32)
J3 = (151 e1^3/96)
J4 = (1097 e1^4/512)
fp = mu + J1 sin[2 mu] + J2 sin[4 mu] + J3 sin[6 mu] + J4 sin[8 mu]
// println["fp = " + format[fp, "degrees", 5]]
// Now go to lat/long
C1 = datum.eprime^2 cos[fp]^2
T1 = tan[fp]^2
// This is the same as rho in the forward conversion formulas, but
// calculated for fp instead of lat.
R1 = datum.a (1 - datum.e^2) / (1 - datum.e^2 sin[fp]^2)^(3/2)
// This is the same as nu in the forward conversion formulas above, but
// calculated for fp instead of lat.
N1 = datum.a / (1 - datum.e^2 sin[fp]^2)^(1/2)
D = x / (N1 k0/m)
// Coefficients for latitude
Q1 = N1 tan[fp]/R1
Q2 = (D^2/2)
Q3 = (5 + 3 T1 + 10 C1 - 4 C1^2 - 9 datum.eprime^2) D^4/24
Q4 = (61 + 90 T1 + 298 C1 + 45 T1^2 - 3 C1^2 - 252 datum.eprime^2) D^6/720
// println["$Q1\t$Q2\t$Q3\t$Q4"]
lat = fp - Q1 (Q2 - Q3 + Q4)
// println["lat= " + format[lat, "degrees", 5]]
// Coefficients for longitude
Q5 = D
Q6 = (1 + 2 T1 + C1) D^3/6
Q7 = (5 - 2 C1 + 28 T1 - 3C1^2 + 8 datum.eprime^2 + 24 T1^2) D^5/120
// Get central meridian
long0 = UTMZoneToLong[zone]@1
long = long0 + (Q5 - Q6 + Q7) / cos[fp]
// println["long= " + format[long, "degrees", 5]]
return [lat,long]
}
// Convert a set of coordinates from Lat/Long to UTM. The datum is an
// object of type Datum (see Datum.frink)
// Equation numbers listed in comments are for reference to Snyder,
// _Map Projections, A Working Manual_
// Returns: [easting, northing, zone]
LatLongToUTM[lat, long, datum is Datum] :=
{
k0 = 0.9996 // Scale along lat
e = datum.e
// Calculate the meridional arc
// TODO: Add more terms
M = datum.a ((1 - e^2/4 - 3 e^4/64 - 5 e^6/256) lat -
(3 e^2/8 + 3 e^4/32 + 45 e^6/1024) sin[2 lat] +
(15 e^4/256 + 45 e^6/1024) sin[4 lat] -
(35 e^6/3072) sin[6 lat])
ep2 = e^2/(1-e^2) // 8-12
N = datum.a/(1 - e^2 sin[lat]^2)^(1/2) // 4-20
T = tan[lat]^2 // 8-13
C = ep2 cos[lat]^2 // 8-14
long0 = centralMeridianLongitude[long]
A = (long - long0) cos[lat] // 8-15
// Calculate (false) easting
// Eq. 8-9
x = k0 N ( A + (1-T+C) A^3/6 + (5 - 18T + T^2 + 72C - 58 ep2) A^5 / 120)
// Eq. 8-10
y = k0 * (M + N tan[lat] (A^2/2 + (5 - T + 9C + 4C^2) A^4/24 +
(61 - 58T + T^2 + 600C - 330 ep2) A^6/720))
return [x/m+500000, y/m, LatLongToUTMZone[lat, long]]
}
// This converts a UTM zone to a longitude triplet.
// The zone is a string like "13E". The letters are currently ignored.
// Returns:
// [longW, longCenter, longE] representing the west, center, and east meridian
// of a zone.
UTMZoneToLong[zone] :=
{
[zoneNum] = zone =~ %r/^\s*(\d+)/ // Parse out number
zoneNum = parseInt[zoneNum]
longL = -180 degrees + 6 degrees (zoneNum-1)
return [longL, longL + 3 degrees, longL + 6 degrees]
}
// This finds the central meridian of the nearest UTM zone for a given
// longitude.
centralMeridianLongitude[long] :=
{
long1 = (long/degrees) + 180
west = floor[long1/6]*6 - 180 // West side of zone.
center = west + 3
return center degrees
}
// This converts a lat/long to the UTM zone in which it is located.
// This may give unpredictable results if taken exactly on the border between
// two zones, so don't do that.
// This also returns a zone letter for latitude zones C-W. If the latitude
// is outside this zone, this will return a question mark for the zone letter,
// as a warning that things are getting really unsafe for UTM coordinates.
LatLongToUTMZone[lat, long] :=
{
long1 = ((long/degrees)+180)/6
zoneNum= ceil[long1]
// Now find latitude letter. Some letters like I and O are not used.
letters = ["C","D","E","F","G","H","J","K","L","M","N","P","Q","R","S","T","U","V","W"]
if (lat < 72 degrees) and (lat > -80 degrees)
{
lat1 = ((lat/degrees) + 80)/8 // C starts at 80 south, bands 8 deg tall
latband = floor[lat1]
zoneLetter = letters@latband
} else
zoneLetter = "?"
return "$zoneNum$zoneLetter"
}
"CoordinateConversions included successfully"
Download or view CoordinateConversions.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19970 days, 10 hours, 24 minutes ago.