Obtain terrain elevation for geographic coordinates in Europe using Pandas
There is a simple and reliable way in Python to retrieve terrain elevation, commonly referred to as altitude, for WGS84 latitude and longitude coordinates over Europe, without relying on an API or paid services.
The dataset used is European Digital Elevation Model (EU-DEM), version 1.1, provided by EU’s Copernicus Programme. It consists of Rasters in GeoTIFF format divided into 1000x1000km tiles with 25m resolution. It has an aproximaste accuracy of ~7m RMSE.
The only libraries necessary are pandas, NumPy, Pyproj and Rasterio. Since the GeoTIFF files are large, the process was designed to minimize memory consumption. Tiles or chunks of the rasters that contain input points are read sequentially from storage only once. Simple and efficient functions were created in Numpy to generate the bounding box coordinates of a regular 2d grid, and to intersect points with the grid through an index. This means that spatial indexes, like a R-tree, or handling geometries with Shapely, Pygeos or Geopandas is not required.
The process has only three steps:
- Download and unzip the raster tiles from Copernicus website. Make sure you follow the terms and conditions and the dataset license. The tiles don’t need to be consecutive, you can only download the areas you need.
- Get the script from this Github repository (direct link to Jupyter notebook).
- Instantiate the object
CopernicusDEMwith the geotiff file paths and call theget_elevationmethod over a Pandas dataframe with latitude and longitude columns.
That’s it! There are other advanced options and parameters documented in the notebook.
Performance wise takes approximately 1m30s and 3500MB of ram to retrieve elevation values for 18 million points evenly spread across two tiles.
It’s also possible to create a distributed version of this script with Dask for Big Data use cases. But that’s another article.
