Base SAS has some GIS-type functions (GIS = Geographic Information System), the most useful one is beyond any doubt the GEODIST function that returns the geodetic distance between two latitude and longitude coordinates.
Remember that Euclidean distances won’t work on earth or any other (near-)spherical surface. On earth you need the geodesic distance (the great-circle distance). See this article about how to calculate geodesic distance using GEODIST or PROC DISTANCE.
However, there is no function for the MIDPOINT which is the half-way point along a great circle path between two points / positions. I present a little program below to calculate (the latitude-longitude coordinates of) this midpoint. The formula used is for calculating the midpoint based on a spherical earth (ignoring ellipsoidal effects) – which is accurate enough for most purposes… It’s hence a formula from spherical trigonometry.
This program can be useful for geographers, transport economists, transport scientists, shipping companies, airline companies and the like.
The program is using positions of two objects; objects can be cars, vessels (ships), airplanes, persons, …
PROC DATASETS library=work nolist;
delete positions_two_objects_in_degrees / memtype=DATA; run;
delete same_position_in_degrees_radians / memtype=DATA; run;
delete add_geodist_between_object1_n_o2 / memtype=DATA; run;
delete add_midpoint / memtype=DATA; run;
delete DS_with_all_variables / memtype=DATA; run;
QUIT;
data work.positions_two_objects_in_degrees;
latitude1=80; longitude1=-121; latitude2=-80; longitude2=23; output;
latitude1=45; longitude1=-175; latitude2=-80; longitude2=70; output;
run;
/* conversion from degrees to radians */
/* The reason for converting degrees to radians in below program is */
/* because the COS and SIN functions require an argument specified in radians */
/* and, also, the ATAN2 returns the result in radians */
data work.same_position_in_degrees_radians(drop=i);
set work.positions_two_objects_in_degrees;
array degrees{4} latitude1 longitude1 latitude2 longitude2;
array radians{4} latRadia1 longRadia1 latRadia2 longRadia2;
pi = constant('pi');
do i = 1 to dim(degrees);
radians(i) = (degrees(i)*pi)/180; /* conversion from degrees to radians */
end;
run;
data work.add_geodist_between_object1_n_o2;
set work.same_position_in_degrees_radians;
/* conversion from miles to nautical miles: */
/* for an approximate result, divide the length value in miles by 1,151 */
Distance = geodist(latRadia1,longRadia1,latRadia2,longRadia2,'RM') / 1.151;
run;
data work.add_midpoint(drop=B:);
set work.add_geodist_between_object1_n_o2;
/* What are the coordinates (latitude, longitude) of the Midpoint?: */
/* This is the half-way point along a great circle path between the two points. */
Bx = cos(latRadia2) * cos(longRadia2-longRadia1);
By = cos(latRadia2) * sin(longRadia2-longRadia1);
latRMid = atan2(sin(latRadia1) + sin(latRadia2),
sqrt( (cos(latRadia1)+Bx)*(cos(latRadia1)+Bx) + By*By ) );
longRMid = longRadia1 + atan2(By, cos(latRadia1) + Bx);
/* FOR GEODIST : */
/* If the LATITUDE value is expressed in degrees, it must be between 90 and –90. */
/* If the LATITUDE value is expressed in radians, it must be between pi/2 and –pi/2. */
/* If the LONGITUDE value is expressed in degrees, it must be between 180 and –180. */
/* If the LONGITUDE value is expressed in radians, it must be between pi and –pi. */
latRMid =max(-pi/2,latRMid) ; latRMid =min(pi/2,latRMid) ;
longRMid=max(-pi ,longRMid); longRMid=min(pi ,longRMid);
latDMid =(latRMid *180)/pi; /* conversion from radians to degrees */
longDMid=(longRMid*180)/pi; /* conversion from radians to degrees */
run;
data work.DS_with_all_variables(replace=yes repempty=no);
LENGTH result $ 7;
set work.add_midpoint;
Distance_Position1_to_MidPoint = geodist(latRadia1,longRadia1,latRMid,longRMid,'RM') / 1.151;
Distance_Position2_to_MidPoint = geodist(latRadia2,longRadia2,latRMid,longRMid,'RM') / 1.151;
if abs(Distance_Position1_to_MidPoint-Distance_Position2_to_MidPoint)/Distance > 0.025
then result='WRONG';
else result='CORRECT';
run;
PROC DATASETS library=work nolist;
modify DS_with_all_variables (label='Data set with coordinates and midpoint and distances');
label result = 'midpoint formula seems CORRECT (wanted) or WRONG (not wanted)';
label latitude1 = 'latitude of object_1 in degrees';
label longitude1 = 'longitude of object_1 in degrees';
label latitude2 = 'latitude of object_2 in degrees';
label longitude2 = 'longitude of object_2 in degrees';
label latRadia1 = 'latitude of object_1 in radians';
label longRadia1 = 'longitude of object_1 in radians';
label latRadia2 = 'latitude of object_2 in radians';
label longRadia2 = 'longitude of object_2 in radians';
label latRMid = 'latitude of the midpoint in radians';
label longRMid = 'longitude of the midpoint in radians';
label latDMid = 'latitude of the midpoint in degrees';
label longDMid = 'longitude of the midpoint in degrees';
label Distance = 'distance between both objects in nautical miles';
label Distance_Position1_to_MidPoint = 'distance between object_1 and MidPoint in nautical miles';
label Distance_Position2_to_MidPoint = 'distance between object_2 and MidPoint in nautical miles';
label pi = 'the constant pi (the ratio of the circumference of any circle to the diameter of that circle)';
run;
QUIT;
/* end of program */
SAS has a number of other features to support mapping and analytics related to geospatial data:
Enjoy reading and who knows, it may come in handy one day.
