diff --git a/Project.toml b/Project.toml
index 3aa9534..b4ad60b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "WorldOceanAtlasTools"
 uuid = "04f20302-f1b9-11e8-29d9-7d841cb0a64a"
 authors = ["Benoit Pasquier <pasquieb@uci.edu>"]
-version = "0.3.4"
+version = "0.3.5"
 
 [deps]
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
diff --git a/src/functions.jl b/src/functions.jl
index b542daa..e00c448 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -86,6 +86,19 @@ end
 get_unit(ds, tracer, field) = convert_to_Unitful(ds[WOA_varname(tracer, field)].attrib["units"])
 
 function mean_and_variance_gridded_3d_field(grid::OceanGrid, field3D, lat, lon, depth)
+    χ_3D, σ²_3D, n_3D = raw_mean_and_variance_gridded_3d_field(grid, field3D, lat, lon, depth)
+    # Enforce μ = 0 and σ = ∞ where no observations
+    # (Instead of NaNs)
+    println("  Setting μ = 0 and σ² = ∞ where no obs")
+    χ_3D[findall(n_3D .== 0)] .= 0.0
+    σ²_3D[findall(n_3D .== 0)] .= Inf
+    println("  Setting a realistic minimum for σ²")
+    meanχ = mean(χ_3D, weights(n_3D))
+    σ²_3D .= max.(σ²_3D, 1e-4meanχ^2)
+    return χ_3D, σ²_3D
+end
+
+function raw_mean_and_variance_gridded_3d_field(grid::OceanGrid, field3D, lat, lon, depth)
     fieldvec, latvec, lonvec, depthvec, CI = filter_gridded_3D_field(field3D, lat, lon, depth)
     println("  Averaging data over each grid box")
     χ_3D = zeros(size(grid))
@@ -103,15 +116,7 @@ function mean_and_variance_gridded_3d_field(grid::OceanGrid, field3D, lat, lon,
     end
     χ_3D .= χ_3D ./ n_3D               # μ = Σᵢ μᵢ / n
     σ²_3D .= σ²_3D ./ n_3D .- χ_3D.^2  # σ² = Σᵢ μᵢ² / n - μ²
-    # Enforce μ = 0 and σ = ∞ where no observations
-    # (Instead of NaNs)
-    println("  Setting μ = 0 and σ² = ∞ where no obs")
-    χ_3D[findall(n_3D .== 0)] .= 0.0
-    σ²_3D[findall(n_3D .== 0)] .= Inf
-    println("  Setting a realistic minimum for σ²")
-    meanχ = mean(χ_3D, weights(n_3D))
-    σ²_3D .= max.(σ²_3D, 1e-4meanχ^2)
-    return χ_3D, σ²_3D
+    return χ_3D, σ²_3D, n_3D 
 end
 
 function convert_to_SI_unit!(χ_3D, σ²_3D, ds, tracer, field)
@@ -129,6 +134,14 @@ function fit_to_grid(grid::OceanGrid, product_year, tracer, period, resolution,
     return χ_3D, σ²_3D
 end
 
+function raw_to_grid(grid::OceanGrid, product_year, tracer, period, resolution, field)
+    ds = WOA_Dataset(product_year, tracer, period, resolution)
+    field3D, lat, lon, depth = get_gridded_3D_field(ds, tracer, field)
+    χ_3D, σ²_3D, n_3D = raw_mean_and_variance_gridded_3d_field(grid, field3D, lat, lon, depth)
+    convert_to_SI_unit!(χ_3D, σ²_3D, ds, tracer, field)
+    return χ_3D, n_3D
+end
+
 #==================================
 Helper functions
 ==================================#