diff --git a/Project.toml b/Project.toml
index d1eb704cc..2a6b43ec6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -97,7 +97,8 @@ StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92"
 
 [targets]
-test = ["Colors", "Distributions", "FileIO", "FilePathsBase", "Gaston", "GeometryBasics", "Gtk", "ImageMagick", "Images", "InspectDR", "LibGit2", "OffsetArrays", "PGFPlotsX", "PlotlyJS", "PlotlyBase", "PyPlot", "PlotlyKaleido", "HDF5", "RDatasets", "StableRNGs", "StaticArrays", "StatsPlots", "Test", "TestImages", "UnicodePlots", "VisualRegressionTests"]
+test = ["Colors", "Distributions", "FileIO", "FilePathsBase", "Gaston", "GeometryBasics", "Gtk", "ImageMagick", "Images", "InspectDR", "LibGit2", "OffsetArrays", "PGFPlotsX", "PlotlyJS", "PlotlyBase", "PyPlot", "PlotlyKaleido", "HDF5", "RDatasets", "StableRNGs", "StaticArrays", "StatsPlots", "Test", "TestImages", "UnicodePlots", "Unitful", "VisualRegressionTests"]
diff --git a/src/Plots.jl b/src/Plots.jl
index 5bc8c2093..c17a62707 100644
--- a/src/Plots.jl
+++ b/src/Plots.jl
@@ -1,7 +1,5 @@
 module Plots
 
-using Pkg
-
 if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
     @eval Base.Experimental.@optlevel 1
 end
@@ -9,9 +7,12 @@ if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@max_m
     @eval Base.Experimental.@max_methods 1
 end
 
+using Pkg
+
 const _plots_project = Pkg.Types.read_project(normpath(@__DIR__, "..", "Project.toml"))
 const _current_plots_version = _plots_project.version
 const _plots_compats = _plots_project.compat
+
 function _check_compat(sim::Module)
     sim_str = string(sim)
     haskey(_plots_compats, sim_str) || return nothing
@@ -28,26 +29,42 @@ function _check_compat(sim::Module)
     end
 end
 
-using Reexport
-
-using Dates, Printf, Statistics, Base64, LinearAlgebra, Random, Unzip
+using Dates, Printf, Statistics, Base64, LinearAlgebra, Random
 using SparseArrays
-
-using FFMPEG
-
-@reexport using RecipesBase
-import RecipesBase: plot, plot!, animate, is_explicit, grid
 using Base.Meta
-@reexport using PlotUtils
+using Requires
+using Reexport
+using Unzip
+@reexport using RecipesBase
 @reexport using PlotThemes
+@reexport using PlotUtils
+
+import RecipesBase: plot, plot!, animate, is_explicit, grid
+import RecipesPipeline
+import RecipesPipeline:
+    inverse_scale_func,
+    datetimeformatter,
+    AbstractSurface,
+    group_as_matrix, # for StatsPlots
+    dateformatter,
+    timeformatter,
+    needs_3d_axes,
+    DefaultsDict,
+    scale_func,
+    is_surface,
+    Formatted,
+    reset_kw!,
+    SliceIt,
+    Surface,
+    pop_kw!,
+    Volume,
+    is3d
 import UnicodeFun
 import StatsBase
 import Downloads
 import Showoff
-import JSON
 import JLFzf
-
-using Requires
+import JSON
 
 #! format: off
 export
@@ -142,56 +159,19 @@ ignorenan_extrema(x::AbstractArray{<:AbstractFloat}) = NaNMath.extrema(x)
 ignorenan_extrema(x) = Base.extrema(x)
 
 # ---------------------------------------------------------
-
-# to cater for block matrices, Base.transpose is recursive.
-# This makes it impossible to create row vectors of String and Symbol with the transpose operator.
-# This solves this issue, internally in Plots at least.
-
-# commented out on the insistence of the METADATA maintainers
-
-#Base.transpose(x::Symbol) = x
-#Base.transpose(x::String) = x
-
-# ---------------------------------------------------------
-
 import Measures
-
 include("plotmeasures.jl")
-
 using .PlotMeasures
 import .PlotMeasures: Length, AbsoluteLength, Measure, width, height
 # ---------------------------------------------------------
 
-import RecipesPipeline
-import RecipesPipeline:
-    SliceIt,
-    DefaultsDict,
-    Formatted,
-    AbstractSurface,
-    Surface,
-    Volume,
-    is3d,
-    is_surface,
-    needs_3d_axes,
-    group_as_matrix, # for StatsPlots
-    reset_kw!,
-    pop_kw!,
-    scale_func,
-    inverse_scale_func,
-    dateformatter,
-    datetimeformatter,
-    timeformatter
-
-# Use fixed version of Plotly instead of the latest one for stable dependency
-# Ref: https://github.com/JuliaPlots/Plots.jl/pull/2779
-const _plotly_min_js_filename = "plotly-2.6.3.min.js"
-
 include("types.jl")
 include("utils.jl")
 include("colorbars.jl")
 include("axes.jl")
 include("args.jl")
 include("components.jl")
+include("legend.jl")
 include("consts.jl")
 include("themes.jl")
 include("plot.jl")
@@ -206,55 +186,20 @@ include("backends.jl")
 include("output.jl")
 include("ijulia.jl")
 include("fileio.jl")
+
+# Use fixed version of Plotly instead of the latest one for stable dependency
+# Ref: https://github.com/JuliaPlots/Plots.jl/pull/2779
+const _plotly_min_js_filename = "plotly-2.6.3.min.js"
+const CURRENT_BACKEND = CurrentBackend(:none)
+const PLOTS_SEED = 1234
+
 include("init.jl")
-include("legend.jl")
 
 include("backends/plotly.jl")
-include("backends/gr.jl")
 include("backends/web.jl")
+include("backends/gr.jl")
 
-const PlotOrSubplot = Union{Plot,Subplot}
 include("shorthands.jl")
