diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp
index c0623bd9c5bdb3ab33cc73f58f6879f117068b53..1ccb616f93e3ec8ee063513b946d5d9ffc32596e 100644
--- a/GeoLib/Raster.cpp
+++ b/GeoLib/Raster.cpp
@@ -105,6 +105,27 @@ Raster* Raster::getRasterFromSurface(Surface const& sfc, double cell_size, doubl
 	return new Raster(n_cols, n_rows, ll[0], ll[1], cell_size, z_vals, z_vals+n_cols*n_rows ,-9999);
 }
 
+double Raster::getValueAtPoint(const GeoLib::Point &pnt)
+{
+	const double* coords (pnt.getCoords());
+
+	if (coords[0]>=_ll_pnt[0] && coords[0]<(_ll_pnt[0]+(_cell_size*_n_cols)) && 
+		coords[1]>=_ll_pnt[1] && coords[1]<(_ll_pnt[1]+(_cell_size*_n_rows)))
+	{
+		int cell_x = static_cast<int>(floor((coords[0] - _ll_pnt[0])/_cell_size));
+		int cell_y = static_cast<int>(floor((coords[1] - _ll_pnt[1])/_cell_size));
+
+		// if node outside of raster use raster boundary values
+		cell_x = (cell_x < 0) ?  0 : ((cell_x > static_cast<int>(_n_cols)) ? (_n_cols-1) : cell_x);
+		cell_y = (cell_y < 0) ?  0 : ((cell_y > static_cast<int>(_n_rows)) ? (_n_rows-1) : cell_y);
+
+		size_t index = cell_y*_n_cols+cell_x;
+		if (fabs(_raster_data[index] - _no_data_val) > std::numeric_limits<float>::epsilon())
+			return _raster_data[index];
+	}
+	return _no_data_val;
+}
+
 void Raster::writeRasterAsASC(std::ostream &os) const
 {
 	// write header
diff --git a/GeoLib/Raster.h b/GeoLib/Raster.h
index 6e4c7acad9c30d5cb1a1c7964848f265b5616504..2cf4c3f440358c664b0cac354d382cfe5e734994 100644
--- a/GeoLib/Raster.h
+++ b/GeoLib/Raster.h
@@ -96,6 +96,11 @@ public:
 	 */
 	const_iterator end() const { return _raster_data + _n_rows*_n_cols; }
 
+	/**
+	 * Returns the raster value at the position of the given point.
+	 */
+	double getValueAtPoint(const GeoLib::Point &pnt);
+
 	~Raster();
 
 	/**