diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
index 1d3aec7..f018c9d 100644
--- a/.github/workflows/CompatHelper.yml
+++ b/.github/workflows/CompatHelper.yml
@@ -9,7 +9,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: [1.8.0]
+        julia-version: ["lts"]
         julia-arch: [x86]
         os: [ubuntu-latest]
     steps:
diff --git a/.github/workflows/RunTests.yml b/.github/workflows/RunTests.yml
index 8f04439..5ab90f3 100644
--- a/.github/workflows/RunTests.yml
+++ b/.github/workflows/RunTests.yml
@@ -7,13 +7,13 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ["1.8"]
+        julia-version: ["1", "1.9", "lts"]
         julia-arch: [x64]
         os: [ubuntu-latest]
 
     steps:
       - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
           arch: ${{ matrix.julia-arch }}
diff --git a/Project.toml b/Project.toml
index f4d1292..2b0eb98 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,17 +1,15 @@
 name = "SpatialBoundaries"
 uuid = "8d2ba62a-3d23-4a2b-b692-6b63e9267be2"
 authors = ["Tanya Strydom <tanya.strydom@umontreal.ca>", "Timothée Poisot <timothee.poisot@umontreal.ca>"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-SpeciesDistributionToolkit = "72b53823-5c0b-4575-ad0e-8e97227ad13b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 VoronoiDelaunay = "72f80fcb-8c52-57d9-aff0-40c1a3526986"
 
 [compat]
-SpeciesDistributionToolkit = "0.0.1, 0.0.2, 0.0.3"
-StatsBase = "0.33"
+StatsBase = "0.34"
 VoronoiDelaunay = "0.4"
-julia = "1.8"
+julia = "1.9"
diff --git a/docs/Project.toml b/docs/Project.toml
index af3b078..129ade2 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,11 +1,6 @@
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 NeutralLandscapes = "71847384-8354-4223-ac08-659a5128069f"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-SpeciesDistributionToolkit = "72b53823-5c0b-4575-ad0e-8e97227ad13b"
-StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
-
-[compat]
-SpeciesDistributionToolkit = "0.0"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
diff --git a/docs/make.jl b/docs/make.jl
index 1265fdb..56111c4 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -15,6 +15,7 @@ makedocs(;
     sitename = "SpatialBoundaries",
     format = Documenter.HTML(),
     modules = [SpatialBoundaries],
+    warnonly = true,
     pages = [
         "Home" => [
             "Introduction" => "vignettes/introduction.md",
@@ -23,18 +24,12 @@ makedocs(;
         "Vignettes" => [
             "Finding boundaries" => "vignettes/boundaries.md",
             "Triangulation wombling" => "vignettes/triangulation.md",
-            "SDM Layers" => "vignettes/simplesdmlayers.md",
         ],
     ],
 )
 
-@info "Prepare to remove"
-run(`find . -type f -name ".tif" -delete`)
-run(`find . -type f -name ".tif"`)
-
 @info "Prepare to deploy"
 deploydocs(;
-    deps = Deps.pip("pygments", "python-markdown-math"),
     repo = "github.com/PoisotLab/SpatialBoundaries.jl.git",
     push_preview = true,
     devbranch = "main",
diff --git a/docs/src/vignettes/boundaries.jl b/docs/src/vignettes/boundaries.jl
index 3822b84..7b9bb6c 100644
--- a/docs/src/vignettes/boundaries.jl
+++ b/docs/src/vignettes/boundaries.jl
@@ -6,8 +6,7 @@
 # three-patch landscape.
 
 using SpatialBoundaries
-using StatsPlots
-default(; dpi=500, size=(600, 600), aspectratio=1, c=:batlow, frame=:box)
+using CairoMakie
 
 # Let's create a landscape with two values:
 
@@ -26,7 +25,7 @@ W = wombling(A);
 
 # Let's look at the rate of change:
 
-heatmap(W.m; c=:nuuk)
+heatmap(W.m; colormap = :Greys)
 
 # Picking the boundaries is done by passing the wombling output to the
 # `boundaries` function, with a specific threshold giving the proportion of
@@ -36,9 +35,11 @@ heatmap(W.m; c=:nuuk)
 thresholds = LinRange(0.0, 0.2, 200)
 patches = [length(boundaries(W, t)) for t in thresholds]
 
-plot(thresholds, log1p.(patches), aspectratio=:none)
-xaxis!("Threshold", (0., 0.2))
-yaxis!("log(boundary patches + 1)", (0., 9.))
+f, ax, plt = lines(thresholds, log1p.(patches), color=:black, linewidth=2)
+tightlimits!(ax)
+ax.xlabel = "Threshold"
+ax.ylabel = "log(boundary patches + 1)"
+f
 
 # Let's eyeball this as 0.01, and see how the patches are distributed.
 
@@ -49,11 +50,11 @@ yaxis!("log(boundary patches + 1)", (0., 9.))
 
 b = similar(W.m)
 
-for t in reverse(LinRange(0.0, 1.0, 200))
+for t in reverse(LinRange(0.0, 1.0, 300))
     b[boundaries(W, t)] .= t
 end
 
-heatmap(b, c=:tofino, clim=(0,1))
+heatmap(b, colormap=:tofino, colorrange=(0,1))
 
 # This also suggests that we will get well delineated patches for low values of
 # the threshold.
@@ -63,7 +64,8 @@ B = boundaries(W, 0.01);
 # In the following figure, cells identified as candidate boundaries are marked
 # in white:
 
-heatmap(A)
-scatter!([(reverse(x.I)) for x in B], leg=false, msw=0, c=:white)
+f, ax, hm = heatmap(A)
+scatter!(ax, [b.I for b in B], color=:white)
+f
 
 # We can see that the boundaries of the patches have been well identified!
diff --git a/docs/src/vignettes/introduction.jl b/docs/src/vignettes/introduction.jl
index e97a331..4730a3f 100644
--- a/docs/src/vignettes/introduction.jl
+++ b/docs/src/vignettes/introduction.jl
@@ -42,11 +42,7 @@
 
 using SpatialBoundaries
 using NeutralLandscapes
-using StatsPlots
-
-# We will set a few options for the default plots:
-
-default(; dpi=500, size=(600, 600), aspectratio=1, c=:davos, frame=:box)
+using CairoMakie
 
 # The landscape generation is done using the `NeutralLandscapes` package, and we
 # will pick a 500x500 grid:
@@ -76,7 +72,7 @@ W = wombling(landscape);
 
 # Let's have a look at the rate of change:
 
-heatmap(W.m, c=:tokyo, clim=(0, maximum(W.m)))
+heatmap(W.m, colormap=:navia, colorrange=(0, maximum(W.m)))
 
 # The rate of change informs us on the potential for there to be a boundary
 # (zone of change) within a window. Cells with a high rate of change are
@@ -87,7 +83,7 @@ heatmap(W.m, c=:tokyo, clim=(0, maximum(W.m)))
 # for instance, an angle of 180° means that the value is smaller in the South,
 # and larger in the North:
 
-heatmap(W.θ, c=:romaO, clim=(0., 360.))
+heatmap(W.θ, colormap=:romaO, colorrange=(0., 360.))
 
 # The direction of change is *not* the direction the boundary would be if you
 # were to draw it on the landscape but rather the direction the rate of change
diff --git a/docs/src/vignettes/simplesdmlayers.jl b/docs/src/vignettes/simplesdmlayers.jl
deleted file mode 100644
index 62f307e..0000000
--- a/docs/src/vignettes/simplesdmlayers.jl
+++ /dev/null
@@ -1,143 +0,0 @@
-# # Integration with SimpleSDMLayers
-
-# Below is an example using the various functions within `SpatialBoundaries` to
-# estimate boundaries for (*i.e.* patches of) wooded areas on the Southwestern
-# islands of the Hawaiian Islands using landcover data from the EarthEnv project
-# [@Tuanmu2014Glo1km] as well as integrating some functionality from
-# `SimpleSDMLayers` [@Dansereau2021SimJl] for easier work with the spatial
-# nature of the input data. The `SpatialBoundaries` package works really well
-# with `SimpleSDMLayers`, so that you can (i) apply wombling and boundaries
-# finding to a `SimpleSDMLayer` object, and (ii) convert the output of a
-# `Womble` object to a *pair* of `SimpleSDMLayer` corresponding to the rate and
-# direction of change.
-
-# Because there are four different layers in the EarthEnv database that
-# represent different types of woody cover we will use the overall mean wombling
-# value. As the data are arranged in a matrix *i.e.* a lattice this example will
-# focus on lattice wombling, however for triangulation wombling the
-# implementation of functions and workflow would look similar with the exception
-# that the input data would be structured differently (as three vectors of $x$,
-# $y$, $z$) and the output data would be typed as `TriangulationWomble` objects.
-
-using SpatialBoundaries
-using SpeciesDistributionToolkit
-using CairoMakie
-import Plots # Required for radial histograms
-
-# First we can start by defining the extent of the Southwestern islands of
-# Hawaii, which can be used to restrict the extraction of the various landcover
-# layers from the EarthEnv database. We do the actual layers querying using
-# `SimpleSDMLayers`:
-
-hawaii = (left = -160.2, right = -154.5, bottom = 18.6, top = 22.5)
-dataprovider = RasterData(EarthEnv, LandCover)
-landcover_classes = SimpleSDMDatasets.layers(dataprovider)
-landcover = [SimpleSDMPredictor(dataprovider; layer=class, full=true, hawaii...) for class in landcover_classes]
-
-# We can remove all the areas that contain 100% water from the landcover data as
-# our question of interest is restricted to the terrestrial realm. We do this by
-# using the "Open Water" layer to mask over each of the landcover layers
-# individually:
-
-ow_index = findfirst(isequal("Open Water"), landcover_classes)
-not_water = landcover[ow_index] .!== 0x64
-lc = [mask(not_water, layer) for layer in landcover]
-
-# As layers one through four of the EarthEnv data are concerned with data on
-# woody cover (*i.e.* "Evergreen/Deciduous Needleleaf Trees", "Evergreen
-# Broadleaf Trees", "Deciduous Broadleaf Trees", and "Mixed/Other Trees") we
-# will work with only these layers. For a quick idea we of what the raw
-# landcover data looks like we can sum these four layers and plot the total
-# woody cover:
-
-classes_with_trees = findall(contains.(landcover_classes, "Trees"))
-tree_lc = convert(Float32, reduce(+, lc[classes_with_trees]))
-heatmap(tree_lc; colormap=:linear_kbgyw_5_98_c62_n256)
-
-# Although we have previously summed the four landcover layers for the actual
-# wombling part we will apply the wombling function to each layer before we
-# calculate the overall mean wombling value. We will do this using `broadcast`,
-# this will apply `wombling` in an element-wise fashion to the four different
-# woody cover layers. This will give as a vector containing four `LatticeWomble`
-# objects (since the input data was in the form of a matrix).
-
-wombled_layers = wombling.(lc[classes_with_trees]);
-
-# As we are interested in overall woody cover for Southwestern islands we can
-# take the `wombled_layers` vector and use them with the `mean` function to get
-# the overall mean wombling value of the rate and direction of change for woody
-# cover. This will 'flatten' the four wombled layers into a single
-# `LatticeWomble` object.
-
-wombled_mean = mean(wombled_layers);
-
-# From the `wombled_mean` object we can 'extract' the layers for both the mean
-# rate and direction of change. For ease of plotting we will also convert these
-# layers to `SimpleSDMPredictor` type objects. It is also possible to call these
-# matrices directly from the `wombled_mean` object using `wombled_mean.m` or
-# `wombled_mean.θ` for the rate and direction respectively.
-
-rate, direction = SimpleSDMPredictor(wombled_mean)
-
-# Lastly we can identify candidate boundaries using the `boundaries`. Here we
-# will use a thresholding value (t) of 0.1 and save these candidate boundary
-# cells as `b`. Note that we are now working with a `SimpleSDMResponse` object
-# and this is simply for ease of plotting.
-
-b = similar(rate)
-b.grid[boundaries(wombled_mean, 0.1; ignorezero = true)] .= 1.0
-
-# We will overlay the identified boundaries (in green) over the rate of change
-# (in levels of grey):
-
-heatmap(rate, colormap=[:grey95, :grey5])
-heatmap!(b, colormap=[:transparent, :green])
-current_figure()
-
-# For this example we will plot the direction of change as radial plots to get
-# an idea of the prominent direction of change. Here we will plot *all* the
-# direction values from `direction` for which the rate of change is greater than
-# zero (so as to avoid denoting directions for a slope that does not exist) as
-# well as the `direction` values from only candidate cells using the same
-# masking principle as what we did for the rate of change. It is of course also
-# possible to forgo the radial plots and plot the direction of change in the
-# same manner as the rate of change should one wish.
-
-# Before we plot let us create our two 'masked layers'. For all direction values
-# for which there is a corresponding rate of change greater than zero we can use
-# `rate` as a masking layer but first replace all zero values with 'nothing'.
-# For the candidate boundary cells we can simply mask `direction` with `b` as we
-# did for the rate of change.
-
-direction_all = mask(replace(rate, 0 => nothing), direction)
-
-direction_candidate = mask(b, direction)
-
-# Because stephist() requires a vector of radians for plotting we must first
-# collect the cells and convert them from degrees to radians. Then we can start
-# by plotting the direction of change of *all* cells.
-
-Plots.stephist(
-         deg2rad.(values(direction_all));
-         proj=:polar,
-         lab="",
-         c=:teal,
-         nbins = 36,
-         yshowaxis=false,
-         normalize = false)
-
-# Followed by plotting the direction of change only for cells that are
-# considered as candidate boundary cells.
-
-Plots.stephist(
-        deg2rad.(values(direction_candidate));
-        proj=:polar,
-        lab="",
-        c=:red,
-        nbins = 36,
-        yshowaxis=false,
-        normalize = false)
-
-# End
-
-rm(joinpath(SimpleSDMLayers._layers_assets_path, "EarthEnv"); force=true, recursive=true) #hide
\ No newline at end of file
diff --git a/docs/src/vignettes/triangulation.jl b/docs/src/vignettes/triangulation.jl
index 836e4d6..b31a565 100644
--- a/docs/src/vignettes/triangulation.jl
+++ b/docs/src/vignettes/triangulation.jl
@@ -1,20 +1,17 @@
 # # Triangulation wombling
 
 using SpatialBoundaries
-using StatsPlots
-
-# Plot defaults
-
-default(; dpi=500, size=(600, 600), aspectratio=1, c=:davos, frame=:box)
+using CairoMakie
+using StatsBase
 
 # Get some points at random
 
 n = 500
-x = rand(n)
+x = 2rand(n)
 y = rand(n)
-z = [(x[i]<=0.5)&(y[i]<=0.5) ? rand() : rand().+1.2 for i in eachindex(x)]
+z = [(x[i] <= 0.5) & (y[i] <= 0.5) ? rand() : rand() .+ 1.2 for i in eachindex(x)]
 
-scatter(x, y, marker_z = z, lab="")
+scatter(x, y; color = z)
 
 # Womble
 
@@ -22,20 +19,18 @@ W = wombling(x, y, z)
 
 # Get the rate of change
 
-scatter(x, y, c=:lightgrey, msw=0.0, lab="", m=:square, ms=3)
-scatter!(W.x, W.y, marker_z = log1p.(W.m), lab="")
+scatter(x, y; color = :lightgrey)
+scatter!(W.x, W.y; color = log1p.(W.m))
+current_figure()
 
 # Angle histogram
 
-stephist(
-    deg2rad.(sort(vec(W.θ)));
-    proj=:polar,
-    lab="",
-    c=:teal,
-    fill=(0, 0.2, :teal),
-    nbins=100,
-)
+f = Figure()
+ax = PolarAxis(f[1, 1])
+h = fit(Histogram, deg2rad.(values(W.θ)); nbins = 100)
+stairs!(ax, h.edges[1], h.weights[[axes(h.weights, 1)..., 1]]; color = :teal, linewidth = 2)
+current_figure()
 
 # Show the rotation with a color
 
-scatter(W.x, W.y, marker_z = W.θ, c=:vik, clim=(0, 360))
\ No newline at end of file
+scatter(W.x, W.y; color = W.θ, colormap = :vikO, colorrange = (0, 360))
\ No newline at end of file
diff --git a/src/SpatialBoundaries.jl b/src/SpatialBoundaries.jl
index 01c654f..caa1b11 100644
--- a/src/SpatialBoundaries.jl
+++ b/src/SpatialBoundaries.jl
@@ -5,7 +5,6 @@ using VoronoiDelaunay
 using LinearAlgebra
 using Statistics
 using StatsBase
-using SpeciesDistributionToolkit
 
 include(joinpath("types", "Womble.jl"))
 export Womble, TriangulationWomble, LatticeWomble
@@ -13,7 +12,6 @@ export Womble, TriangulationWomble, LatticeWomble
 include(joinpath("lib", "rateofchange.jl"))
 
 include(joinpath("lib", "wombling.jl"))
-include(joinpath("extensions", "SpeciesDistributionToolkit.jl"))
 export wombling
 
 include(joinpath("lib", "boundaries.jl"))
diff --git a/src/extensions/SpeciesDistributionToolkit.jl b/src/extensions/SpeciesDistributionToolkit.jl
deleted file mode 100644
index 12c867f..0000000
--- a/src/extensions/SpeciesDistributionToolkit.jl
+++ /dev/null
@@ -1,37 +0,0 @@
-"""
-    wombling(layer::T; convert_to::Type=Float64) where {T <: SimpleSDMLayer}
-
-Performs a lattice wombling on a `SimpleSDMLayer`.
-"""
-function wombling(layer::T; convert_to::Type = Float64) where {T <: SimpleSDMLayer}
-    try
-        global nan = convert(convert_to, NaN)
-    catch
-        throw(ArgumentError("The type given as `convert_to` must have a `NaN` value."))
-    end
-
-    # Get the values for x and y
-    y = collect(SimpleSDMLayers.longitudes(layer))
-    x = collect(SimpleSDMLayers.latitudes(layer))
-
-    # Get the grid
-    z = convert(Matrix{Union{Nothing, convert_to}}, layer.grid)
-    replace!(z, nothing => nan)
-    return wombling(x, y, convert(Matrix{convert_to}, z))
-end
-
-# We want to make sure that the layers are returned without NaN values, and
-# adding this method makes the code easier to write
-Base.isnan(::Nothing) = false
-
-function SimpleSDMLayers.SimpleSDMPredictor(W::T) where {T <: LatticeWomble}
-    rate = SimpleSDMLayers.SimpleSDMPredictor(W.m, extrema(W.y)..., extrema(W.x)...)
-    direction = SimpleSDMLayers.SimpleSDMPredictor(W.θ, extrema(W.y)..., extrema(W.x)...)
-    rate.grid[findall(isnan, rate.grid)] .= nothing
-    direction.grid[findall(isnan, direction.grid)] .= nothing
-    return (rate, direction)
-end
-
-function SimpleSDMLayers.SimpleSDMResponse(W::T) where {T <: LatticeWomble}
-    return convert.(SimpleSDMLayers.SimpleSDMResponse, SimpleSDMLayers.SimpleSDMPredictor(W))
-end
diff --git a/test/Project.toml b/test/Project.toml
index b0487c8..e68d3b1 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,3 @@
 [deps]
 NeutralLandscapes = "71847384-8354-4223-ac08-659a5128069f"
-SpeciesDistributionToolkit = "72b53823-5c0b-4575-ad0e-8e97227ad13b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/runtests.jl b/test/runtests.jl
index 820f76f..b0ae9a2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,7 +7,6 @@ tests = [
     "Triangulation wombling" => "03_delaunay.jl",
     "Direction of change" => "04_direction.jl",
     "Overall mean wombling" => "05_mean.jl",
-    "REQUIRE: SimpleSDMLayers" => "R1_simplesdmlayers.jl",
 ]
 
 for test in tests
diff --git a/test/units/R1_simplesdmlayers.jl b/test/units/R1_simplesdmlayers.jl
deleted file mode 100644
index 4debfe9..0000000
--- a/test/units/R1_simplesdmlayers.jl
+++ /dev/null
@@ -1,29 +0,0 @@
-module SBTestSimpleSDMLayers
-
-using SpatialBoundaries
-using SpeciesDistributionToolkit
-using Test
-
-precipitation = SimpleSDMPredictor(
-    RasterData(WorldClim2, BioClim),
-    layer=12;
-    left = -80.0,
-    right = -56.0,
-    bottom = 44.0,
-    top = 62.0,
-)
-
-W = wombling(precipitation)
-
-@test isa(W, LatticeWomble)
-
-Lt, Ld = SimpleSDMPredictor(W)
-@test isa(Lt, SimpleSDMPredictor)
-
-Rt, Rd = SimpleSDMResponse(W)
-@test isa(Rt, SimpleSDMResponse)
-
-@test !any(map(isnan, Rt.grid))
-@test !any(map(isnan, Rd.grid))
-
-end