Who needs GIS when you can build eye-catching 3D topography maps with Python?
Cartography twitter is currently awash with some eye catching topography maps of individual countries or regions. This is primarily due to expansion in readily available data sources and a lot of bored geography nerds who are keeping themselves busy during a pandemic. There are lots of tools and methods that can be used to generate the beautiful map of Italy shared above however in this article I am going to walk you through an unconventional approach using Python and hopefully you will come away convinced as I am that if something is worth doing, it is worth doing in Python.
There are a number of data sources that are available, all with varying degrees of resolution. On a purely aesthetic level, the higher the resolution of the image, the better the final image will look. From an accuracy point of you, lower resolution datasets can also exaggerate features and make relatively flat regions look quite hilly.
While they may not seem all that different, the image on the left uses much higher resolution data (1) compared to the image on the right (2) and hence areas that are relatively flat, for example the Po valley in northern Italy appears much hillier than it actually is in the image on the right. So going forward we will be using the higher resolution data which comes from the USGS (1). The data is free to use and further information on the data usage can be found at the following link. The data is split into 30 x 20 degree tif files which cover different parts of the globe, luckily Italy sits within one of those so we don’t need to worry about combining different tif files, although I will explore such a problem in a future article. It is very simple to download the data, simply click on the tile(s) you want and follow the download instructions. Then copy the data to wherever you like to do your programming.
First thing to do is open and read the data using rasterio, this is a relatively trivial process but I have included the code anyway. Then it is important to plot the data to get a sense for what it actually looks like and what you are dealing with. There are various functions out there to summarize data however I think they are rarely a good substitute for a plot. Some things to note, rasterio.open creates a rasterio dataset object, which contains a 2D array (Latitude x Longitude) of elevation values as well as information on the projection used and extent of the image, this will be important later. The read method reads the data in the rasterio dataset object into a 2D numpy array of (in this case) elevation values. The 2D array is used for plotting while the rasterio dataset object will be used later.
There is an immediate issue that has to be overcome, the dataset contains information for most of Southern Europe and some of North Africa when in reality, all we want is Italy. Luckily, rasterio has a useful method to clip rasters based upon a georeferenced shape, e.g. a polygon. In this case I have used the Italy polygon which is part of the NaturalEarth (3) dataset to clip the raster using the rasterio.mask.mask function. I loaded the NaturalEarth shapefile and extracted the Italy geometry using geopandas and then used it with the mask function. The mask function takes rasterio dataset object and returns the parts of the raster which are inside the polygon provided.
Now we have elevation values corresponding solely to those in Italy however there is still an issue. By default rasterio.mask.mask will fill all values which are not within the Italy polygon with 0. While this is sensible these 0’s will make the plotting tricky because they act as an anchor at the bottom of the colourmap and if there is a large gap between 0 and the minimum elevation value in the Italy data then you get the map seen above, with real data in one half of the colourmap and the other half of the colourmap absent because there is a gap between 0 and the minimum of the real data. For example, if the gap between the 0’s and the minimum value is equal to the gap between the minimum and maximum altitude then only half of the colourmap will actually be used on real data.
To get around this the mask function allows you explicitly set what value will be applied to values not within the Italy polygon with the nodatakeyword. I have included a function below which solves this problem. In the function below I pass the Italy GeoDataFrame and the rasterio dataset object. You’ll notice that the mask function is called twice, first it is called as above and values not within the Italy polygon are returned as 0’s. On the second occasion, the nodata argument is used and the values not part of the Italy polygon are set to be 1 greater than the maximum value in the Italy topography dataset (as calculated with the first mask). The result is that we now have a dataset with no natural gaps in values and the plot is starting to take shape. Also being returned by this function is the value_range variable. This corresponds to the gap between the smallest and largest value in the array and is needed when constructing colourmaps later.
Now we need to construct an appropriate colourmap. To be absolutely clear, the purpose of this article is to show you how to generate interesting but crucially, beautiful topography maps. If you are planning on creating something that military units are going to use in battle then I would suggest you use a more quantitative colourmap. In this example though we are going to build a colourmap based on the Italian flag. I have constructed one below using the colours in the Italian flag. Just using the three colours produces a colourmap with too much white in the middle so I have added an extra green and red at positions 2 and 4 in this colourmap to minimize the dominance of the white colour.
We still need to deal with the values that are not part of the Italy elevation dataset created earlier. We set those to be 1 greater than the maximum value in the elevation data. The solution is to build a colourmap with enough colours so that each unique value within the Italian elevation data has it’s own colour and then replace the last colour (red side) in the colourmap with our background colour. For example, consider a scenario where the minimum elevation value is 10, the maximum elevation value is 100 and we set the non Italian values as 101. If we create a colourmap with 91 colours, replace the 91st colour with our background colour, those non Italian values will be mapped to the 91st colour, which is our background colour.
Now, armed with our new colourmap and our clipped data we can plot the data.
While I think this still looks pretty good, it is still 2D and topography is 3D. So the final thing to do is add hillshade to mimic light shining on the topography. A hillshade is a 3D representation of a surface and are generally rendered in greyscale. The darker and lighter colors represent the shadows and highlights that you would visually expect to see in a terrain model. Hillshades are often used as an underlay in a map, to make the data appear more 3-Dimensional and thus visually interesting. We will use the earthpy hillshade function to generate our hillshade data. There are two parameters that can be tuned and they give significantly different results depending on your dataset. The first is the azimuth value which can range from 0–360 degrees and relates to where the light source is shining from. 0 degrees corresponds to a light source pointing due north. The second is the altitude that the light source is at and these values can range from 1–90. Below are a few examples highlighting how changing the two values can give vastly different results.
In the end I decided on an azimuth value of 180 so that I got a nice shadow on the southern side of the Alps and an altitude value of 1. I would encourage you to play around with these values yourself though. Finally it is time to plot the finished product. The main Italy topography data is plotted first and the hillshade is plotted on top with a small alpha value.
Zooming in allows you to really see the detail of the image!
And there we have it, a very visually pleasing topography map of Italy, ready to be framed and hung on your wall. This method can be applied to any country or region however additional steps will be required if the data for a country is spanning multiple tif files. I will be covering this problem in a future post so subscribe to ensure you see it.
(1) —High resolution data source – Danielson, J.J., and Gesch, D.B., 2011, Global multi-resolution terrain elevation data 2010 (GMTED2010): U.S. Geological Survey Open-File Report 2011–1073, 26 p. doi:10.5066/F7J38R2N Data is free to use, further information on the use can be found here –https://topotools.cr.usgs.gov/GMTED_viewer/gmted2010_fgdc_metadata.html
(2) —Low resolution data source – doi:10.7289/V5C8276M
(3) — NaturalEarth — https://www.naturalearthdata.com/