[GIS] Extracting values from raster under line using GDAL

Question

How can I find the cells which lie on a straight line between two points?
I.e. how can I find the indices of the cells (in grey) which are intersected by the line defined by the two points (in blue):

I want to do this in order to read off values from the underlying DTM and thus extract the height profile along the line….

Using the geotransform of the raster, it is easy to convert the two points to cell offsets, so my first thought was to create many points along the line and convert each of these. However the problem here is creating points at the right distance to successfully capture each intersected cell (while also not creating a ridiculous number of points so that it would be incredibly inefficient)

Looking for ideas that are either programmatic (python/C++) or theory (e.g. no arcgis tools! gdal ok…)

Answer 1

You can use gdal_rasterize command as base of the procedure. I tested my approach with this situation:

The code is:

import os

#extent of base raster
extent = "-te 354971.3488602247089148 4473058.4587739501148462  355341.2924763654009439 4473428.4023900907486677 "
attribute_name = " -a id "
#resolution of base raster
resolution = " -tr 73.9887 73.9887 "
input_vector_line = " -l test_line /home/zeito/pyqgis_data/test_line.shp "
output_raster = " /home/zeito/pyqgis_data/test_line_raster.tif"

cmd = "gdal_rasterize -at " + \
                     extent + \
             attribute_name + \
                 resolution + \
          input_vector_line + \
              output_raster

os.system(cmd)

where the option -at is very important because "Enables the ALL_TOUCHED rasterization option so that all pixels touched by lines or polygons will be updated, not just those on the line render path, or whose center point is within the polygon. Defaults to disabled for normal rendering rules" (http://www.gdal.org/gdal_rasterize.html).

After running the code at the Python Console of QGIS I got:

where the line was adequately rasterized with a constant value of 1. Now, with map algebra, you can extract the height profile along the line.

You can use the GDAL options to get raster extent, resolution, etc, for automatizing all procedure.