Choropleths

A choropleth is a thematic map where areas are colored in proportion to a variable such as population density. GMT lets us plot choropleth maps but the process is not straightforward because we need to put the color information in the header of a multi-segment file. To facilitate this, some tools were added to the Julia wrapper.

(But see also the more extended example in tutorials)

Countries (Europe DCW)

Load packages needed to download data and put into a DataFrame

using GMT, DataFrames, CSV

paises = CSV.File(HTTP.get("https://raw.githubusercontent.com/tillnagel/unfolding/master/data/data/countries-population-density.csv").body, delim=';') |> DataFrame;

Extract the Europe countries from the DCW file

eu = coast(DCW="=EU+z", dump=true);

Generate two vector with country codes (2 chars codes) and the population density for European countries.

codes, vals = GMT.mk_codes_values(paises[!, 2], paises[!, 3], region="eu");

Create a Categorical CPT that allows plotting the choropleth. Note that we need to limit the CPT range to leave out the very high density states like Monaco otherwise all others would get the same color.

C = cpt4dcw(codes, vals, range=[10, 500, 10]);

Show the result

plot(eu, cmap=C, fill="+z", proj=:guess, R="-76/36/26/84", title="Population density",
	colorbar=true, show=true)

Countries (World)

While GMT come with the world countries in its DCW data base, those countries polygons have a somewhat too high resolution for world maps and hence produce files unnecessarily, which also makes them slower to create. An alternative is to use a much lower resolution file that can be found here but whose download link is here. And we can read it directly into GMT (and wait a bit while it gets downloaded).

countries = gmtread("/vsicurl/https://github.com/datasets/geo-countries/raw/master/data/countries.geojson");

The country polygons have attributes like:

Like in the Europe case above, we load the world population into a DataFrame

pop = CSV.File(HTTP.get("https://raw.githubusercontent.com/tillnagel/unfolding/master/data/data/countries-population-density.csv").body, delim=';') |> DataFrame;

The contents of that file look like this:

julia> pop
248×3 DataFrame
 Row │ Country Name                       Country Code  2010
     │ String                             String3       Float64?
─────┼───────────────────────────────────────────────────────────
   1 │ Arab World                         ARB            26.2725
   2 │ Caribbean small states             CSS            16.9939
   3 │ East Asia & Pacific (all income …  EAS            90.5122

So we must now link the countries in countries with density in 2010 via the attribute ISO_A3 in countries and the third column in pop. To make a global choropleth map with the population density we will use first the polygonlevels function to get as the population density in the same order as the country polygons are stored in the countries GMTdataset.

zvals = polygonlevels(countries, string.(pop[!,2]), pop[!,3], att="ISO_A3");

And now we are ready to make the map

# First create the colormap. Limit the maximum range to 500 otherwise states
# like Monaco would the all the `reds` and the rest of the world would be all blues.
C = makecpt(range=(0,500,10));
plot(countries, region=(-180,180,-60,85), level=zvals, cmap=C, proj=:guess,
     title="World population density", colorbar=true, show=true)

Covid rate of infection in Portugal 2021

First, download the Portuguese district polygons shape file from this Github repo Next load it with:

using GMT, DataFrames, CSV

Pt = gmtread("C:\\programs\\compa_libs\\covid19pt\\extra\\mapas\\concelhos\\concelhos.shp");

Download and load a CSV file from same repo with rate of infection per district. Load it into a DataFrame to simplify data extraction.

incidence = CSV.read("C:\\programs\\compa_libs\\covid19pt\\data_concelhos_incidencia.csv", DataFrame);

Get the rate of incidence in number of infected per 100_000 habitants for the last reported week.

r = collect(incidence[end, 2:end]);

But the damn polygon names above are all uppercase, Ghrrr. We will have to take care of that.

ids = names(incidence)[2:end];

Each of the Pt datasets have attributes (e.g., Pt[1].attrib) and the one that is common with the names in ids is the Pt[1].attrib["NAME2"] (the conselho name). But the names in dataconcelhosincidencia.csv (from which the ids are derived) and the concelhos.shp (that we read into Pt) do not use the same case (one is full upper case) so we need to use the nocase=true below. The comparison is made inside the next call to the polygonlevels() function that takes care to return the numerical vector that we need in plot's level option.

zvals = polygonlevels(Pt, ids, r, att="NAME_2", nocase=true);

Create a Colormap to paint the district polygons

C = makecpt(range=(0,1500,10), inverse=true, cmap=:oleron, hinge=240, bg=:o);

Get the date for the data being represented to use in title

date = incidence[end,1];

And finaly do the plot

plot(Pt, level=zvals, cmap=C, pen=0.5, region=(-9.75,-5.9,36.9,42.1), proj=:Mercator,
     title="Infected / 100.000 habitants " * date)

colorbar!(pos=(anchor=:MR,length=(12,0.6), offset=(-2.4,-4)), color=C,
          axes=(annot=100,), show=true)