On 30/08/2010 17:19, Ewan Richardson wrote:
> Has anyone written any libraries for Geographic functions within rev?
>
>
> Im thinking of a library for dealing with latitude and longitude, for
> example determining distance between points, bearings radian calculations
> ect.
Ewan,
Sorry for the delayed reply, just catching up with my mailing lists.
I did a project some years ago involving an embedded computer tracking GPS on
a moving vehicle. Looking at the code today, I see I started trying to do it
all by the book - and subsequently realised that since I was only going to be
interested in fine detail when the vehicle was close to any of the waypoints,
and the device was only going to used in Western Europe, I didn't actually
care about the true shape of the Earth etc!
However, for what it's worth, I find these two fragments of code; the first
looks jolly useful, but it appears that I don't actually use it in the final
system. The second is what's actually used, and produces a 'good enough'
result... if nothing else, the references in the comments might be useful.
Also a couple of utilities for working NMEA (the default 'language' of most
GPS units) - trivial but might save a few minutes.
I present these directly pasted from scripts, no checking, no warranty...
Ben
on initialiseGeoConstants
-- an array giving the approximage length in metres of one second of
longitude, at a given degree of latitude
-- source:
http://www.gis.unbc.ca/2003/courses/geog205/lectures/thegraticule/bottomframe.html
constant kSrcInfo = "0 111.32,10 109.64,20 104.65,30 96.49,40 85.39,50
71.70,60 55.80,70 38.19,80 19.39,90 0.00"
put empty into aTemp
put empty into bTemp
repeat for each item tDataPoint in kSrcInfo
put (((word 2 of tDataPoint) * 1000) / 3600) into aTemp[word 1 of
tDataPoint] -- km/degree to m/sec
put ((word 2 of tDataPoint) * 1000) into bTemp[word 1 of tDataPoint] --
km/degree to m/degree
end repeat
-- now try to fill in the degrees between, by crude linear interpolation
repeat with d = 1 to 89
get aTemp[d]
if it <> empty then next repeat
put trunc(d / 10) * 10 into d1
put d1 + 10 into d2
--
get (aTemp[d1] - aTemp[d2])
multiply it by (d2 - d)
put aTemp[d2] + (it / 10) into aTemp[d]
--
get (bTemp[d1] - bTemp[d2])
multiply it by (d2 - d)
put bTemp[d2] + (it / 10) into bTemp[d]
end repeat
set the customProperties["uLongSecLengthAtLatitude"] of this stack to aTemp
set the customProperties["uLongDegLengthAtLatitude"] of this stack to bTemp
end initialiseGeoConstants
--
-- %MACRO distance(lat1, lon1, lat2, lon2);
--
--%* calculates distance between two locations on earth; %* depends crucially
on the accurcy of the calculations!;
--
--%* input: latitude and longitude in radians; %* output: distance in meters,
as a data step expression; %* version: 26 June 1996, Willem Dekker;
--
--%Let Mile=1853.245; %* Nautical mile in meters; %Let
DegRad=180/3.141592653589793238; %* Degrees per Radian; %LET
MeterRad=60*&Mile*&DegRad; %* meters per radian;
--
--(&MeterRad * 2 * arsin(sqrt(2 - 2 * cos(&lat1) *cos(&lat2) * cos(&lon1 -
&lon2) - 2 * sin(&lat1) * sin(&lat2) ) / 2)); %mend;
-- source: http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
--
-- Presuming a spherical Earth with radius R (see below), and the locations of
the two points in spherical coordinates (longitude, latitude) are lon1, lat1
and lon2, lat2, then the Haversine Formula (Sinott 1984) will give
mathematically and computationally exact results. The intermediate result c is
the great circle distance in radians. The great circle distance d will be in
the same units as R. We provide the inverse tangent version of the Haversine
Formula:
--
-- dlon = lon2 - lon1
--
-- dlat = lat2 - lat1
--
-- a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
--
-- c = 2 * atan2( sqrt(a), sqrt(1-a) )
--
-- d = R * c
--
-- For R, use 3957 miles (6367 km), which is the radius of curvature at 30°
latitude; the use of other possible values results in very small differences.
Most computers require the arguments of trigonometric functions to be
expressed in radians. To convert lon1, lat1 and lon2, lat2 from degrees,
minutes, and seconds to radians, first convert them to decimal degrees. To
convert decimal degrees to radians, multiply the number of degrees by: pi/180
= 0.017453293 radians/degree. Inverse trigonometric functions return results
expressed in radians. To express c in decimal degrees, multiply the number of
radians by 180/pi = 57.295780 degrees/radian. (But be sure to multiply the
number of radians by R to get d.)
--
function rangeBetweenPoints tLat1, tLon1, tLat2, tLon2
constant kRadius = 6367000 -- radius of curvature at 30° latitude; in metres
local kDegreesToRadians
put pi / 180 into kDegreesToRadians
local a, c, iDeltaLat, iDeltaLon
multiply tLat1 by kDegreesToRadians
multiply tLon1 by kDegreesToRadians
multiply tLat2 by kDegreesToRadians
multiply tLon2 by kDegreesToRadians
put tLat2 - tLat1 into iDeltaLat
put tLon2 - tLon1 into iDeltaLon
put (sin(iDeltaLat / 2) ^ 2) + cos(tLat1) * cos(tLat2) * (sin(iDeltaLon /
2) ^ 2) into a
put 2 * atan2(sqrt(a), sqrt(1 - a)) into c
-- put round(kRadius * c) & return & "lat1:" && tLat1 & return & "lon1:" &&
tLon1 & return & "lat2:" && tLat2 & return & "lon2:" && tLon2 & return &
"Dlat" && iDeltaLat & return & "Dlon" && iDeltaLon & return & "a" && a &
return & "c" && c
return round(kRadius * c)
end rangeBetweenPoints
function nmeaChecksum tString --> tSum
local iSum
put 0 into iSum
repeat for each char c in tString
put iSum bitXor chartonum(c) into iSum
end repeat
return format("*%02X", iSum)
end nmeaChecksum
function nmeaLatLongToDecimal tNMEAlatitude, tHemisphere --> tDecimalLat
get offset(".", tNMEAlatitude)
if it = 0 then
cancelNonRevMessages
breakpoint
return empty
end if
put char (it - 2) to -1 of tNMEAlatitude into tMinutes
put char 1 to (it - 3) of tNMEAlatitude into tDecimalLat
add (tMinutes / 60) to tDecimalLat
-- put (tDecimalLat * 60) + tMinutes into tDecimalLat -- express in minutes
-- multiply tDecimalLat by 60 -- actually, seconds
if tHemisphere is in "SW" then put "-" before tDecimalLat
return tDecimalLat
end nmeaLatLongToDecimal
_______________________________________________
use-revolution mailing list
[email protected]
Please visit this url to subscribe, unsubscribe and manage your subscription
preferences:
http://lists.runrev.com/mailman/listinfo/use-revolution