-
-# ---------------------------------------------------------
-
-const CURRENT_BACKEND = CurrentBackend(:none)
-const PLOTS_SEED = 1234
-
-using SnoopPrecompile
-
-@precompile_setup begin
-    n = length(_examples)
-    imports = sizehint!(Expr[], n)
-    examples = sizehint!(Expr[], 10n)
-    for i in setdiff(1:n, _backend_skips[:gr])
-        _examples[i].external && continue
-        (imp = _examples[i].imports) === nothing || push!(imports, imp)
-        func = gensym(string(i))
-        push!(examples, quote
-            $func() = begin  # evaluate each example in a local scope
-                # @show $i  # debug
-                $(_examples[i].exprs)
-                if $i == 1  # only for one example
-                    fn = tempname()
-                    pl = current()
-                    gui(pl)
-                    savefig(pl, "$fn.png")
-                    savefig(pl, "$fn.pdf")
-                end
-                nothing
-            end
-            $func()
-        end)
-    end
-    withenv("GKSwstype" => "nul") do
-        @precompile_all_calls begin
-            eval.(imports)
-            gr()
-            eval.(examples)
-            # eventually eval for another backend ...
-        end
-    end
-end
+include("precompilation.jl")
 
 end
diff --git a/src/animation.jl b/src/animation.jl
index 324a4f3b0..44cc99804 100644
--- a/src/animation.jl
+++ b/src/animation.jl
@@ -1,3 +1,5 @@
+using FFMPEG
+
 "Represents an animation object"
 struct Animation
     dir::String
diff --git a/src/backends/pgfplotsx.jl b/src/backends/pgfplotsx.jl
index e38753c9a..93df2d579 100644
--- a/src/backends/pgfplotsx.jl
+++ b/src/backends/pgfplotsx.jl
@@ -3,9 +3,7 @@ import UUIDs: uuid4
 import Latexify
 import Contour
 
-# FIXME: cannot use `import PGFPlotsX: ...` with @require
-const Options = PGFPlotsX.Options
-const Table = PGFPlotsX.Table
+import .PGFPlotsX: Options, Table
 
 Base.@kwdef mutable struct PGFPlotsXPlot
     is_created::Bool = false
diff --git a/src/init.jl b/src/init.jl
index 07df4c42a..26baddf6e 100644
--- a/src/init.jl
+++ b/src/init.jl
@@ -1,6 +1,6 @@
-using REPL
-using Scratch
 using RelocatableFolders
+using Scratch
+using REPL
 
 const plotly_local_file_path = Ref{Union{Nothing,String}}(nothing)
 const BACKEND_PATH_GASTON = @path joinpath(@__DIR__, "backends", "gaston.jl")
@@ -21,13 +21,30 @@ _plots_defaults() =
         Dict{Symbol,Any}()
     end
 
-function __init__()
+function _plots_theme_defaults()
     user_defaults = _plots_defaults()
     if haskey(user_defaults, :theme)
         theme(pop!(user_defaults, :theme); user_defaults...)
     else
         default(; user_defaults...)
     end
