Arctic sea ice map with MapLibre GL JS

The ice that covers the Arctic Ocean shrinks every summer and is smallest in September. Satellites have measured it since 1979, and this map plays every September since then over the ice of 1979, so the orange is the ice that's gone. The ice has pulled back from the coasts of Alaska and Siberia, and in a few summers it opened the sea routes along them. Slide to any year.

Sea ice that year Sea ice in 1979 that's gone

Build it yourself

A Node.js script turns NSIDC's shapefiles into one GeoJSON file, with a MultiPolygon of ice for each September, and a page draws the year you pick over 1979 on a dark globe. It's the map above without the play button. It was tested with MapLibre GL JS 6.15 in Chromium and Node.js 22.

What you need

  • Node.js 18 or later.
  • A browser key with the Map tiles product, made as API keys shows. Add http://localhost:* to its allowed origins, so it works from your computer while you build. NSIDC's data needs no key.

Collect the ice

Make a folder for the map, and in it a file called sea-ice.mjs:

sea-ice.mjs
import { writeFile } from "node:fs/promises";

const base = "https://noaadata.apps.nsidc.org/NOAA/G02135/north/monthly";

// NSIDC's polar stereographic projection, turned back into longitude and latitude.
const a = 6378273;
const f = 1 / 298.279411123064;
const e = Math.sqrt(2 * f - f * f);
const sin70 = Math.sin((70 * Math.PI) / 180);
const tc = Math.tan((10 * Math.PI) / 180) / ((1 - e * sin70) / (1 + e * sin70)) ** (e / 2);
const mc = Math.cos((70 * Math.PI) / 180) / Math.sqrt(1 - (e * sin70) ** 2);
const k = tc / (a * mc);
function lonLat([x, y]) {
  const t = Math.hypot(x, y) * k;
  let lat = Math.PI / 2 - 2 * Math.atan(t);
  for (let i = 0; i < 5; i++) lat = Math.PI / 2 - 2 * Math.atan(t * ((1 - e * Math.sin(lat)) / (1 + e * Math.sin(lat))) ** (e / 2));
  const lon = (Math.atan2(x, -y) * 180) / Math.PI - 45;
  return [lon < -180 ? lon + 360 : lon, (lat * 180) / Math.PI];
}

const area = (ring) => ring.reduce((sum, [x, y], i) => sum + x * ring[(i + 1) % ring.length][1] - ring[(i + 1) % ring.length][0] * y, 0) / 2;
function inside([x, y], ring) {
  let found = false;
  for (let i = 0, j = ring.length - 1; i < ring.length; j = i++) {
    const [[xi, yi], [xj, yj]] = [ring[i], ring[j]];
    if (yi > y !== yj > y && x < ((xj - xi) * (y - yi)) / (yj - yi) + xi) found = !found;
  }
  return found;
}

function geographic(ring, hole) {
  const points = [];
  for (const point of ring.map(lonLat)) {
    const previous = points.at(-1);
    if (previous) point[0] = previous[0] + ((((point[0] - previous[0]) % 360) + 540) % 360) - 180;
    points.push(point);
  }
  if (inside([0, 0], ring)) points.push([points.at(-1)[0], 90], [points[0][0], 90], points[0]);
  const shift = 360 * Math.round(points.reduce((sum, [lon]) => sum + lon, 0) / points.length / 360);
  const result = points.map(([lon, lat]) => [Math.round((lon - shift) * 100) / 100, Math.round(lat * 100) / 100]);
  return area(result) > 0 === hole ? result.reverse() : result;
}

// Each ring of the shapefile in a zip that NSIDC stores without compressing.
function rings(zip) {
  let shp;
  for (let i = 0; zip.readUInt32LE(i) === 0x04034b50; ) {
    const size = zip.readUInt32LE(i + 18);
    const nameLength = zip.readUInt16LE(i + 26);
    const start = i + 30 + nameLength + zip.readUInt16LE(i + 28);
    if (zip.toString("utf8", i + 30, i + 30 + nameLength).endsWith(".shp")) shp = zip.subarray(start, start + size);
    i = start + size;
  }
  const found = [];
  for (let i = 100; i < shp.length; i += 8 + shp.readInt32BE(i + 4) * 2) {
    const record = i + 8;
    if (shp.readInt32LE(record) !== 5) continue;
    const parts = shp.readInt32LE(record + 36);
    const count = shp.readInt32LE(record + 40);
    const points = record + 44 + parts * 4;
    for (let part = 0; part < parts; part++) {
      const end = part + 1 < parts ? shp.readInt32LE(record + 48 + part * 4) : count;
      const ring = [];
      for (let p = shp.readInt32LE(record + 44 + part * 4); p < end; p++) ring.push([shp.readDoubleLE(points + p * 16), shp.readDoubleLE(points + p * 16 + 8)]);
      found.push(ring);
    }
  }
  return found;
}

function polygons(found) {
  const result = found.filter((ring) => area(ring) < 0).map((ring) => [ring]);
  for (const hole of found.filter((ring) => area(ring) > 0)) {
    const [[x1, y1], [x2, y2]] = hole;
    const length = Math.hypot(x2 - x1, y2 - y1);
    const point = [(x1 + x2) / 2 + ((y2 - y1) / length) * 1000, (y1 + y2) / 2 - ((x2 - x1) / length) * 1000];
    result.find(([outside]) => inside(point, outside)).push(hole);
  }
  return result.map(([outside, ...holes]) => [geographic(outside, false), ...holes.map((hole) => geographic(hole, true))]);
}

