diff --git a/docs/src/advanced.md b/docs/src/advanced.md
index 88d31d0..3da4068 100644
--- a/docs/src/advanced.md
+++ b/docs/src/advanced.md
@@ -198,7 +198,7 @@ saveas("advanced013b") # hide
 ```
 ![](assets/advanced013b.png)
 
-The [`hist()`](@ref) function compute also 2D histograms by passing two vectors (with the same lengths), e.g.: 
+The [`hist()`](@ref) function compute also 2D histograms by passing two vectors (with the same lengths), e.g.:
 ```@example abc
 x = randn(10_000)
 y = randn(10_000)
@@ -272,6 +272,53 @@ saveas("advanced014f") # hide
 ![](assets/advanced014f.png)
 
 
+## Interpolation of 2D scattered data
+The `dgrid3d()` function allows to interpolate 2D scattered data onto a 2D regular grid, e.g.:
+```@example abc
+x = (rand(200) .- 0.5) .* 3;
+y = (rand(200) .- 0.5) .* 3;
+z = exp.(-(x.^2 .+ y.^2));
+
+# Interpolate on a 20x30 regular grid with splines
+gx, gy, gz = dgrid3d(x, y, z, "20,30 splines")
+
+@gsp "set size ratio -1" "set xyplane at 0" xlab="X" ylab="Y" :-
+@gsp :-  x  y  z "w p t 'Scattered data' lc pal"
+@gsp :- gx gy gz "w l t 'Interpolation on a grid' lc pal"
+saveas("advanced015a") # hide
+```
+![](assets/advanced015a.png)
+
+
+!!! warn
+    The `splines` algorithm may be very slow on large datasets.  An alternative option is to use a smoothing kernel, such as `gauss`.
+
+
+The interpolated data in scarcely sampled regions are poorly constrained, i.e. they are actually *extrapolated values*.  By using the `extra=false` keyword all extrapolated values are set to `NaN`:
+
+```@example abc
+x = randn(2000) .* 0.5;
+y = randn(2000) .* 0.5;
+rsq = x.^2 + y.^2;
+z = exp.(-rsq) .* sin.(y) .* cos.(2 * rsq);
+
+@gsp "set size ratio -1" palette(:balance, smooth=true) "set view map" "set pm3d" :-
+@gsp :- "set multiplot layout 1,3" xr=[-2,2] yr=[-2,2] :-
+@gsp :- 1 tit="Scattered data"  x  y  z "w p notit lc pal"
+
+# Show extrapolated values
+gx, gy, gz = dgrid3d(x, y, z, "40,40 gauss 0.1,0.1")
+@gsp :- 2 tit="Interpolation on a grid\\n(extrapolated values are shown)"  gx gy gz "w l notit lc pal"
+
+# Hide exrapolated values
+gx, gy, gz = dgrid3d(x, y, z, "40,40 gauss 0.1,0.1", extra=false)
+@gsp :- 3 tit="Interpolation on a grid\\n(extrapolated values are hidden)" gx gy gz "w l notit lc pal"
+save(term="pngcairo size 1000,400 fontscale 1.0", output="assets/advanced015b.png")  # hide
+```
+![](assets/advanced015b.png)
+
+
+
 ## Animations
 
 The [Multiplot](@ref) capabilities can also be used to stack plots one above the other in order to create an animation, as in the following example:
diff --git a/docs/src/api.md b/docs/src/api.md
index 14c5b2a..636327e 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -13,6 +13,7 @@ The list of **Gnuplot.jl** exported symbols is as follows:
 boxxy
 contourlines
 dataset_names
+dgrid3d
 gpexec
 gpmargins
 gpranges
