I have a task that requires conversion of Latitude, Longitude extracted from a JPEG file (using ExifTool) to Lambert conformal conic (LCC) projection Easting and Northing. A search of the web didn't come up with much in the way of programs, so I researched the formulas and made my own. I've found it convenient to package functions in stand alone scripts then instantiate them from other scripts communicating via StdIn and StdOut. So what I've loaded here is two scripts, a LambertPipe (feed in latitude,longitude and get back Easting and Northing), and a LambertPipeTest which is a bare bones script to drive the LambertPipe. In real use, LambertPipe is driven by a much larger script that does more useful things.
Understand that there is no explanation of the LCC projection, you will need to research it separately if you don't know what it is. If you do so, you will realise the script here is specific to GDA94 used by the NSW government, and you will have to replace the parameters with ones appropriate to your usage.
These are not supposed to be shining examples of well written code. I'm publishing them because someone else might have the same problem, this could help them out.
LambertPipe.vbs
Notes:
1. this can be run from the command prompt if you wish.
2. The lat/lon input is very flexible, included in the output is a standardised lat/lon suitable for coordinates for KML files for input to Google Earth
3. There is also a reverse Easting, Northing to lat/lon conversion
Option Explicit
'=========================================================================
' This program is designed to be instantiated from another VBscript by:
' Set ExecObj = WshShell.Exec("cscript /B <path>\lambertpipe.vbs")
' input via ExecObj.StdIn and output via ExecObj.StdOut.
' The program prompts for input by writing ? to StdOut and expects as input
' EITHER a line of the form D<tab><Latitude><tab><Longitude>
' where Lat/Lon format is flexible
' The output is a line <Easting><tab><Northing><tab><Latitude><tab><Longitude>
' where Lat/Lon are in a format suitable for Google Earth KML files
' OR a line of the form L<tab><Easting><tab><Northing>, the output is a
' line <Latitude><tab><Longitude>
' The program is then ready for another line of input (no prompt given)
' The program terminates if it receives an empty line
'=========================================================================
'
' Configure a NSW GDA94 Lambert Conic Conformal Projection
' --------------------------------------------------------
'
' falseEasting & falseNorthing are just an offset in meters added to the final
' coordinate calculated.
'
' longOfOrigin & LatOfOrigin are the "center" latitude and longitude of the
' area being projected. All coordinates will be calculated in meters relative
' to this point on the earth.
'
' firstStdParallel & secondStdParallel are the two lines of longitude (that
' is they run east-west) that define where the "cone" intersects the earth.
' Simply put they should bracket the area being projected.
'
' find out more from the web
'
' Calculation method based on:
' European Petroleum Survey Group Guidance Note Number 7, part 2
' "Coordinate Conversions and Transformations including Formulas"
' Revised October 2004
'
'=========================================================================
'
' Constants from "NSW GDA94 Lambert Conformal Conic Projection Policy"
'
const StrFirstStdParallel = "-30 45 00" ' Latitude of 1st std parallel
const StrSecondStdParallel = "-35 45 00" ' Latitude of 2nd std parallel
const StrLatOfOrigin = "-33 15 00" ' Latitude of Origin
const StrLongOfOrigin = "147 00 00" ' Longitude of Origin
const StrFalseEasting = "9 300 000 m"
const StrFalseNorthing = "4 500 000 m"
'
' Constants from Geocentric Datum of Australia Technical Manual
' Reference Ellipsoid: Geodetic Reference System 1980 (GRS80) ellipsoid with a
' semi-major axis (a) of 6 378 137 metres exactly
' and an inverse flattening (l/f) of 298.257 222 101.
'
Const a = 6378137.0
Const InvF = 298.257222101
'
Dim StdIn,StdOu, work, rec
'
' Calculation variables
'
Dim PIon4, DtoR
dim f ' flattening
dim e2 ' eccentricity squared = 2f - f*f
dim e ' eccentricity
'
dim falseEasting, falseNorthing
dim phi1, phi2, phio, lamdao
dim m1, m2, t1, t2, t0, n, bigF, rbigF
'
dim lat, lon, phi, lamda, t, r, theta, lccEasting, lccNorthing
dim rRev, thetaRev, tRev, phiTrial, lv
'
PIon4 = Atn(1) ' Calculate the value of pi/4.
DtoR = PIon4 /45 ' degrees to radians
'
' Calculation of Derived Ellipsoid Parameters
'
f = 1/Cdbl(InvF)
e2 = 2 * f - f * f
e = Sqr(e2)
'
' Calculation of Lambert Projection variables related to NSW GDA94
'
falseEasting = StrFalse2Dec (StrFalseEasting)
falseNorthing = StrFalse2Dec (StrFalseNorthing)
phi1 = deg2rad(DegStr2DegDec(StrFirstStdParallel))
phi2 = deg2rad(DegStr2DegDec(StrSecondStdParallel))
phio = deg2rad(DegStr2DegDec(StrLatOfOrigin))
lamdao = deg2rad(DegStr2DegDec(StrLongOfOrigin))
m1 = cos(phi1) / sqr(( 1 - e2*sin(phi1)*sin(phi1)))
m2 = cos(phi2) / sqr(( 1 - e2*sin(phi2)*sin(phi2)))
t1 = tan((PIon4)-(phi1/2)) / (( ( 1 - e*sin(phi1) ) / ( 1 + e*sin(phi1) )) ^ (e/2))
t2 = tan((PIon4)-(phi2/2)) / (( ( 1 - e*sin(phi2) ) / ( 1 + e*sin(phi2) )) ^ (e/2))
t0 = tan((PIon4)-(phio/2)) / (( ( 1 - e*sin(phio) ) / ( 1 + e*sin(phio) )) ^ (e/2))
n = (log(m1)-log(m2)) / (log(t1)-log(t2))
bigF = m1/(n*(t1 ^ n))
rbigF = a*bigF*(t0 ^ n)
'
' Calculation related to a particular input
'
Set StdIn = WScript.StdIn
Set StdOu = WScript.StdOut
StdOu.Write "?"
rec = StdIn.ReadLine
Do while len(rec) > 0
work = split(rec,chr(9))
if work(0) = "D" then
'
' Input Lat, Lon as a tab separated pair
' Output Easting, Northing, Lat, Lon as a tab separated line
'
lat = DegStr2DegDec(trim(work(1)))
lon = DegStr2DegDec(trim(work(2)))
phi = deg2rad(lat) ' Latitude to convert
lamda = deg2rad(lon) ' Longitude to convert
t = tan((PIon4)-(phi /2)) / (( ( 1 - e*sin(phi ) ) / ( 1 + e*sin(phi ) )) ^ (e/2))
r = a*bigF*(t ^ n)
theta = n*(lamda - lamdao)
'
lccEasting = falseEasting + r*sin(theta)
lccNorthing = falseNorthing + rbigF - r*cos(theta)
StdOu.Write lccEasting & chr(9) & lccNorthing & chr(9) & lat & chr(9) & lon & VbNewLine
elseif work(0) = "L" then
'
' Reverse calculation
' Input Easting, Northing as a tab separated pair
' Output Lat, Lon as a tab separated pair
'
lccEasting = Cdbl(trim(work(1)))
lccNorthing = Cdbl(trim(work(2)))
rRev = sqr((lccEasting - falseEasting) ^ 2 + (rbigF - (lccNorthing - falseNorthing)) ^ 2) * sgn(n)
thetaRev = atn((lccEasting - falseEasting)/(rbigF - (lccNorthing - falseNorthing)))
tRev = (rRev/(a * bigF)) ^ (1/n)
phi = Cdbl(2 * (PIon4 - atn(tRev)))
phiTrial = Cdbl(10) ' ridiculous value
for lv = 0 to 9 ' do the loop until convergence or 10 times
phiTrial = phi
phi = 2 * (PIon4 - atn(tRev * (((1 - e * sin(phiTrial))/(1 + e * sin(phiTrial))) ^ (e / 2))))
if phiTrial/phi = 1 then exit for
next
lamda = thetaRev / n + lamdao
StdOu.Write rad2deg(phi) & chr(9) & rad2deg(Lamda) & VbNewLine
end if
'
rec = StdIn.ReadLine
loop
'
' ================================================
'
Function deg2rad (deg)
deg2rad = deg * DtoR
End Function
Function rad2deg (rad)
rad2deg = rad / DtoR
End Function
'
' Convert a string representing false northing/easting to a number
'
Function StrFalse2Dec (FalseStr)
' no checks, convenience conversion of known string constant to a double
Dim i, x, z
x=""
for i=1 to len(FalseStr)
z = mid(FalseStr,i,1)
if z => "0" and z <= "9" then x = x & z
next
StrFalse2Dec = Cdbl(x)
End Function
Function DegStr2DegDec (DegStr)
'
' Convert a string representing Latitude, Longitude to decimal degrees (S,W => -ve)
' input can be one to three numbers separated by any non numeric, representing
' degrees, optional minutes, optional seconds. Only the last (or only) number
' can have decimal places. South or West can be designated EITHER by an S or W
' as the first or last character, OR the degrees preceded immediately by -
'
Dim Degs,Mins,Secs
Dim LocStr, Wstr, Signed, Positive, DecPt, Found, i, z
LocStr = DegStr + Chr(0) ' ensures final number is captured
Degs=0:Mins=0:Secs=0
Signed = False
Positive = True
DecPt = False
Found = 0
Wstr = ""
for i = 1 to len(LocStr)
z = mid(LocStr,i,1)
if (z => "0" and z <= "9") or z = "." then
if len(Wstr) = 0 then
if Found = 3 then
DegStr2DegDec = 999 ' An error - too many numbers
Exit Function
end if
if DecPt then
DegStr2DegDec = 996 ' only last num can be non integer
Exit Function
end if
Found = Found + 1
if z = "." then Wstr "0"
end if
Wstr = Wstr & z
if z = "." then DecPt = True
else
if len(Wstr) > 0 then
if not IsNumeric(Wstr) then
DegStr2DegDec = 998 ' An error - can't process a number
Exit Function
end if
if Found = 1 then Degs = Cdbl(Wstr)
if Found = 2 then Mins = Cdbl(Wstr)
if Found = 3 then Secs = Cdbl(Wstr)
Wstr = ""
end if
' not part of a number
if z = "-" then
z = mid(LocStr,i+1,1)
if not Signed and Found = 0 then
if (z => "0" and z <= "9") or z = "." then
Signed = True
Positive = False
else
DegStr2DegDec = 997 ' - must be start of a number
Exit Function
end if
end if ' ignore hyphen in middle, might be a number separator
end if
if (ucase(z) = "S" or ucase(z) = "W") and (i = 1 or i = len(LocStr)) then
if not Signed then
Positive = False
Signed = True
else
DegStr2DegDec = 995 ' An error - too many signers
Exit Function
end if
end if
' ignore all other separators
end if
next
if Found = 0 then
DegStr2DegDec = 994 ' An error - no numbers
Exit Function
end if
if Degs => 360 or (Degs > 180 and Signed) or Mins => 60 or Secs => 60 then
DegStr2DegDec = 993 ' An error - something out of range
Exit Function
end if
if Positive then
DegStr2DegDec = Degs + Mins/60 + Secs/3600
else
DegStr2DegDec = -Degs - Mins/60 - Secs/3600
end if
End Function
LambertPipeTest.vbs
Demonstrates how to run LambertPipe as a dependent script. The two scripts must be in the same folder.
Option Explicit
'
' testbed for lambertpipe.vbs
' follow the prompts
'
Dim fso, ProgFolder, PipeCommand, WshShell, oExec, CharIn, Rec
Set fso = CreateObject("Scripting.FileSystemObject")
' find out where we are
ProgFolder = fso.GetParentFolderName(Wscript.ScriptFullName)
PipeCommand = "cscript /B " & ProgFolder & "\lambertpipe.vbs"
Set WshShell = CreateObject("Wscript.Shell")
Set oExec = WshShell.Exec(PipeCommand)
' do this to synch with the other program
CharIn = oExec.StdOut.Read(1)
do until CharIn = "?":CharIn = oExec.StdOut.Read(1):loop
' a normal "loop until end"
Rec = GetInput()
do while len(Rec) > 0
' WriteLine doesn't work!
oExec.StdIn.Write Rec & VbNewLine
Rec = oExec.StdOut.ReadLine
Wscript.echo Rec
Rec = GetInput()
loop
' tell lambertpipe we are finished with it
oExec.StdIn.Write VbNewLine
Set oExec = Nothing
Function GetInput()
Dim Inp
Do
Wscript.StdOut.WriteLine "D for Lat/Lon to Lambert,"
Wscript.StdOut.WriteLine "L for Lambert to Lat/Lon,"
Wscript.StdOut.Write "Q to Quit>"
Inp = Wscript.StdIn.ReadLine
If Inp = "Q" then
GetInput = ""
Exit Function
End if
If Inp = "D" or Inp = "L" then Exit Do
loop
If Inp = "D" then
Wscript.StdOut.Write "Latitude >"
Inp = Inp & Chr(9) & Wscript.StdIn.ReadLine
Wscript.StdOut.Write "Longitude>"
else
Wscript.StdOut.Write "Easting >"
Inp = Inp & Chr(9) & Wscript.StdIn.ReadLine
Wscript.StdOut.Write "Northing>"
end if
GetInput = Inp & Chr(9) & Wscript.StdIn.ReadLine
End Function