+end
+
+function _plots_plotly_defaults()
+    if get(ENV, "PLOTS_HOST_DEPENDENCY_LOCAL", "false") == "true"
+        global plotly_local_file_path[] =
+            joinpath(@get_scratch!("plotly"), _plotly_min_js_filename)
+        isfile(plotly_local_file_path[]) || Downloads.download(
+            "https://cdn.plot.ly/$(_plotly_min_js_filename)",
+            plotly_local_file_path[],
+        )
+        use_local_plotlyjs[] = true
+    end
+    use_local_dependencies[] = use_local_plotlyjs[]
+end
+
+function __init__()
+    _plots_theme_defaults()
 
     insert!(
         Base.Multimedia.displays,
@@ -77,6 +94,8 @@ function __init__()
         include(BACKEND_PATH_PLOTLYJS)
     end
 
+    _plots_plotly_defaults()
+
     @require PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" begin
         include(BACKEND_PATH_PYPLOT)
     end
@@ -96,19 +115,6 @@ function __init__()
         end
     end
 
-    if get(ENV, "PLOTS_HOST_DEPENDENCY_LOCAL", "false") == "true"
-        global plotly_local_file_path[] =
-            joinpath(@get_scratch!("plotly"), _plotly_min_js_filename)
-        isfile(plotly_local_file_path[]) || Downloads.download(
-            "https://cdn.plot.ly/$(_plotly_min_js_filename)",
-            plotly_local_file_path[],
-        )
-
-        use_local_plotlyjs[] = true
-    end
-
-    use_local_dependencies[] = use_local_plotlyjs[]
-
     @require ImageInTerminal = "d8c32880-2388-543b-8c61-d9f865259254" begin
         if get(ENV, "PLOTS_IMAGE_IN_TERMINAL", "false") == "true" &&
            ImageInTerminal.ENCODER_BACKEND[] == :Sixel
@@ -153,9 +159,14 @@ function __init__()
         # Lists of tuples and GeometryBasics.Points
         # --------------------------------------------------------------------
         @recipe f(v::AVec{<:GeometryBasics.Point}) = RecipesPipeline.unzip(v)
-        @recipe f(p::GeometryBasics.Point) = [p]# Special case for 4-tuples in :ohlc series
+        @recipe f(p::GeometryBasics.Point) = [p]  # Special case for 4-tuples in :ohlc series
         @recipe f(xyuv::AVec{<:Tuple{R1,R2,R3,R4}}) where {R1,R2,R3,R4} =
             get(plotattributes, :seriestype, :path) === :ohlc ?
             OHLC[OHLC(t...) for t in xyuv] : RecipesPipeline.unzip(xyuv)
     end
+
+    @require Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" begin
+        include("unitful.jl")
+        @reexport using .UnitfulRecipes
+    end
 end
diff --git a/src/plotmeasures.jl b/src/plotmeasures.jl
index af7bb16ec..7c9278493 100644
--- a/src/plotmeasures.jl
+++ b/src/plotmeasures.jl
@@ -1,7 +1,9 @@
 module PlotMeasures
-import Measures
-import Measures:
+
+import ..Measures
+import ..Measures:
     Length, AbsoluteLength, Measure, BoundingBox, mm, cm, inch, pt, width, height, w, h
+
 const BBox = Measures.Absolute2DBox
 export BBox, BoundingBox, mm, cm, inch, px, pct, pt, w, h
 
@@ -15,4 +17,5 @@ Base.:*(m1::AbsoluteLength, m2::Length{:pct}) = AbsoluteLength(m1.value * m2.val
 Base.:*(m1::Length{:pct}, m2::AbsoluteLength) = AbsoluteLength(m2.value * m1.value)
 Base.:/(m1::AbsoluteLength, m2::Length{:pct}) = AbsoluteLength(m1.value / m2.value)
 Base.:/(m1::Length{:pct}, m2::AbsoluteLength) = AbsoluteLength(m2.value / m1.value)
+
 end
diff --git a/src/precompilation.jl b/src/precompilation.jl
new file mode 100644
index 000000000..1bba37c73
--- /dev/null
+++ b/src/precompilation.jl
@@ -0,0 +1,37 @@
+using SnoopPrecompile
+
+if get(ENV, "PLOTS_PRECOMPILE", "true") == "true"
+    @precompile_setup begin
+        n = length(_examples)
+        imports = sizehint!(Expr[], n)
+        examples = sizehint!(Expr[], 10n)
+        for i in setdiff(1:n, _backend_skips[:gr])
+            _examples[i].external && continue
+            (imp = _examples[i].imports) === nothing || push!(imports, imp)
+            func = gensym(string(i))
+            push!(examples, quote
+                $func() = begin  # evaluate each example in a local scope
+                    # @show $i  # debug
+                    $(_examples[i].exprs)
+                    if $i == 1  # only for one example
+                        fn = tempname()
+                        pl = current()
+                        gui(pl)
+                        savefig(pl, "$fn.png")
+                        savefig(pl, "$fn.pdf")
+                    end
+                    nothing
+                end
+                $func()
+            end)
+        end
+        withenv("GKSwstype" => "nul") do
+            @precompile_all_calls begin
+                eval.(imports)
+                gr()
+                eval.(examples)
+                # eventually eval for another backend ...
+            end
+        end
+    end
+end
diff --git a/src/types.jl b/src/types.jl
index f6b5cbb07..80b7b00e1 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -109,6 +109,7 @@ mutable struct Plot{T<:AbstractBackend} <: AbstractPlot{T}
 end
 
 struct PlaceHolder end
+const PlotOrSubplot = Union{Plot,Subplot}
 
 # -----------------------------------------------------------
 
diff --git a/src/unitful.jl b/src/unitful.jl
new file mode 100644
index 000000000..f9ac78b77
--- /dev/null
+++ b/src/unitful.jl
@@ -0,0 +1,293 @@
+# previously https://github.com/jw3126/UnitfulRecipes.jl
+# authors: Benoit Pasquier (@briochemc) - David Gustavsson (@gustaphe) - Jan Weidner (@jw3126)
+
+module UnitfulRecipes
+
+using ..Unitful: Quantity, unit, ustrip, Unitful, dimension, Units
+using ..RecipesBase
+export @P_str
+
+import ..locate_annotation, ..PlotText, ..Subplot, ..AVec, ..AMat
+
+const MissingOrQuantity = Union{Missing,<:Quantity}
+
+#==========
+Main recipe
+==========#
+
+@recipe function f(::Type{T}, x::T) where {T<:AbstractArray{<:MissingOrQuantity}}
+    axisletter = plotattributes[:letter]   # x, y, or z
+    clims_types = (:contour, :contourf, :heatmap, :surface)
+    if axisletter === :z && get(plotattributes, :seriestype, :nothing) ∈ clims_types
+        u = get(plotattributes, :zunit, unit(eltype(x)))
+        ustripattribute!(plotattributes, :clims, u)
+        append_unit_if_needed!(plotattributes, :colorbar_title, u)
+    end
+    fixaxis!(plotattributes, x, axisletter)
+end
+
+function fixaxis!(attr, x, axisletter)
+    # Attribute keys
+    axislabel = Symbol(axisletter, :guide) # xguide, yguide, zguide
+    axislims = Symbol(axisletter, :lims)   # xlims, ylims, zlims
+    axisticks = Symbol(axisletter, :ticks) # xticks, yticks, zticks
+    err = Symbol(axisletter, :error)       # xerror, yerror, zerror
+    axisunit = Symbol(axisletter, :unit)   # xunit, yunit, zunit
+    axis = Symbol(axisletter, :axis)       # xaxis, yaxis, zaxis
+    # Get the unit
+    u = pop!(attr, axisunit, unit(eltype(x)))
+    # If the subplot already exists with data, get its unit
+    sp = get(attr, :subplot, 1)
+    if sp ≤ length(attr[:plot_object]) && attr[:plot_object].n > 0
+        label = attr[:plot_object][sp][axis][:guide]
+        if label isa UnitfulString
+            u = label.unit
+        end
+        # If label was not given as an argument, reuse
+        get!(attr, axislabel, label)
+    end
+    # Fix the attributes: labels, lims, ticks, marker/line stuff, etc.
+    append_unit_if_needed!(attr, axislabel, u)
+    ustripattribute!(attr, axislims, u)
+    ustripattribute!(attr, axisticks, u)
+    ustripattribute!(attr, err, u)
+    if axisletter === :y
+        ustripattribute!(attr, :ribbon, u)
+        ustripattribute!(attr, :fillrange, u)
+    end
+    fixaspectratio!(attr, u, axisletter)
+    fixmarkercolor!(attr)
+    fixmarkersize!(attr)
+    fixlinecolor!(attr)
+    # Strip the unit
+    ustrip.(u, x)
+end
+
+# Recipe for (x::AVec, y::AVec, z::Surface) types
+@recipe function f(x::AVec, y::AVec, z::AMat{T}) where {T<:Quantity}
+    u = get(plotattributes, :zunit, unit(eltype(z)))
+    ustripattribute!(plotattributes, :clims, u)
+    z = fixaxis!(plotattributes, z, :z)
+    append_unit_if_needed!(plotattributes, :colorbar_title, u)
+    x, y, z
+end
+
+# Recipe for vectors of vectors
+@recipe function f(::Type{T}, x::T) where {T<:AVec{<:AVec{<:MissingOrQuantity}}}
+    axisletter = plotattributes[:letter]   # x, y, or z
+    map(x -> fixaxis!(plotattributes, x, axisletter), x)
+end
+
+# Recipe for bare units
+@recipe function f(::Type{T}, x::T) where {T<:Units}
+    primary := false
+    Float64[] * x
+end
+
+# Recipes for functions
+@recipe f(f::Function, x::T) where {T<:AVec{<:MissingOrQuantity}} = x, f.(x)
+@recipe f(x::T, f::Function) where {T<:AVec{<:MissingOrQuantity}} = x, f.(x)
+@recipe f(x::T, y::AVec, f::Function) where {T<:AVec{<:MissingOrQuantity}} = x, y, f.(x', y)
+@recipe f(x::AVec, y::T, f::Function) where {T<:AVec{<:MissingOrQuantity}} = x, y, f.(x', y)
+@recipe function f(
+    x::T1,
+    y::T2,
+    f::Function,
+) where {T1<:AVec{<:MissingOrQuantity},T2<:AVec{<:MissingOrQuantity}}
+    x, y, f.(x', y)
+end
+@recipe function f(f::Function, u::Units)
+    uf = UnitFunction(f, [u])
+    recipedata = RecipesBase.apply_recipe(plotattributes, uf)
+    _, xmin, xmax = recipedata[1].args
+    return f, xmin * u, xmax * u
+end
+
+"""
+```julia
+UnitFunction
+```
+A function, bundled with the assumed units of each of its inputs.
+
+```julia
+f(x, y) = x^2 + y
+uf = UnitFunction(f, u"m", u"m^2")
+uf(3, 2) == f(3u"m", 2u"m"^2) == 7u"m^2"
+```
+"""
+struct UnitFunction <: Function
+    f::Function
+    u::Vector{Units}
+end
+(f::UnitFunction)(args...) = f.f((args .* f.u)...)
+
+#===============
+Attribute fixing
+===============#
+# Aspect ratio
+function fixaspectratio!(attr, u, axisletter)
+    aspect_ratio = get!(attr, :aspect_ratio, :auto)
+    if aspect_ratio in (:auto, :none)
+        # Keep the default behavior (let Plots figure it out)
+        return
+    end
+    if aspect_ratio === :equal
+        aspect_ratio = 1
+    end
+    #=======================================================================================
+    Implementation example:
+
+    Consider an x axis in `u"m"` and a y axis in `u"s"`, and an `aspect_ratio` in `u"m/s"`.
+    On the first pass, `axisletter` is `:x`, so `aspect_ratio` is converted to `u"m/s"/u"m"
+    = u"s^-1"`. On the second pass, `axisletter` is `:y`, so `aspect_ratio` becomes
+    `u"s^-1"*u"s" = 1`. If at this point `aspect_ratio` is *not* unitless, an error has been
+    made, and the default aspect ratio fixing of Plots throws a `DimensionError` as it tries
+    to compare `0 < 1u"m"`.
+    =======================================================================================#
+    if axisletter === :y
+        attr[:aspect_ratio] = aspect_ratio * u
+        return
+    end
+    if axisletter === :x
+        attr[:aspect_ratio] = aspect_ratio / u
+        return
+    end
+    return
+end
+
+# Markers / lines
+function fixmarkercolor!(attr)
+    u = ustripattribute!(attr, :marker_z)
+    ustripattribute!(attr, :clims, u)
+    u == Unitful.NoUnits || append_unit_if_needed!(attr, :colorbar_title, u)
+end
+fixmarkersize!(attr) = ustripattribute!(attr, :markersize)
+fixlinecolor!(attr) = ustripattribute!(attr, :line_z)
+
+# strip unit from attribute[key]
+ustripattribute!(attr, key) =
+    if haskey(attr, key)
+        v = attr[key]
+        u = unit(eltype(v))
+        attr[key] = ustrip.(u, v)
+        return u
+    else
+        return Unitful.NoUnits
+    end
+# If supplied, use the unit (optional 3rd argument)
+function ustripattribute!(attr, key, u)
+    if haskey(attr, key)
+        v = attr[key]
+        if eltype(v) <: Quantity
+            attr[key] = ustrip.(u, v)
+        end
+    end
+    u
+end
+
+#=======================================
+Label string containing unit information
+=======================================#
+
+abstract type AbstractProtectedString <: AbstractString end
+struct ProtectedString{S} <: AbstractProtectedString
+    content::S
+end
+struct UnitfulString{S,U} <: AbstractProtectedString
+    content::S
+    unit::U
+end
+# Minimum required AbstractString interface to work with Plots
+const S = AbstractProtectedString
+Base.iterate(n::S) = iterate(n.content)
+Base.iterate(n::S, i::Integer) = iterate(n.content, i)
+Base.codeunit(n::S) = codeunit(n.content)
+Base.ncodeunits(n::S) = ncodeunits(n.content)
+Base.isvalid(n::S, i::Integer) = isvalid(n.content, i)
+Base.pointer(n::S) = pointer(n.content)
+Base.pointer(n::S, i::Integer) = pointer(n.content, i)
+"""
+    P_str(s)
+
+Creates a string that will be Protected from recipe passes.
+
+Example:
+```julia
+julia> plot([0,1]u"m", [1,2]u"m/s^2", xlabel=P"This label will NOT display units")
+
+julia> plot([0,1]u"m", [1,2]u"m/s^2", xlabel="This label will display units")
+```
+"""
+macro P_str(s)
+    return ProtectedString(s)
+end
+
+#=====================================
+Append unit to labels when appropriate
+=====================================#
+
+function append_unit_if_needed!(attr, key, u::Unitful.Units)
+    label = get(attr, key, nothing)
+    append_unit_if_needed!(attr, key, label, u)
+end
+# dispatch on the type of `label`
+append_unit_if_needed!(attr, key, label::ProtectedString, u) = nothing
+append_unit_if_needed!(attr, key, label::UnitfulString, u) = nothing
+function append_unit_if_needed!(attr, key, label::Nothing, u)
+    attr[key] = UnitfulString(string(u), u)
+end
+function append_unit_if_needed!(attr, key, label::S, u) where {S<:AbstractString}
+    if !isempty(label)
+        attr[key] =
+            UnitfulString(S(format_unit_label(label, u, get(attr, :unitformat, :round))), u)
+    end
+end
+
+#=============================================
+Surround unit string with specified delimiters
+=============================================#
+format_unit_label(l, u, f::Nothing)                    = string(l, ' ', u)
+format_unit_label(l, u, f::Function)                   = f(l, u)
+format_unit_label(l, u, f::AbstractString)             = string(l, f, u)
+format_unit_label(l, u, f::NTuple{2,<:AbstractString}) = string(l, f[1], u, f[2])
+format_unit_label(l, u, f::NTuple{3,<:AbstractString}) = string(f[1], l, f[2], u, f[3])
+format_unit_label(l, u, f::Char)                       = string(l, ' ', f, ' ', u)
+format_unit_label(l, u, f::NTuple{2,Char})             = string(l, ' ', f[1], u, f[2])
+format_unit_label(l, u, f::NTuple{3,Char})             = string(f[1], l, ' ', f[2], u, f[3])
+format_unit_label(l, u, f::Bool)                       = f ? format_unit_label(l, u, :round) : format_unit_label(l, u, nothing)
+
+const UNIT_FORMATS = Dict(
+    :round => ('(', ')'),
+    :square => ('[', ']'),
+    :curly => ('{', '}'),
+    :angle => ('<', '>'),
+    :slash => '/',
+    :slashround => (" / (", ")"),
+    :slashsquare => (" / [", "]"),
+    :slashcurly => (" / {", "}"),
+    :slashangle => (" / <", ">"),
+    :verbose => " in units of ",
+)
+
+format_unit_label(l, u, f::Symbol) = format_unit_label(l, u, UNIT_FORMATS[f])
+
+#==============
+Fix annotations
+===============#
+locate_annotation(
+    sp::Subplot,
+    x::MissingOrQuantity,
+    y::MissingOrQuantity,
+    label::PlotText,
+) = (ustrip(x), ustrip(y), label)
+locate_annotation(
+    sp::Subplot,
+    x::MissingOrQuantity,
+    y::MissingOrQuantity,
+    z::MissingOrQuantity,
+    label::PlotText,
+) = (ustrip(x), ustrip(y), ustrip(z), label)
+locate_annotation(sp::Subplot, rel::NTuple{N,<:MissingOrQuantity}, label) where {N} =
+    locate_annotation(sp, ustrip.(rel), label)
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index be6bd8c8f..9f7992769 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,3 +1,4 @@
+import Unitful: m, s, cm, DimensionError
 import Plots: PLOTS_SEED, Plot, with
 import GeometryBasics
 import ImageMagick
@@ -10,6 +11,7 @@ using FilePathsBase
 using LaTeXStrings
 using RecipesBase
 using TestImages
+using Unitful
 using FileIO
 using Plots
 using Dates
@@ -36,6 +38,7 @@ for name in (
     "components",
     "shorthands",
     "recipes",
+    "unitful",
     "hdf5plots",
     "pgfplotsx",
     "plotly",
diff --git a/test/test_defaults.jl b/test/test_defaults.jl
index bd7cbc8cb..d96f82152 100644
--- a/test/test_defaults.jl
+++ b/test/test_defaults.jl
@@ -1,8 +1,5 @@
-using Plots, Test
-using Plots.Colors
-
 const PLOTS_DEFAULTS = Dict(:theme => :wong2, :fontfamily => :palantino)
-Plots.__init__()
+Plots._plots_theme_defaults()
 
 @testset "Loading theme" begin
     pl = plot(1:5)
@@ -11,7 +8,7 @@ Plots.__init__()
 end
 
 empty!(PLOTS_DEFAULTS)
-Plots.__init__()
+Plots._plots_theme_defaults()
 
 @testset "default" begin
     default(fillrange = 0)
diff --git a/test/test_misc.jl b/test/test_misc.jl
index 19ed05f91..878d9801b 100644
--- a/test/test_misc.jl
+++ b/test/test_misc.jl
@@ -11,7 +11,7 @@ end
     @test Plots.plotly_local_file_path[] === nothing
     temp = Plots.use_local_dependencies[]
     withenv("PLOTS_HOST_DEPENDENCY_LOCAL" => true) do
-        Plots.__init__()
+        Plots._plots_plotly_defaults()
         @test Plots.plotly_local_file_path[] isa String
         @test isfile(Plots.plotly_local_file_path[])
         @test Plots.use_local_dependencies[] = true
diff --git a/test/test_recipes.jl b/test/test_recipes.jl
index 7c03e29c2..7ba088ce3 100644
--- a/test/test_recipes.jl
+++ b/test/test_recipes.jl
@@ -43,14 +43,17 @@ end
     @test length(sticks) == 1
 end
 
+# NOTE: the following test seems to trigger these deprecated warnings:
+# WARNING: importing deprecated binding Colors.RGB1 into PlotUtils.
+# WARNING: importing deprecated binding Colors.RGB1 into Plots.
 @testset "framestyle axes" begin
     pl = plot(-1:1, -1:1, -1:1)
     sp = pl.subplots[1]
     defaultret = Plots.axis_drawing_info_3d(sp, :x)
-    for letter in [:x, :y, :z]
-        for fr in [:box :semi :origin :zerolines :grid :none]
+    for letter in (:x, :y, :z)
+        for framestyle in [:box :semi :origin :zerolines :grid :none]
             prevha = UInt64(0)
-            push!(sp.attr, :framestyle => fr)
+            push!(sp.attr, :framestyle => framestyle)
             ret = Plots.axis_drawing_info_3d(sp, letter)
             ha = hash(string(ret))
             @test ha != prevha
diff --git a/test/test_unitful.jl b/test/test_unitful.jl
new file mode 100644
index 000000000..50e94b2e0
--- /dev/null
+++ b/test/test_unitful.jl
@@ -0,0 +1,341 @@
+
+# Some helper functions to access the subplot labels and the series inside each test plot
+xguide(pl, idx = length(pl.subplots)) = pl.subplots[idx].attr[:xaxis].plotattributes[:guide]
+yguide(pl, idx = length(pl.subplots)) = pl.subplots[idx].attr[:yaxis].plotattributes[:guide]
+zguide(pl, idx = length(pl.subplots)) = pl.subplots[idx].attr[:zaxis].plotattributes[:guide]
+xseries(pl, idx = length(pl.series_list)) = pl.series_list[idx].plotattributes[:x]
+yseries(pl, idx = length(pl.series_list)) = pl.series_list[idx].plotattributes[:y]
+zseries(pl, idx = length(pl.series_list)) = pl.series_list[idx].plotattributes[:z]
+
+macro isplot(ex) # @isplot macro to streamline tests
+    :(@test $(esc(ex)) isa Plot)
+end
+
+@testset "heatmap" begin
+    x = (1:3)m
+    @isplot heatmap(x * x', clims = (1, 7)) # unitless
+    @isplot heatmap(x * x', clims = (2m^2, 8m^2)) # units
+    @isplot heatmap(x * x', clims = (2e6u"mm^2", 7e-6u"km^2")) # conversion
+    @isplot heatmap(1:3, (1:3)m, x * x', clims = (1m^2, 7e-6u"km^2")) # mixed
+end
+
+@testset "plot(y)" begin
+    y = rand(3)m
+
+    @testset "no keyword argument" begin
+        @test yguide(plot(y)) == "m"
+        @test yseries(plot(y)) ≈ ustrip.(y)
+    end
+
+    @testset "ylabel" begin
+        @test yguide(plot(y, ylabel = "hello")) == "hello (m)"
+        @test yguide(plot(y, ylabel = P"hello")) == "hello"
+        @test yguide(plot(y, ylabel = "")) == ""
+        pl = plot(y; ylabel = "hello")
+        plot!(pl, -y)
+        @test yguide(pl) == "hello (m)"
+    end
+
+    @testset "yunit" begin
+        @test yguide(plot(y, yunit = cm)) == "cm"
+        @test yseries(plot(y, yunit = cm)) ≈ ustrip.(cm, y)
+    end
+
+    @testset "ylims" begin # Using all(lims .≈ lims) because of uncontrolled type conversions?
+        @test all(ylims(plot(y, ylims = (-1, 3))) .≈ (-1, 3))
+        @test all(ylims(plot(y, ylims = (-1m, 3m))) .≈ (-1, 3))
+        @test all(ylims(plot(y, ylims = (-100cm, 300cm))) .≈ (-1, 3))
+        @test all(ylims(plot(y, ylims = (-100cm, 3m))) .≈ (-1, 3))
+    end
+
+    @testset "keyword combinations" begin
+        @test yguide(plot(y, yunit = cm, ylabel = "hello")) == "hello (cm)"
+        @test yseries(plot(y, yunit = cm, ylabel = "hello")) ≈ ustrip.(cm, y)
+        @test all(ylims(plot(y, yunit = cm, ylims = (-1, 3))) .≈ (-1, 3))
+        @test all(ylims(plot(y, yunit = cm, ylims = (-1, 3))) .≈ (-1, 3))
+        @test all(ylims(plot(y, yunit = cm, ylims = (-100cm, 300cm))) .≈ (-100, 300))
+        @test all(ylims(plot(y, yunit = cm, ylims = (-100cm, 3m))) .≈ (-100, 300))
+    end
+end
+
+@testset "plot(x,y)" begin
+    x, y = randn(3)m, randn(3)s
+
+    @testset "no keyword argument" begin
+        @test xguide(plot(x, y)) == "m"
+        @test xseries(plot(x, y)) ≈ ustrip.(x)
+        @test yguide(plot(x, y)) == "s"
+        @test yseries(plot(x, y)) ≈ ustrip.(y)
+    end
+
+    @testset "labels" begin
+        @test xguide(plot(x, y, xlabel = "hello")) == "hello (m)"
+        @test xguide(plot(x, y, xlabel = P"hello")) == "hello"
+        @test yguide(plot(x, y, ylabel = "hello")) == "hello (s)"
+        @test yguide(plot(x, y, ylabel = P"hello")) == "hello"
+        @test xguide(plot(x, y, xlabel = "hello", ylabel = "hello")) == "hello (m)"
+        @test xguide(plot(x, y, xlabel = P"hello", ylabel = P"hello")) == "hello"
+        @test yguide(plot(x, y, xlabel = "hello", ylabel = "hello")) == "hello (s)"
+        @test yguide(plot(x, y, xlabel = P"hello", ylabel = P"hello")) == "hello"
+    end
+
+    @testset "unitformat" begin
+        args = (x, y)
+        kwargs = (:xlabel => "hello", :ylabel => "hello")
+        @test yguide(plot(args...; kwargs..., unitformat = nothing)) == "hello s"
+        @test yguide(
+            plot(
+                args...;
+                kwargs...,
+                unitformat = (l, u) -> string(u, " is the unit of ", l),
+            ),
+        ) == "s is the unit of hello"
+        @test yguide(plot(args...; kwargs..., unitformat = ", dear ")) == "hello, dear s"
+        @test yguide(plot(args...; kwargs..., unitformat = (", dear ", " esq."))) ==
+              "hello, dear s esq."
+        @test yguide(
+            plot(args...; kwargs..., unitformat = ("well ", ", dear ", " esq.")),
+        ) == "well hello, dear s esq."
+        @test yguide(plot(args...; kwargs..., unitformat = '?')) == "hello ? s"
+        @test yguide(plot(args...; kwargs..., unitformat = ('<', '>'))) == "hello <s>"
+        @test yguide(plot(args...; kwargs..., unitformat = ('A', 'B', 'C'))) == "Ahello BsC"
+        @test yguide(plot(args...; kwargs..., unitformat = false)) == "hello s"
+        @test yguide(plot(args...; kwargs..., unitformat = true)) == "hello (s)"
+        @test yguide(plot(args...; kwargs..., unitformat = :round)) == "hello (s)"
+        @test yguide(plot(args...; kwargs..., unitformat = :square)) == "hello [s]"
+        @test yguide(plot(args...; kwargs..., unitformat = :curly)) == "hello {s}"
+        @test yguide(plot(args...; kwargs..., unitformat = :angle)) == "hello <s>"
+        @test yguide(plot(args...; kwargs..., unitformat = :slash)) == "hello / s"
+        @test yguide(plot(args...; kwargs..., unitformat = :slashround)) == "hello / (s)"
+        @test yguide(plot(args...; kwargs..., unitformat = :slashsquare)) == "hello / [s]"
+        @test yguide(plot(args...; kwargs..., unitformat = :slashcurly)) == "hello / {s}"
+        @test yguide(plot(args...; kwargs..., unitformat = :slashangle)) == "hello / <s>"
+        @test yguide(plot(args...; kwargs..., unitformat = :verbose)) ==
+              "hello in units of s"
+    end
+end
+
+@testset "With functions" begin
+    x, y = randn(3), randn(3)
+    @testset "plot(f, x) / plot(x, f)" begin
+        f(x) = x^2
+        @test plot(f, x * m) isa Plot
+        @test plot(x * m, f) isa Plot
+        g(x) = x * m # If the unit comes from the function only then it throws
+        @test_throws DimensionError plot(x, g) isa Plot
+        @test_throws DimensionError plot(g, x) isa Plot
+    end
+    @testset "plot(x, y, f)" begin
+        f(x, y) = x * y
+        @test plot(x * m, y * s, f) isa Plot
+        @test plot(x * m, y, f) isa Plot
+        @test plot(x, y * s, f) isa Plot
+        g(x, y) = x * y * m # If the unit comes from the function only then it throws
+        @test_throws DimensionError plot(x, y, g) isa Plot
+    end
+    @testset "plot(f, u)" begin
+        f(x) = x^2
+        pl = plot(x * m, f.(x * m))
+        @test plot!(pl, f, m) isa Plot
+        @test_throws DimensionError plot!(pl, f, s) isa Plot
+        pl = plot(f, m)
+        @test xguide(pl) == string(m)
+        @test yguide(pl) == string(m^2)
+        g(x) = exp(x / (3m))
+        @test plot(g, u"m") isa Plot
+    end
+end
+
+@testset "More plots" begin
+    @testset "data as $dtype" for dtype in
+                                  [:Vectors, :Matrices, Symbol("Vectors of vectors")]
+        if dtype == :Vectors
+            x, y, z = randn(10), randn(10), randn(10)
+        elseif dtype == :Matrices
+            x, y, z = randn(10, 2), randn(10, 2), randn(10, 2)
+        else
+            x, y, z = [rand(10), rand(20)], [rand(10), rand(20)], [rand(10), rand(20)]
+        end
+
+        @testset "One array" begin
+            @test plot(x * m) isa Plot
+            @test plot(x * m, ylabel = "x") isa Plot
+            @test plot(x * m, ylims = (-1, 1)) isa Plot
+            @test plot(x * m, ylims = (-1, 1) .* m) isa Plot
+            @test plot(x * m, yunit = u"km") isa Plot
+            @test plot(x * m, xticks = (1:3) * m) isa Plot
+        end
+
+        @testset "Two arrays" begin
+            @test plot(x * m, y * s) isa Plot
+            @test plot(x * m, y * s, xlabel = "x") isa Plot
+            @test plot(x * m, y * s, xlims = (-1, 1)) isa Plot
+            @test plot(x * m, y * s, xlims = (-1, 1) .* m) isa Plot
+            @test plot(x * m, y * s, xunit = u"km") isa Plot
+            @test plot(x * m, y * s, ylabel = "y") isa Plot
+            @test plot(x * m, y * s, ylims = (-1, 1)) isa Plot
+            @test plot(x * m, y * s, ylims = (-1, 1) .* s) isa Plot
+            @test plot(x * m, y * s, yunit = u"ks") isa Plot
+            @test plot(x * m, y * s, yticks = (1:3) * s) isa Plot
+            @test scatter(x * m, y * s) isa Plot
+            if dtype ≠ Symbol("Vectors of vectors")
+                @test scatter(x * m, y * s, zcolor = z * (m / s)) isa Plot
+            end
+        end
+
+        @testset "Three arrays" begin
+            @test plot(x * m, y * s, z * (m / s)) isa Plot
+            @test plot(x * m, y * s, z * (m / s), xlabel = "x") isa Plot
+            @test plot(x * m, y * s, z * (m / s), xlims = (-1, 1)) isa Plot
+            @test plot(x * m, y * s, z * (m / s), xlims = (-1, 1) .* m) isa Plot
+            @test plot(x * m, y * s, z * (m / s), xunit = u"km") isa Plot
+            @test plot(x * m, y * s, z * (m / s), ylabel = "y") isa Plot
+            @test plot(x * m, y * s, z * (m / s), ylims = (-1, 1)) isa Plot
+            @test plot(x * m, y * s, z * (m / s), ylims = (-1, 1) .* s) isa Plot
+            @test plot(x * m, y * s, z * (m / s), yunit = u"ks") isa Plot
+            @test plot(x * m, y * s, z * (m / s), zlabel = "z") isa Plot
+            @test plot(x * m, y * s, z * (m / s), zlims = (-1, 1)) isa Plot
+            @test plot(x * m, y * s, z * (m / s), zlims = (-1, 1) .* (m / s)) isa Plot
+            @test plot(x * m, y * s, z * (m / s), zunit = u"km/hr") isa Plot
+            @test plot(x * m, y * s, z * (m / s), zticks = (1:2) * m / s) isa Plot
+            @test scatter(x * m, y * s, z * (m / s)) isa Plot
+        end
+
+        @testset "Unitful/unitless combinations" begin
+            mystr(x::Array{<:Quantity}) = "Q"
+            mystr(x::Array) = "A"
+            @testset "plot($(mystr(xs)), $(mystr(ys)))" for xs in [x, x * m],
+                ys in [y, y * s]
+
+                @test plot(xs, ys) isa Plot
+            end
+            @testset "plot($(mystr(xs)), $(mystr(ys)), $(mystr(zs)))" for xs in [x, x * m],
+                ys in [y, y * s],
+                zs in [z, z * (m / s)]
+
+                @test plot(xs, ys, zs) isa Plot
+            end
+        end
+    end
+
+    @testset "scatter(x::$(us[1]), y::$(us[2]))" for us in collect(
+        Iterators.product(fill([1, u"m", u"s"], 2)...),
+    )
+        x, y = rand(10) * us[1], rand(10) * us[2]
+        @test scatter(x, y) isa Plot
+        @test scatter(x, y, markersize = x) isa Plot
+        @test scatter(x, y, line_z = x) isa Plot
+    end
+
+    @testset "contour(x::$(us[1]), y::$(us[2]))" for us in collect(
+        Iterators.product(fill([1, u"m", u"s"], 2)...),
+    )
+        x, y = (1:0.01:2) * us[1], (1:0.02:2) * us[2]
+        z = x' ./ y
+        @test contour(x, y, z) isa Plot
+        @test contourf(x, y, z) isa Plot
+    end
+
+    @testset "ProtectedString" begin
+        y = rand(10) * u"m"
+        @test plot(y, label = P"meters") isa Plot
+    end
+end
+
+@testset "Comparing apples and oranges" begin
+    x1 = rand(10) * u"m"
+    x2 = rand(10) * u"cm"
+    x3 = rand(10) * u"s"
+    pl = plot(x1)
+    pl = plot!(pl, x2)
+    @test yguide(pl) == "m"
+    @test yseries(pl) ≈ ustrip.(x2) / 100
+    @test_throws DimensionError plot!(pl, x3) # can't place seconds on top of meters!
+end
+
+@testset "Bare units" begin
+    pl = plot(u"m", u"s")
+    @test xguide(pl) == "m"
+    @test yguide(pl) == "s"
+    @test iszero(length(pl.series_list[1].plotattributes[:y]))
+    hline!(pl, [1u"hr"])
+    @test yguide(pl) == "s"
+end
+
+@testset "Inset subplots" begin
+    x1 = rand(10) * u"m"
+    x2 = rand(10) * u"s"
+    pl = plot(x1)
+    pl = plot!(x2, inset = bbox(0.5, 0.5, 0.3, 0.3), subplot = 2)
+    @test yguide(pl, 1) == "m"
+    @test yguide(pl, 2) == "s"
+end
+
+@testset "Missing values" begin
+    x = 1:5
+    y = [1.0 * u"s", 2.0 * u"s", missing, missing, missing]
+    pl = plot(x, y)
+    @test yguide(pl, 1) == "s"
+end
+
+@testset "Errors" begin
+    x = rand(10) * u"mm"
+    ex = rand(10) * u"μm"
+    y = rand(10) * u"s"
+    ey = rand(10) * u"ms"
+    pl = plot(x, y, xerr = ex, yerr = ey)
+    @test pl isa Plot
+    @test xguide(pl) == "mm"
+    @test yguide(pl) == "s"
+end
+
+@testset "Ribbon" begin
+    x = rand(10) * u"mm"
+    y = rand(10) * u"s"
+    ribbon = rand(10) * u"ms"
+    pl = plot(x, y, ribbon = ribbon)
+    @test pl isa Plot
+    @test xguide(pl) == "mm"
+    @test yguide(pl) == "s"
+end
+
+@testset "Fillrange" begin
+    x = rand(10) * u"mm"
+    y = rand(10) * u"s"
+    fillrange = rand(10) * u"ms"
+    pl = plot(x, y, fillrange = fillrange)
+    @test pl isa Plot
+    @test xguide(pl) == "mm"
+    @test yguide(pl) == "s"
+end
+
+@testset "Aspect ratio" begin
+    testfile = tempname() * ".png"
+    pl = plot((1:10)u"m", (1:10)u"dm"; aspect_ratio = :equal)
+    savefig(pl, testfile) # Force a render, to make it evaluate aspect ratio
+    @test abs(-(ylims(pl)...)) > 50
+    pl = plot((1:10)u"m", (1:10)u"dm"; aspect_ratio = 2)
+    savefig(pl, testfile)
+    @test 25 < abs(-(ylims(pl)...)) < 50
+    pl = plot((1:10)u"m", (1:10)u"s"; aspect_ratio = 1u"m/s")
+    savefig(pl, testfile)
+    @test 7.5 < abs(-(ylims(pl)...)) < 12.5
+    @test_throws DimensionError savefig(
+        plot((1:10)u"m", (1:10)u"s"; aspect_ratio = :equal),
+        testfile,
+    )
+end
+
+# https://github.com/jw3126/UnitfulRecipes.jl/issues/60
+@testset "Start with empty plot" begin
+    pl = plot()
+    plot!(pl, (1:3)m)
+    @test yguide(pl) == "m"
+end
+
+# https://github.com/jw3126/UnitfulRecipes.jl/issues/79
+@testset "Annotate" begin
+    pl = plot([0, 1]u"s", [0, 1]u"m")
+    annotate!(pl, [0.25]u"s", [0.5]u"m", text("annotation"))
+    @test show(devnull, pl) isa Nothing
+end