const csv = await (await fetch(`${base}/data/N_09_extent_v4.0.csv`)).text();
const features = [];
for (const row of csv.trim().split("\n").slice(1)) {
  const [year, , , , extent] = row.split(",").map((field) => field.trim());
  if (!(Number(extent) > 0)) continue;
  const response = await fetch(`${base}/shapefiles/shp_extent/09_Sep/extent_N_${year}09_polygon_v4.0.zip`);
  if (!response.ok) throw new Error(`NSIDC returned ${response.status} for ${year}`);
  features.push({
    type: "Feature",
    properties: { year: Number(year), extent: Number(extent) },
    geometry: { type: "MultiPolygon", coordinates: polygons(rings(Buffer.from(await response.arrayBuffer()))) },
  });
}
await writeFile("sea-ice.geojson", JSON.stringify({ type: "FeatureCollection", features }));
console.log(`Saved ${features.length} Septembers, ${features[0].properties.year} to ${features.at(-1).properties.year}`);
  • NSIDC publishes a zipped shapefile for every month since late 1978, and a CSV of each September's extent: the area of ocean where at least 15% of the surface is ice.
  • The shapefiles are in NSIDC's own polar stereographic projection, in meters from the North Pole, so lonLat turns each point back into longitude and latitude with the projection's inverse formulas.
  • NSIDC stores the files in each zip without compressing them, so the script finds the .shp file by its header and reads its polygon records as they are, with no libraries.
  • A shapefile winds each polygon's outside clockwise and its holes the other way, and doesn't say which hole belongs to which polygon. polygons puts each hole with the polygon around a point just outside the hole's first edge, on the ice side.
  • The ice around the North Pole doesn't fit in longitude and latitude as it is, since its edge goes all the way around the world. geographic lets its longitudes run on past 180°, then closes the ring along the pole itself, at latitude 90°. MapLibre handles longitudes past 180° as long as they're within one world of the map.

Run it:

Terminal
node sea-ice.mjs

It downloads about 2 MB, then prints Saved 48 Septembers, 1979 to 2026, or however many there are when you run it, and writes sea-ice.geojson.

Make the page

Next to sea-ice.geojson, make index.html, with your browser key in place of mapricot_pk_…:

index.html
<!doctype html>
<html lang="en">
  <head>
    <meta charset="utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1">
    <title>Arctic sea ice</title>
    <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/maplibre-gl@6.15.0/dist/maplibre-gl.css">
    <script type="importmap">
      {
        "imports": {
          "maplibre-gl": "https://cdn.jsdelivr.net/npm/maplibre-gl@6.15.0/+esm"
        }
      }
    </script>
  </head>
  <body style="margin: 0">
    <label style="position: absolute; z-index: 1; top: 10px; left: 10px; padding: 4px 8px; background: white; font: 16px sans-serif">
      September <span id="year">2026</span>
      <input id="slider" type="range" min="1979" max="2026" value="2026">
    </label>
    <div id="map" style="height: 100vh"></div>
    <script type="module">
      import * as maplibregl from "maplibre-gl";

      const map = new maplibregl.Map({
        container: "map",
        style: "https://api.mapricot.com/tiles/v1/styles/dark.json?key=mapricot_pk_…",
      });
      map.on("style.load", () => {
        map.setProjection({ type: "globe" });
        map.jumpTo({ center: [-45, 84.5], zoom: -0.15 });
      });
      await map.once("load");

      const labels = map.getStyle().layers.find((layer) => layer.type === "symbol").id;
      map.addSource("ice", { type: "geojson", data: "sea-ice.geojson", attribution: "Sea ice: NSIDC Sea Ice Index" });
      map.addLayer({ id: "1979", type: "fill", source: "ice", filter: ["==", ["get", "year"], 1979], paint: { "fill-color": "#c25a1d" } }, labels);
      map.addLayer({ id: "ice", type: "fill", source: "ice", filter: ["==", ["get", "year"], 2026], paint: { "fill-color": "#ffffff" } }, labels);

      const slider = document.querySelector("#slider");
      slider.addEventListener("input", () => {
        document.querySelector("#year").textContent = slider.value;
        map.setFilter("ice", ["==", ["get", "year"], Number(slider.value)]);
      });
    </script>
  </body>
</html>
  • The globe is set once the style loads, since a style decides its own projection, and then the view moves to above the pole. A flat map can't show the pole, so it has to be a globe before the view goes there.
  • The ice goes in under the style's first labels, so place names stay on top.
  • Both layers draw the same source, filtered to one year each: 1979 in orange, then the chosen year in white on top, so the orange left showing is the ice that's gone.
  • The slider sets the white layer's filter, so it redraws without loading the file again. When a new September comes out, raise the slider's max and value to it.

Open it

Serve the folder with npx serve or python3 -m http.server, and open the address it prints. Opened straight from the file, the map stays blank, as Open the page explains. You should see the Arctic Ocean from above, with this year's ice in white, and orange where the ice was in 1979, mostly along the coasts of Alaska and Siberia. Slide back to 2012 for the least ice on record.

NSIDC has the same files for every month and for the Antarctic, where S takes the place of N in the file names. The Antarctic uses a different projection, so lonLat needs its south polar formulas.

Data

May we use an analytics cookie? It shows us how people find and use Mapricot.