What is the roundest country?

We can define roundness in many ways. For example, as you may know, the circle is the shape that given a fixed perimeter maximizes the area. This definition has many problems. One of the problems is that countries generally have chaotic perimeters (also known as borders), so they tend to be much longer than they seem to be.

For that reason, we have to define roundness some other way. Given a country, I will represent it as a plane region, more precisely a compact set \(C \subset \mathbb{R}^2\), and I will define its roundness as

\[ roundness(C) = \max_{ x \in \mathbb{R}^2, r \in \mathbb{R}_{>0}
} \frac{ area(C \cap D(x,r)) }{ \max \{ area(D(x,r)), area(C) \} } \]

where \(D(x, r)\) is the disk of center \(x\) and radius \(r\).

Note that:

• For all \(C\), \(0 \leq roundness(C) \leq 1 \).
• If \(C\) is a circle, then taking \(D(x,r) = C\) maximizes the formula, so \(roundness(C) = 1\).
• If \(roundness(C) = 1\), then \(C\) is a circle, because it overlaps completely with \(D(x,r)\), for some \(x\) and \(r\).

Let’s see the results and then we will see the code I used to get them. You can scroll to see all the results, and sort by any column you like.

Dataset

First off, we need something that tells us the shape of each country. For that, I used the Natural Earth sovereign states dataset. It contains a .shp file, that for each country stores a sequence of points \(\{(x_i, y_i)\}_{0 \leq i < N} \) and you have to connect the dots between them to get the shape of the country.

It also has information about the parts of each country (the \(C\) we talked about before isn’t necessarily connected), it tells us which \((x_i, y_i)\) belong to each part.

Then, a country will be represented as a finite union of polygons, one polygon per part of the country.

To read the .shp file I used the pyshp library.

There’s a small problem too that I have to talk about. We cannot use those \((x_i, y_i)\) directly to draw the countries in a plane because they are in WGS84 format (a very accurate longitude-latitude system), so we would be doing an equirectangular projection of the Earth, which really distorts a lot of countries.

So what, for each country, I transform its points from the equirectangular projection to an azimutal projection centered in some point that depends on the country in question.

That point is obtained by computing the midpoint of two random points from the border of the country. This is not the best (if the country is not convex, then the midpoint isn’t necessarily inside it), but it works pretty well.

In order to do all the transformations I used the pyproj library.

Algorithm

For each country, what we have to do is to find the circle that maximizes the formula at the top of the article. In order to do that, I use a regular gradient descent algorithm. I have two initial guesses and make them a little better each iteration.

But let’s talk about inscribed and circumscribed circles first. Given a geometrical shape (for example a rectangle), there is a circle that fits inside (resp. outside) and maximizes (resp. minimizes) its area. So, for each country we compute a frame that contains it (a rectangle) and compute the inscribed and circumscribed circles for that rectangle. Those two circles will be the initial guesses.

As we said, we apply a gradient descent algorithm to those circles. In order to use gradient descent, we have to define the neighbourhood of a circle, that is the circles to which we will try to move to make our guess better.

The neighbourhood of a circle will be the circles that we get by moving it a little up, down, left or right, and also making it a little bigger or smaller.

As always, the gradient method has the problem of local maxima, because it can get stuck in them. I don’t care about that, given that our dataset is prone to having only one maxima.

I also run the gradient descent algorithm with a couple more of guesses, that come from the inscribed and circumscribed circles of the frames of the largest parts of the country.

To compute areas of polygons, and intersection between polygons and circles, I used the shapely library, that is impressively good.

Code

If you want to see the details of the implementation, feel free to take a look at the code. It outputs code to generate the table in this post, it is easy to change it to get the .png files.

Epilogue

The idea of this post came from a great post that analyzes the rectangularness of countries. We both use gradient descent, to find an approximated solution to our problems, although the code is completely different. We used the same dataset though, so I should make the same clarification: all geopolitical decisions were outsourced to Natual Earth.

Edit (28 Jul 2016): Updated list and included a fix for projections. Now I use an azimuthal projection centered in each country, which allows a more precise representation. All the details can be seen in the gist, which has been updated accordingly. Also added an index for polygons to make intersection computation faster.

◑