diff --git a/docs/Project.toml b/docs/Project.toml
index 1dda9289..0698cb8f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -18,6 +18,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
+NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
diff --git a/docs/src/ebm_melt/ebm_melt.md b/docs/src/ebm_melt/ebm_melt.md
index ea2a56df..847042bb 100644
--- a/docs/src/ebm_melt/ebm_melt.md
+++ b/docs/src/ebm_melt/ebm_melt.md
@@ -20,6 +20,7 @@ using CoordRefSystems
 using CairoMakie
 using LinearAlgebra
 using MLStyle
+using NearestNeighbors
 using NetCDF
 using OrdinaryDiffEq
 Point3D = Point3{Float64}
@@ -61,9 +62,11 @@ p_sph = map(point(s)) do p
 end
 
 # Note: You can instead use an algebraic parameterization, rather than nearest-neighbor interpolation.
-# Note: You can alternatively set a value to 0.0 if the distance to the nearest-neighbor is greater than some threshold.
+lat, lon = lat[:], lon[:]
+ll = hcat(lat, lon)'
+kdt = KDTree(ll)
 sit_sph_idxs = map(p_sph) do p
-  argmin(map(i -> sqrt((lat[i] - p[1])^2 + (lon[i] - p[2])^2), eachindex(lat)))
+  nn(kdt, p)[1]
 end
 
 sit_sph = map(sit_sph_idxs, p_sph) do i, p