diff --git a/src/backends/gr.jl b/src/backends/gr.jl
index 7aee8768..e7ab5209 100644
--- a/src/backends/gr.jl
+++ b/src/backends/gr.jl
@@ -920,11 +920,8 @@ function gr_display(sp::Subplot{GRBackend}, w, h, viewport_canvas)
             outside_ticks = true
             for ax in (sp[:xaxis], sp[:yaxis])
                 v = series[ax[:letter]]
-                if length(v) > 1 && diff(collect(extrema(diff(v))))[1] > 1e-6*std(v)
-                    @warn("GR: heatmap only supported with equally spaced data.")
-                end
             end
-            x, y = heatmap_edges(series[:x], sp[:xaxis][:scale]), heatmap_edges(series[:y], sp[:yaxis][:scale])
+            x, y = heatmap_edges(series[:x], sp[:xaxis][:scale], series[:y], sp[:yaxis][:scale], size(series[:z]))
             xy_lims = x[1], x[end], y[1], y[end]
             expand_extrema!(sp[:xaxis], x)
             expand_extrema!(sp[:yaxis], y)
@@ -1317,23 +1314,31 @@ function gr_display(sp::Subplot{GRBackend}, w, h, viewport_canvas)
         elseif st == :heatmap
             zmin, zmax = clims
             if !ispolar(sp)
-                xmin, xmax, ymin, ymax = xy_lims
-                m, n = length(x), length(y)
                 GR.setspace(zmin, zmax, 0, 90)
-                grad = isa(series[:fillcolor], ColorGradient) ? series[:fillcolor] : cgrad()
-                colors = [plot_color(grad[clamp((zi-zmin) / (zmax-zmin), 0, 1)], series[:fillalpha]) for zi=z]
-                rgba = map(c -> UInt32( round(UInt, alpha(c) * 255) << 24 +
-                                        round(UInt,  blue(c) * 255) << 16 +
-                                        round(UInt, green(c) * 255) << 8  +
-                                        round(UInt,   red(c) * 255) ), colors)
-                w, h = length(x), length(y)
-                GR.drawimage(xmin, xmax, ymax, ymin, w, h, rgba)
+                x, y = heatmap_edges(series[:x], sp[:xaxis][:scale], series[:y], sp[:yaxis][:scale], size(series[:z]))
+                w, h = length(x) - 1, length(y) - 1
+                z_normalized = map(x -> GR.jlgr.normalize_color(x, zmin, zmax), z)
+                z_normalized = map(x -> isnan(x) ? 256/255 : x, z_normalized) # results in color index = 1256 -> transparent
+                colors = Int32[round(Int32, 1000 + _i * 255) for _i in z_normalized]
+                GR.nonuniformcellarray(x, y, w, h, colors)
             else
-                h, w = length(x), length(y)
-                z = reshape(z, h, w)
-                colors = Int32[round(Int32, 1000 + _i * 255) for _i in z']
-                GR.setwindow(-1, 1, -1, 1)
-                GR.polarcellarray(0, 0, 0, 360, 0, 1, w, h, colors)
+                phimin, phimax = 0.0, 360.0 # nonuniform polar array is not yet supported in GR.jl
+                nx, ny = length(series[:x]), length(series[:y])
+                z_normalized = map(x -> GR.jlgr.normalize_color(x, zmin, zmax), z)
+                colors = Int32[round(Int32, 1000 + _i * 255) for _i in z_normalized]
+                xmin, xmax, ymin, ymax = xy_lims 
+                rmax = data_lims[4]
+                GR.setwindow(-rmax, rmax, -rmax, rmax)
+                if ymin > 0 
+                    @warn "'ymin[1] > 0' (rmin) is not yet supported."
+                end
+                if series[:y][end] != ny
+                    @warn "Right now only the maximum value of y (r) is taken into account."
+                end
+                # GR.polarcellarray(0, 0, phimin, phimax, ymin, ymax, nx, ny, colors)
+                GR.polarcellarray(0, 0, phimin, phimax, 0, ymax, nx, ny, colors)
+                # Right now only the maximum value of y (r) is taken into account. 
+                # This is certainly not perfect but nonuniform polar array is not yet supported in GR.jl
             end
 
         elseif st in (:path3d, :scatter3d)
diff --git a/src/utils.jl b/src/utils.jl
index 0c2b4b80..324aa2cb 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -342,8 +342,11 @@ const _scale_base = Dict{Symbol, Real}(
     :ln => ℯ,
 )
 
-function _heatmap_edges(v::AVec)
+function _heatmap_edges(v::AVec, isedges::Bool = false)
   length(v) == 1 && return v[1] .+ [-0.5, 0.5]
+  if isedges return v end 
+  # `isedges = true` means that v is a vector which already describes edges 
+  # and does not need to be extended.
   vmin, vmax = ignorenan_extrema(v)
   extra_min = (v[2] - v[1]) / 2
   extra_max = (v[end] - v[end - 1]) / 2
@@ -351,9 +354,25 @@ function _heatmap_edges(v::AVec)
 end
 
 "create an (n+1) list of the outsides of heatmap rectangles"
-function heatmap_edges(v::AVec, scale::Symbol = :identity)
+function heatmap_edges(v::AVec, scale::Symbol = :identity, isedges::Bool = false)
   f, invf = scalefunc(scale), invscalefunc(scale)
-  map(invf, _heatmap_edges(map(f,v)))
+  map(invf, _heatmap_edges(map(f,v), isedges))
+end
+
+function heatmap_edges(x::AVec, xscale::Symbol, y::AVec, yscale::Symbol, z_size::Tuple{Int, Int})
+    nx, ny = length(x), length(y)
+    # ismidpoints = z_size == (ny, nx) # This fails some tests, but would actually be 
+    # the correct check, since (4, 3) != (3, 4) and a missleading plot is produced.
+    ismidpoints = prod(z_size) == (ny * nx) 
+    isedges = z_size == (ny - 1, nx - 1)
+    if !ismidpoints && !isedges
+        error("""Length of x & y does not match the size of z. 
+                Must be either `size(z) == (length(y), length(x))` (x & y define midpoints)
+                or `size(z) == (length(y)+1, length(x)+1))` (x & y define edges).""")
+    end
+    x, y = heatmap_edges(x, xscale, isedges), 
+           heatmap_edges(y, yscale, isedges)
+    return x, y
 end
 
 function convert_to_polar(theta, r, r_extrema = ignorenan_extrema(r))