diff --git a/src/data/dem.cpp b/src/data/dem.cpp
index 9ee68cfd..722042ba 100644
--- a/src/data/dem.cpp
+++ b/src/data/dem.cpp
@@ -9,6 +9,21 @@
 	((samples) * (samples) * 2)
 
 
+static qreal interpolate(qreal dx, qreal dy, qreal p0, qreal p1, qreal p2,
+  qreal p3)
+{
+	return p0 * (1 - dx) * (1 - dy) + p1 * dx * (1 - dy) + p2 * dy * (1 - dx)
+	  + p3 * dx * dy;
+}
+
+static qreal value(int col, int row, int samples, const QByteArray data)
+{
+	int pos = ((samples - 1 - row) * samples + col) * 2;
+	qint16 val = qFromBigEndian(*((const qint16*)(data.constData() + pos)));
+
+	return (val == -32768) ? NAN : val;
+}
+
 static qreal height(const Coordinates &c, const QByteArray data)
 {
 	int samples;
@@ -20,13 +35,17 @@ static qreal height(const Coordinates &c, const QByteArray data)
 	else
 		return NAN;
 
-	int lat = int(qRound((c.lat() - int(c.lat())) * (samples - 1)));
-	int lon = int(qRound((c.lon() - int(c.lon())) * (samples - 1)));
+	int row = (int)((c.lat() - (int)c.lat()) * (samples - 1));
+	int col = (int)((c.lon() - (int)c.lon()) * (samples - 1));
+	qreal dx = ((c.lon() - (int)c.lon()) * (samples - 1)) - col;
+	qreal dy = ((c.lat() - (int)c.lat()) * (samples - 1)) - row;
 
-	int pos = ((samples - 1 - lat) * samples + lon) * 2;
-	qint16 val = qFromBigEndian(*((const qint16*)(data.constData() + pos)));
+	qreal p0 = value(col, row, samples, data);
+	qreal p1 = value(col + 1, row, samples, data);
+	qreal p2 = value(col, row + 1, samples, data);
+	qreal p3 = value(col + 1, row + 1, samples, data);
 
-	return (val == -32768) ? NAN : val;
+	return interpolate(dx, dy, p0, p1, p2, p3);
 }
 
 QString DEM::fileName(const Key &key) const