This made me curious. Where would be the furthest place someone could travel to SGF17 from?
proc fcmp outlib=work.func.geo;
function radians(deg) group='Geospatial' label='Convert degrees to radians';
return ((deg/180)*constant('PI'));
endsub;
function degrees(rad) group='Geospatial' label='Convert radians to degrees';
return ((rad/constant('PI'))*180);
endsub;
function disrad(d, option $) group='Geospatial' label='Convert distance to unit';
if option='nm'
then return ((d*180*60)/constant('PI'));
if option='km'
then return ((d*180*60*1.852)/constant('PI'));
if option='miles'
then return ((d*180*60*1.150779)/constant('PI'));
endsub;
function antipodal(lat1, lon1, lat2, lon2) group='Geospatial' ;
tol=1e-9;
diflat=abs(lat1 - lat2);
diflon=abs(lon1 + lon2);
return ((diflat < tol) & (abs(mod(diflon,360)-180) < tol));
endsub;
function gcDistance(lat1, lon1, lat2, lon2) group='Geospatial' label='Great Circle Distance in Radians';
return (2*arsin(sqrt((sin((lat1-lat2)/2))**2 + cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))**2)));
endsub;
subroutine gcIntermediate(lat1, lon1, lat2, lon2, f, lat, lon) group='Geospatial' label='Calculate intermediate points along a fractional point along a Great Circle route';
outargs lat, lon;
if (antipodal(lat1, lon1, lat2, lon2))
then return;
if (lat1=lat2 and lon1=lon2)
then return;
d=gcDistance(lat1,lon1,lat2,lon2);
A=sin((1-f)*d)/sin(d);
B=sin(f*d)/sin(d);
x=A*cos(lat1)*cos(lon1) + B*cos(lat2)*cos(lon2);
y=A*cos(lat1)*sin(lon1) + B*cos(lat2)*sin(lon2);
z=A*sin(lat1) + B*sin(lat2);
lat=atan2(z,sqrt(x**2+y**2));
lon=atan2(y,x);
endsub;
quit;
options cmplib=(work.func);
data map;
set mapsgfk.world(drop=x y);
where cont^=97 and density<2;
x=radians(long);
y=radians(lat);
run;
proc gproject data=map out=right dupok project=none
longmin=0 longmax=4;
id id;
run;
proc gproject data=map out=left dupok project=none
longmin=-4 longmax=0;
id id;
run;
data right;
set right;
segment=segment+5000;
x=x-constant('PI')*2;
run;
data pmap;
set right left;
run;
proc gproject data=pmap out=pmap parmin=mapsgfk.projparm parmentry=world norangecheck radians parmout=work.projparm;
id id;
run;
data orlando_dist;
set mapsgfk.world_cities(rename=(lat=_lat long=_long) where=(city='Orlando' and MapIDName1='Florida'));
_lat=radians(_lat);
_long=radians(_long);
do until(done);
set mapsgfk.world_cities(where=(prxmatch('/^cities_/o',CtType))) end=done;
lat=radians(lat);
long=radians(long);
dist_km=disrad(gcDistance(_lat,_long,lat,long),'km');
geod_km=geodist(_lat,_long,lat,long,'rk');
output;
end;
run;
ods select none;
ods output ExtremeObs=orlando_dist_max(keep=high city_high MapIDName1_high idname_high lat_high long_high _lat_high _long_high rename=(high=distance lat_high=lat1 long_high=lon1 _lat_high=lat2 _long_high=lon2));
proc univariate data=orlando_dist nextrobs=1;
var dist_km;
id city idname MapIDName1 lat long _lat _long;
run;quit;
ods exclude none;
data orlando_dist;
set orlando_dist_max end=done;
f=0;
y=lat1;
x=lon1;
output;
do f=0.01 to 0.99 by 0.01;
call gcIntermediate(lat1,lon1,lat2,lon2,f,y,x);
output;
end;
f=1;
y=lat2;
x=lon2;
output;
run;
data orlando_dist;
length id $15;
set orlando_dist;
if x>0 then x=x-constant('PI')*2;
id=cats('gcD',int(distance));
run;
proc gproject data=orlando_dist out=orlando_dist parmin=mapsgfk.projparm parmentry=world norangecheck radians;
id id;
run;
data anno;
length function color $8 style $12 html $300;
retain xsys ysys '2' hsys '3' line 1 anno_flag 1 when 'a';
set orlando_dist;
by id notsorted;
if first.id then do;
function='move'; output;
end;
else do;
function='draw'; size=.25; color='red'; output;
end;
run;
proc print data=orlando_dist_max label;
title 'Furthest Great Circle Distance to Orlando, FL';
footnote 'Location data from mapsgfk.world_cities';
id CITY_High IDNAME_High;
var distance;
label IDNAME_High='Country' CITY_High='City' distance='Distance (km)';
run;
proc gmap map=pmap data=pmap anno=anno;
title 'Great Circle Route Map';
id id;
choro id / levels=1 nolegend des='';
run;
quit;
Now, obviously, there are a lot of assumptions here about how people are traveling, etc... but, it's just for fun. The map definitely deserves some better annotation as well, but I've been wanting to write something for displaying great circle routes for a while, so that's what I've done.
City
Country
Distance (km)
Geraldton
Australia
18442.7
... View more
![SASGF 2020 [Virtual] Attendee](/html/@A53A708409112314CE528E074D4C8C1F/badge_icons/badge_sasgf20_attendee.png "SASGF 2020 [Virtual] Attendee"){kind=icon}
