I've read several articles that mentioned the north magnetic pole has been moving more in the past few decades, than in the previous few hundred years. And as a Map Guy, I knew I just had to plot this data on a map, and see it for myself! I provide a lot of details below - but if you're not interested in the details, feel free to jump to the final map at the end!

Magnetic North, Maps, and Compasses

Before we start on the map, here's a bit of trivia related to magnetic north. When map-makers draw a map that that will be used for something important (such as navigation), they must indicate the magnetic declination on the map, and that map is only good for a short time (because the magnetic north pole moves). Therefore, to use a map for accurate navigation, it must be a recent map, and the map must indicate the magnetic declination. And then you must account for that magnetic declination in your compass readings (some compasses let you set the magnetic declination, so you don't have to add offsets manually).

My friend Gene-Paul is a forest firefighter, and he sent me this picture of some of the compasses they use. Who knew there were this many kinds of compasses!?! (Thanks Gene-Paul!)

Existing Maps

When you're thinking about creating your own map, it's always best to first check what's already out there. You might find some good things you want to include in your map, and some bad things you might want to avoid in your map. Here's a map of the shifting magnetic north pole, from a NOAA page:

I like that they have used a map that shows the terrain (both above, and below water), but I don't like that the text labels are too small for me to read (text labels for the colored markers, and the longitude lines). Also, if I didn't already know what area the magnetic north pole was in, I would not know that the land area shown in this map was Canada. This gave me a few things to keep in mind when creating my map.

The Data

To create my own map, I first needed the data. I did a few web searches, and soon found a good data source. The National Oceanic and Atmospheric Administration (NOAA) has a web page about the Wandering of the Geomagnetic poles (that's where I got the map above), and at the bottom of that page they have a link to the annual International Geomagnetic Reference Field (IGRF) model data. Here's what the text data looks like - the columns are longitude, latitude, and year.

I downloaded/saved their data file (NP.xy.txt) from the browser window, and imported it into SAS with the following code:

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;

My Preliminary SAS Maps

One of the great things about the new Proc SGmap is that it makes plotting lat/long data as markers on a tile map really easy, with a minimum of code. But when I plot this particular data on an OpenStreetmap tile map, it just doesn't produce a very useful map. The default OpenStreetmaps are not really meant for plotting data that's this close to the north pole.

proc sgmap plotdata=magnetic_pole_locations noautolegend;

openstreetmap;

scatter x=long y=lat / markerattrs=(color=red symbol=circlefilled size=5px);

run;

Therefore, instead of using an OpenStreetmap (above), let's use a special Esri polar map that's specifically designed for plotting data near the north pole. With this change (below), it's looking a lot more like the NOAA map from the original article, and is starting to become useful.

proc sgmap plotdata=magnetic_pole_locations noautolegend;

esrimap url='http://services.arcgisonline.com/arcgis/rest/services/Polar/Arctic_Ocean_Base/MapServer';

scatter x=long y=lat / markerattrs=(color=red symbol=circlefilled size=5px);

run;

Next I want to add some lines (or circles) of latitude. This will help see how close to the true north pole the points are, and will provide good reference lines to help get your brain oriented to the map. I use a data step to loop through some values along the rings, and save them in a dataset. Notice that rather than using variables named lat and long, I use lat_ring and long_ring. By using different variable names, this allows me to combine this dataset with my magnetic_pole_locations dataset (you can only specify one plotdata= dataset in SGmap), and still plot the rings and markers separately.

data rings;

do lat_ring = 70 to 90 by 5;

do long_ring = -180 to 180 by 5;

output;

end;

end;

run;

data my_data; set magnetic_pole_locations rings;

run;

proc sgmap plotdata=my_data noautolegend;

esrimap url='http://services.arcgisonline.com/arcgis/rest/services/Polar/Arctic_Ocean_Base/MapServer';

series x=long_ring y=lat_ring / group=lat_ring lineattrs=(color=gray55) tip=none;

scatter x=long y=lat / markerattrs=(color=red symbol=circlefilled size=5px);

run;

Notice that since Proc SGmap automatically zooms in/out to select a map area that 'fits' the data, plotting these rings zooms that map out a bit, and centers it on the north pole. When studying the location of magnetic north, I think it's good to show the map centered on the north pole like this.

I like having the rings of latitude ... but someone looking at the map might not know exactly which rings these are. Therefore I add a text label at longitude 120 along each ring. Notice that I define an escape character, so I can add the '00ba'x character (degrees symbol) to the end of the latitude number.

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;

proc sgmap plotdata=my_data noautolegend;

esrimap url='http://services.arcgisonline.com/arcgis/rest/services/Polar/Arctic_Ocean_Base/MapServer';

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);

run;

Now that we've got text labels on the latitude rings, let's also add some text labels to the red markers. The easiest and most direct way to add labels would probably be using the datalabel= option on the scatter statement. But I wanted to exercise a lot of control over the labels, and (in my opinion) the easiest and most direct way to do that is by adding another text statement (or statements, in this case). For the first few hundred years, where the data points are very close together (because the magnetic pole didn't move much), I want to label every 50 years with black text. And for the more recent years, I want to label every 10 years with blue text.

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;

proc sgmap plotdata=my_data noautolegend;

esrimap url='http://services.arcgisonline.com/arcgis/rest/services/Polar/Arctic_Ocean_Base/MapServer';

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);

text x=long y=lat text=year_label / textattrs=(color=black size=9pt) position=left;

text x=long y=lat text=recent_label / textattrs=(color=blue size=9pt) position=left;

run;

My Final SAS Map

If you're map-savvy, you might recognize Greenland in the map above, and from Greenland you could intuit where Canada is, etc. But it sure would be easier if we had some labels for the land areas. Rather than relying on an automated technique to calculate the lat/long location for the labels (such as using the %centroid macro), I decided to hard-code the values. This way I can place the labels exactly where I want them.

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;

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";

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);

text x=long y=lat text=year_label / textattrs=(color=black size=9pt) position=left;

text x=long y=lat text=recent_label / textattrs=(color=blue size=9pt) position=left;

run;

Extra Credit

One big question about this data is "how far is the magnetic pole moving each year?" I used the SAS geodist() function to calculate the distance (in miles) between each year's lat/long location, and the previous year's lat/long location.

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;

But how do I show this distance to the user? There's definitely not enough room to place labels on the map for each dot. Therefore I used the tip= option to add custom mouse-over text for each red marker. Note that this is a very new SGmap feature, and you'll need the latest release (Viya 3.5) to use it.

scatter x=long y=lat / markerattrs=(color=red symbol=circlefilled size=5px)

tip=(year lat long yearly_distance);

Here's a screen-capture of what the mouse-over text looks like:

It would be quite cumbersome to mouse-over each of the red markers to see all the distances, or to find the marker for a particular year. Therefore I also created a text table of values under the map (using Proc Print). Here's a portion of the table:

Discussion

Did you like this map? Do you have any of your own projects where you could re-use some of the tricks I've shown here? What other analyses might you use on this type of data? What are some of the factors causing the magnetic north pole to shift? What's your forecast - will the magnetic poles shift more, or less, in the next 100 years? (Feel free to discuss in the comments section.)

Here's a link to the interactive version of the map, if you'd like to see the mouse-over text and full table. And here's a link to the SAS code, if you'd like to see all the details.