Mapping

% Most recent life expectancy for every country (World Bank SP.DYN.LE00.IN), as a choropleth
d = webread('https://api.worldbank.org/v2/country/all/indicator/SP.DYN.LE00.IN?format=json&mrnev=1&per_page=400');
rows = d{2};
codes = {}; vals = [];
for k = 1:numel(rows)
  r = rows{k};
  if ~isempty(r.value) && numel(r.countryiso3code) == 3
    codes{end+1} = r.countryiso3code; vals(end+1) = r.value;
  end
end
worldmap('World')
countrymap(codes, vals)
colormap viridis
title('Life expectancy at birth (years), latest year available')
subtitle('Source: World Bank. Grey: no data')

% Tonnes of CO2 per person, latest year (World Bank EN.GHG.CO2.PC.CE.AR5), zoomed in on Europe
d = webread('https://api.worldbank.org/v2/country/all/indicator/EN.GHG.CO2.PC.CE.AR5?format=json&mrnev=1&per_page=400');
rows = d{2}; codes = {}; vals = [];
for k = 1:numel(rows)
  if ~isempty(rows{k}.value) && numel(rows{k}.countryiso3code) == 3
    codes{end+1} = rows{k}.countryiso3code; vals(end+1) = rows{k}.value;
  end
end
countrymap(codes, vals, 'Europe')
colormap(hot), clim([0 15])
title('CO₂ per person (tonnes), latest year')

% All magnitude 4.5+ earthquakes in the last 30 days, from the US Geological Survey feed
q = webread('https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.geojson');
n = numel(q.features); lat = zeros(n, 1); lon = lat; mag = lat; depth = lat;
for k = 1:n
  f = q.features{k};
  c = f.geometry.coordinates;
  lon(k) = c(1); lat(k) = c(2); depth(k) = c(3); mag(k) = f.properties.mag;
end
worldmap('World')
geoscatter(lat, lon, 4*(mag - 4).^3, depth, 'filled')
colormap(turbo), clim([0 300])
title(sprintf('%d earthquakes of magnitude 4.5+ in the past 30 days', n))
subtitle('Size: magnitude · colour: depth (km) · source: USGS')

% Shortest routes on a sphere are great circles; here interpolated between the end points
names = {'New York', 'Tokyo', 'Cape Town', 'São Paulo', 'Sydney', 'Nairobi'};
lat = [40.64 35.55 -33.97 -23.43 -33.94 -1.32];
lon = [-73.78 139.78 18.60 -46.47 151.18 36.93];
a = [52.31 4.76];                                   % Schiphol
worldmap('World')
for k = 1:numel(lat)
  p1 = [cosd(a(1))*cosd(a(2)), cosd(a(1))*sind(a(2)), sind(a(1))];
  p2 = [cosd(lat(k))*cosd(lon(k)), cosd(lat(k))*sind(lon(k)), sind(lat(k))];
  w = acos(dot(p1, p2)); s = linspace(0, 1, 80)';
  P = (sin((1 - s)*w)*p1 + sin(s*w)*p2) / sin(w);   % slerp along the great circle
  geoplot(asind(P(:, 3)), atan2d(P(:, 2), P(:, 1)))
end
geoscatter([a(1) lat], [a(2) lon], 40, 'k')
km = deg2km(distance(a(1), a(2), lat, lon));
title(sprintf('Longest: %s, %.0f km', names{km == max(km)}, max(km)))

animate(12, 15)                      % 12 seconds, 15 frames per second
worldmap('globe', [30*t - 60, 25])     % centre longitude moves with t
geoscatter([52.37 40.71 35.68 -33.87 -1.29], [4.90 -74.01 139.69 151.21 36.82], 50, 'r')
title('The world turns')