Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Home
- /
- SAS Communities Library
- /
- MidPoint: the half-way point along a great circle path between two poi...

Options

- RSS Feed
- Mark as New
- Mark as Read
- Bookmark
- Subscribe
- Printer Friendly Page
- Report Inappropriate Content

- Article History
- RSS Feed
- Mark as New
- Mark as Read
- Bookmark
- Subscribe
- Printer Friendly Page
- Report Inappropriate Content

Views
4,196

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:

- See for example the beautiful interactive maps in the SAS Visual Analytics Gallery
- Base SAS also has multiple mapping procedures to work with map polygons (like PROC SGMAP, PROC GINSIDE, PROC GPROJECT and so on). There is also the very useful %CENTROID autocall macro to retrieve the centroids of (map) polygons.
- SAS/STAT has several procedures for spatial (regression) models.
- SAS/ETS can do spatial econometrics.

Enjoy reading and who knows, it may come in handy one day.

Comments

06-20-2022
01:38 PM

- Mark as Read
- Mark as New
- Bookmark
- Permalink
- Report Inappropriate Content

06-20-2022
01:38 PM

Registration is open! SAS is returning to Vegas for an AI and analytics experience like no other! Whether you're an executive, manager, end user or SAS partner, SAS Innovate is designed for everyone on your team. Register for just $495 by 12/31/2023.

**If you are interested in speaking, there is still time to submit a session idea. More details are posted on the website. **

Data Literacy is for **all**, even absolute beginners. Jump on board with this free e-learning and boost your career prospects.

Article Labels

Article Tags