diff --git a/docs/src/index.md b/docs/src/index.md
index 1a505ac..a8f2a9d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -21,6 +21,8 @@ The **Gnuplot.jl** package allows easy and fast use of [gnuplot](http://gnuplot.
 
 - enhanced support for contour plots;
 
+- 2D interpolation of scattered data on a regular grid;
+
 - export to a huge number of formats such as `pdf`, `png`, `gif`, ``\LaTeX``, `svg`, etc. (actually all those supported by gnuplot);
 
 - compatibility with Jupyter and Juno;
diff --git a/src/Gnuplot.jl b/src/Gnuplot.jl
index 923eef2..cf3c1fd 100644
--- a/src/Gnuplot.jl
+++ b/src/Gnuplot.jl
@@ -10,7 +10,7 @@ import Base.show
 export session_names, dataset_names, palette_names, linetypes, palette,
     terminal, terminals, test_terminal,
     stats, @gp, @gsp, save, gpexec,
-    boxxy, contourlines, hist, recipe, gpvars, gpmargins, gpranges
+    boxxy, contourlines, dgrid3d, hist, recipe, gpvars, gpmargins, gpranges
 
 
 # ╭───────────────────────────────────────────────────────────────────╮
@@ -2094,6 +2094,105 @@ function contourlines(x::AbstractVector{Float64}, y::AbstractVector{Float64}, z:
 end
 
 
+# --------------------------------------------------------------------
+"""
+    dgrid3d(x, y, z, opts=""; extra=true)
+
+Interpolate non-uniformly spaced 2D data onto a regular grid.
+
+!!! note
+    This feature is not available in *dry* mode and will raise an error if used.
+
+# Arguments:
+- `x`, `y`, `z` (as `AbstractVector{Float64}`): coordinates and values of the function to interpolate;
+- `opts`: interpolation settings (see gnuplot documentation for `dgrid3d`);
+- `extra`: if `true` (default) compute inerpolated values in all regions, even those which are poorly constrained by input data (namely, extrapolated values).  If `false` set these values to `NaN`.
+
+# Return values:
+A tuple with `x` and `y` coordinates on the regular grid (as `Vector{Float64}`), and `z` containing interpolated values (as `Matrix{Float64}`).
+
+# Example
+```julia
+x = (rand(200) .- 0.5) .* 3;
+y = (rand(200) .- 0.5) .* 3;
+z = exp.(-(x.^2 .+ y.^2));
+
+# Interpolate on a 20x30 regular grid with splines
+gx, gy, gz = dgrid3d(x, y, z, "20,30 splines")
+
+@gsp "set size ratio -1" "set xyplane at 0" xlab="X" ylab="Y" :-
+@gsp :-  x  y  z "w p t 'Scattered data' lc pal"
+@gsp :- gx gy gz "w l t 'Interpolation on a grid' lc pal"
+```
+!!! warn
+    The `splines` algorithm may be very slow on large datasets.  An alternative option is to use a smoothing kernel, such as `gauss`:
+
+```julia
+x = randn(2000) .* 0.5;
+y = randn(2000) .* 0.5;
+rsq = x.^2 + y.^2;
+z = exp.(-rsq) .* sin.(y) .* cos.(2 * rsq);
+
+@gsp "set size ratio -1" palette(:balance, smooth=true) "set view map" "set pm3d" :-
+@gsp :- "set multiplot layout 1,3" xr=[-2,2] yr=[-2,2] :-
+@gsp :- 1 tit="Scattered data"  x  y  z "w p notit lc pal"
+
+# Show extrapolated values
+gx, gy, gz = dgrid3d(x, y, z, "40,40 gauss 0.1,0.1")
+@gsp :- 2 tit="Interpolation on a grid\\\\n(extrapolated values are shown)"  gx gy gz "w l notit lc pal"
+
+# Hide exrapolated values
+gx, gy, gz = dgrid3d(x, y, z, "40,40 gauss 0.1,0.1", extra=false)
+@gsp :- 3 tit="Interpolation on a grid\\\\n(extrapolated values are hidden)" gx gy gz "w l notit lc pal"
+```
+"""
+function dgrid3d(x::AbstractVector{Float64},
+                 y::AbstractVector{Float64},
+                 z::AbstractVector{Float64},
+                 opts::String="";
+                 extra=true)
+    c = Gnuplot.gp_write_table("set dgrid3d $opts", x, y, z, is3d=true)
+    gx = Vector{Float64}()
+    gy = Vector{Float64}()
+    gz = Vector{Float64}()
+    ix = 0
+    iy = 0
+    for l in c
+        l = string(strip(l))
+        if l == "# x y z type"
+            ix += 1
+            iy = 1
+        else
+            (l == "" )     &&  continue
+            (l[1] == '#')  &&  continue
+            n = Meta.parse.(split(l)[1:3])
+            (iy == 1)  &&  push!(gx, n[1])
+            (ix == 1)  &&  push!(gy, n[2])
+            if n[3] == :NaN
+                push!(gz, NaN)
+            else
+                push!(gz, n[3])
+            end
+            iy += 1
+        end
+    end
+    gz = collect(reshape(gz, length(gy), length(gx))')
+    if !extra
+        dx = abs(gx[2]-gx[1]) / 2
+        dy = abs(gy[2]-gy[1]) / 2
+        for ix in 1:length(gx)
+            for iy in 1:length(gy)
+                n = length(findall(((gx[ix] - dx) .< x .< (gx[ix] + dx))  .&
+                                   ((gy[iy] - dy) .< y .< (gy[iy] + dy))))
+                (n == 0)  &&  (gz[ix, iy] = NaN)
+            end
+        end
+    end
+    return (gx, gy, gz)
+end
+
+
+
 # ╭───────────────────────────────────────────────────────────────────╮
 # │                        GNUPLOT REPL                               │
 # ╰───────────────────────────────────────────────────────────────────╯