%let name=magnetic_north_pole;

/* 
Set your current-working-directory (to read/write files), if you need to ...
%let rc=%sysfunc(dlgcdir('c:\someplace\public_html')); 
*/
filename odsout '.';

/*
Using data from:
https://www.ngdc.noaa.gov/geomag/GeomagneticPoles.shtml
https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy
*/

/*
filename dataurl url "https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy";
*/
filename dataurl "NP.xy.txt";

data magnetic_pole_locations;
infile dataurl firstobs=1 pad;
format Long Lat comma8.3;
label Lat='Latitude (degrees)';
label Long='Longitude (degrees east)';
input Long Lat Year;
run;

data magnetic_pole_locations; set magnetic_pole_locations;
length year_label recent_label $6;
if mod(year,50)=0 then year_label=trim(left(year))||'a0a0'x;
if mod(year,10)=0 and year>2000 then recent_label=trim(left(year))||'a0a0'x;
run;

data magnetic_pole_locations; set magnetic_pole_locations;
/* Convert long from 0 to 360 values to -180 to 180 values, to use with geodist() */
if long>180 then long=long-360;
label yearly_distance='Distance from previous year (miles)';
format yearly_distance comma8.1;
yearly_distance=geodist(lat, long, lag(lat), lag(long), 'DM');
run;

proc sql noprint;
select min(year) into :minyear separated by ' ' from magnetic_pole_locations;
select max(year) into :maxyear separated by ' ' from magnetic_pole_locations;
quit; run;

data rings;
do lat_ring = 70 to 90 by 5;
 do long_ring = -180 to 180 by 5;
  output;
  end;
 end;
run;

ods escapechar='^';
data rings; set rings;
length label_ring $25;
if long_ring=120 then label_ring=trim(left(lat_ring))||"^{unicode '00ba'x}";
run;

data country_names;
input lat_country long_country label_country $ 24-80;
label_country=trim(left(label_country));
datalines;
68.301721 -153.8469162 Alaska
67.0276347 161.2927574 Russia
66.273992 -104.4765288 Canada
76.2910011 -45.3329711 Greenland
64.7710718 -17.7811442 Iceland
68.7642254  23.4749876 Norway
65.4826228  19.2217949 Sweden
65.436497   29.9555922 Finland
65.985969   60.3296276 Russia
;
run;

data my_data; set magnetic_pole_locations rings country_names;
run;


ODS LISTING CLOSE;
ODS HTML path=odsout body="&name..htm" 
 (title="Magnetic North Pole") 
 style=htmlblue;

ods graphics /
 noscale /* if you don't use this option, the text will be resized */
 imagemap tipmax=2500 
 imagefmt=png imagename="&name"
 width=750px height=900px noborder; 

title1 c=gray33 h=18pt "Shift in magnetic north pole (yearly position, &minyear-&maxyear)";
footnote c=gray h=12pt "Data source: https://www.ngdc.noaa.gov/geomag/GeomagneticPoles.shtml";

ods graphics / width=750 height=900;

proc sgmap plotdata=my_data noautolegend;
esrimap url='http://services.arcgisonline.com/arcgis/rest/services/Polar/Arctic_Ocean_Base/MapServer';
text x=long_country y=lat_country text=label_country / textattrs=(color=gray55 size=14pt) 
 position=center tip=none;
series x=long_ring y=lat_ring / group=lat_ring lineattrs=(color=gray55) tip=none;
text x=long_ring y=lat_ring text=label_ring / textattrs=(color=gray55 size=12pt) 
 position=left tip=none;
scatter x=long y=lat / markerattrs=(color=red symbol=circlefilled size=5px) 
 /* you'll need Viya 3.5 to use the tip= option */
 tip=(year lat long yearly_distance);
text x=long y=lat text=year_label / textattrs=(color=black size=9pt) position=left tip=none;
text x=long y=lat text=recent_label / textattrs=(color=blue size=9pt) position=left tip=none;
run;

proc sort data=magnetic_pole_locations out=table_data;
by descending year;
run;

title1 c=gray33 h=18pt "Shift in magnetic north pole (yearly position)";
footnote 
 link="https://www.ngdc.noaa.gov/geomag/GeomagneticPoles.shtml"
 c=gray h=12pt "Data source: https://www.ngdc.noaa.gov/geomag/GeomagneticPoles.shtml";

proc print data=table_data label noobs
 style(data)={font_size=11pt}
 style(header)={font_size=11pt};
var year lat long yearly_distance;
run;

quit;
ODS HTML CLOSE;
ODS LISTING;