Looking for altitude OpenData

Recently, for a personal project, I looked for opendata for altitude in France. These are not so easy to find, although there are some really good opendata out there, with a very great accuracy. Indeed, OpenStreetMaps provides really accurate opendata for maps (street names and so on), but absolutely nothing concerning the altitude.

I finally found two datasets available for France, one published by IGN and the other one being the SRTM dataset published by the NASA. The one from IGN provides altitude information for the whole French territory, at a resolution of 75m, which is quite good for numerically modeling the terrain. The one from NASA is more complete and provides altitude information for the whole globe surface between -60° and 60° latitude at a resolution of about 90m, obtained by averaging on 30m samples.

IGN BD ALTI opendata

First of all, let’s have a look at the BD ALTI from IGN, which is not as widely used as the SRTM one (at least in a hobbyist context) and does not have much documentation around there.

IGN provides an archive containing the data for the whole French territory. In the archive, there are several folders. The interesting one, in our case, is 1_DONNEES_LIVRAISON_* which contains all the data files.

Altitude data is available in two formats: GeoTiff (Tiff with geolocation metadata) and ASCII Esri grid. I personally used the ASCII Esri Grid files. For more info about this file format, you can refer to the Wikipedia page.

Each file starts with some header informations: number of rows and cols (1000), position of the lower-left corner in X and Y coordinates, size of an elementary cell and value used to indicate missing data. The file then contains a matrix representing the values of each cell, as described in the Wikipedia page.

Position on the map is indicated in X and Y coordinates, and not using the standard latitude and longitude coordinates. This is a different projection system, called Lambert93, and used for official maps of France. Wikipedia contains the relevant information to translate from the standard WGS84 (satellite latitude and longitude coordinates) to the Lambert93 coordinates, and reciprocally but the best way seems to be to use the Proj.4 libray which has a Python binding to handle the conversion.

To find a specific point, you should then iterate over the available data files and check if the position of the point is between xllcorner, yllcorner and xllcorner + ncols * cellsize, yllcorner + nrows * cellsize. Then, you just have to get the value of the matching cell in the matrix description in the file.

An example script is available at https://github.com/Phyks/IGN_OpenData. It takes a rectangle defined by two latitudes and longitudes in input, as well as the data folder to process, and will output the relevant data from the BD ALTI dataset, in XYZ file format, which is equivalent to a CSV file with columns X, Y and altitude (Z), separated by tabulations, and sorted. This is also an example use of the Proj.4 Python binding to handle basic coordinates system transforms.

SRTM opendata

The SRTM dataset exists in multiple versions and different samplings. Have a look at this wiki page to have more information about it. I was interested in SRTMGL3 version 3, which is the third iteration of the dataset sampled at 3 arcseconds (which is about 90m). It is available at http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL3.003/ (before clicking the main folder in this address, please note that this is very heavy and long to process as it contains all the files for the whole world). rv uses it and provides a basic import script to a PostgreSQL database, available here.

All files are named with a convention to describe the tile they represent. For instance, the tile starting at 0° North and 6° East is called N00E006.*.

There is a really nice StackOverflow post to explain how to handle these files and the following will almost be a rewriting of this post.

First of all, take the integer part of the position you want (latitude and longitude), and look for the corresponding file at the previous address. You will download a (zipped) hgt file containing the altitude data. This file is a big-endian binary file that you will have to parse.

First, you have to find the pixels coordinates representing your geographic point. To this purpose, you have to take the floating part of the latitude and longitude of the point you want and convert them to arcseconds (e.g. for 50° 24’ 58.888”, you will have 24 * 60 + 58.888 = 1498.888). As the sampling is done at 3 arcseconds, the pixel coordinates are given by the closest integer to the latitude (resp. longitude) of your point, expressed in arcseconds, divided by 3. Latitude will give you the row index, whereas longitude will give you the col index.

The first row in the file is the northernmost one, so if we are looking at row 500 from the lower edge, we actually have to look at row 1201 - 500 from the beginning of the files (1200 being 3600 / 3, the total number of rows in the data file).

We then have to iterate on cols to find the correct one, using the same principle.

An example script is available here although it was not thoroughly tested. If you fix any problem with it, please report :)