diff --git a/ext/GLMakie_Visualisation.jl b/ext/GLMakie_Visualisation.jl
index 5e07e998..6f130a35 100644
--- a/ext/GLMakie_Visualisation.jl
+++ b/ext/GLMakie_Visualisation.jl
@@ -333,7 +333,7 @@ end
 in-place heatmap for a 2D GeoData object (surface), 
 """
 function heatmap!(x::GeoData, args...; field=:Topography, kwargs...)
-    @assert size(x.z.val,3)==1
+    @assert size(x.depth.val,3)==1
     @assert hasfield(typeof(x.fields), field)
 
     heatmap(x.lon.val[:,1], x.lat.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...)
@@ -345,7 +345,7 @@ end
 in-place heatmap for a 2D GeoData object (surface) with an array `a`. 
 """
 function heatmap!(x::GeoData, a::Array{_T,N}, args...; kwargs...) where{_T,N}
-    @assert size(x.z.val,3)==1
+    @assert size(x.depth.val,3)==1
 
     if N==3
         heatmap!(x.lon.val[:,1], x.lat.val[1,:], ustrip.(a[:,:,1]), args...; kwargs...)
@@ -356,14 +356,23 @@ function heatmap!(x::GeoData, a::Array{_T,N}, args...; kwargs...) where{_T,N}
 end
 
 """
-    heatmap!(x::CartData, args...; field=:Topography, kwargs...)
+    heatmap!(x::CartData, args...; field=:Topography, colorbar=false, kwargs...)
 in-place heatmap for a 2D CartData object (surface)
 """
-function heatmap!(x::CartData, args...; field=:Topography, kwargs...)
+function heatmap!(x::CartData, args...; field=:Topography, colorbar=false, kwargs...)
     @assert size(x.z.val,3)==1
     @assert hasfield(typeof(x.fields), field)
 
-    heatmap!(x.x.val[:,1], x.y.val[1,:], ustrip.(x.fields[field][:,:,1]), args...; kwargs...)
+    data = ustrip.(x.fields[field][:,:,1])
+
+    fig,ax,hm = heatmap!(x.x.val[:,1], x.y.val[1,:], data, args...; kwargs...)
+
+    cb = nothing
+    if colorbar
+        cb = Colorbar(fig[1,2], limits=extrema(filter(!isnan,x.fields[field])), label=field)
+    end
+    
+    return fig_out
 end
 
 """
diff --git a/src/utils.jl b/src/utils.jl
index 5aab2793..fc276bf6 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1216,6 +1216,34 @@ function interpolate_datafields_2D(Original::CartData, New::CartData; Rotate=0.0
     return CartData(New.x.val,New.y.val,Znew, fields_new)
 end
 
+"""
+    interpolate_datafields_2D(Original::GeoData, New::GeoData; Rotate=0.0, Translate=(0,0,0), Scale=(1.0,1.0,1.0))
+
+Interpolates a data field `Original` on a 2D GeoData grid `New`.
+Typically used for horizontal surfaces.
+
+Note: `Original` should have orthogonal coordinates. If it has not, e.g., because it was rotated, you'll have to specify the angle `Rotate` that it was rotated by
+
+"""
+function interpolate_datafields_2D(Original::GeoData, New::GeoData; Rotate=0.0, Translate=(0.0,0.0,0.0), Scale=(1.0,1.0,1.0))
+    if (Rotate!=0.0) || any(Translate .!= (0,0,0)) || any(Scale .!= (1.0,1.0,1.0))
+        Original_r  = rotate_translate_scale(Original, Rotate = -1.0*Rotate, Translate = -1.0.*Translate, Scale=Scale);
+        New_r       = rotate_translate_scale(New, Rotate = -1.0*Rotate, Translate = -1.0.*Translate, Scale=Scale);
+    else
+        Original_r  = Original;
+        New_r       = New;
+    end
+
+    Lon_vec      =  Original_r.lon.val[:,1,1];
+    Lat_vec      =  Original_r.lat.val[1,:,1];
+
+    Lon_new = New_r.lon.val
+    Lat_new = New_r.lat.val
+    Znew, fields_new = GeophysicalModelGenerator.InterpolateDataFields2D_vecs(Lon_vec, Lat_vec, Original_r.depth, Original_r.fields, Lon_new, Lat_new)
+
+    return GeoData(New.lon.val,New.lat.val,Znew, fields_new)
+end
+
 """    
     Surf_interp = interpolate_datafields_2D(V::GeoData, x::AbstractRange, y::AbstractRange;  Lat::Number, Lon::Number)