Geographic Functions, any takers?

[Ben Rubinstein]

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