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 20236 days, 0 hours, 38 minutes ago.