diff --git a/docs/make.jl b/docs/make.jl
index 1c64f13e..104b0426 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -8,6 +8,7 @@ pages=["Home" => "index.md",
         "Interpolation algorithms" => "control.md",
         "Extrapolation" => "extrapolation.md",
         "Convenience Constructors" => "convenience-construction.md",
+        "Knot Iteration" => "iterate.md",
         "Developer documentation" => "devdocs.md",
         "Library" => "api.md"]
 )
diff --git a/docs/src/iterate.md b/docs/src/iterate.md
new file mode 100644
index 00000000..caa2efd5
--- /dev/null
+++ b/docs/src/iterate.md
@@ -0,0 +1,256 @@
+# Knot Iteration
+
+Given an `AbstractInterpolation` `itp` get an iterator over it's knots using
+`knots(itp)`
+
+```julia
+using Interpolations
+itp = interpolate(rand(4), options...)
+kiter = knots(itp) # Iterator over knots
+collect(kiter) # Array of knots [1, 2, 3, 4]
+
+```
+
+For multiple dimensions, the iterator will return tuples of positions
+(ie. `(x, y, ...)`),  with the first coordinate changing the fastest.
+
+```jldoctest iterate-interpolate; setup = :(using Interpolations)
+julia> itp = interpolate(ones(3,3), BSpline(Linear()));
+
+julia> kiter = knots(itp);
+
+julia> collect(kiter)
+3×3 Array{Tuple{Int64,Int64},2}:
+ (1, 1)  (1, 2)  (1, 3)
+ (2, 1)  (2, 2)  (2, 3)
+ (3, 1)  (3, 2)  (3, 3)
+```
+
+The number of elements and size of the iterator can be found as shown:
+
+```jldoctest iterate-interpolate; setup = :(using Interpolations)
+julia> length(kiter)
+9
+
+julia> size(kiter)
+(3, 3)
+
+```
+
+
+## Extrapolated Knots
+
+Given an `AbstractExtrapolation` `etp`, `knots(etp)` will also iterate over the
+the knots with the following behavior.
+
+- For `Throw`, `Flat`, `Line` iterate the knots once
+- For `Periodic` and `Reflect` generate an infinite sequence of knots starting
+  at the first knot.
+
+As `Periodic` and `Reflect` generate infinite sequences of knots, `length` and
+`size` are undefined. For `Throw`, `Flat`, `Line`, `length` and `size` behave as
+expected.
+
+### Periodic
+
+With Periodic boundary condition, knots repeat indefinitely with the first and
+last knot being co-located. (ie. in the example below `etp(2.0) = 1.0` not
+`4.0`).
+
+```jldoctest periodic-demo; setup = :(using Interpolations)
+julia> x = [1.0, 1.5, 1.75, 2.0];
+
+julia> etp = LinearInterpolation(x, x.^2, extrapolation_bc=Periodic());
+
+julia> kiter = knots(etp);
+
+julia> k = Iterators.take(kiter, 6) |> collect
+6-element Array{Float64,1}:
+ 1.0
+ 1.5
+ 1.75
+ 2.0
+ 2.5
+ 2.75
+
+```
+
+Extrapolating to the generated knots `etp.(k)`, confirms that the extrapolated
+knots do map back to the correct inbound knots (ie. `etp(k[1]) == etp(k[4])`).
+
+```jldoctest periodic-demo
+julia> etp.(k)
+6-element Array{Float64,1}:
+ 1.0
+ 2.25
+ 3.0625
+ 1.0
+ 2.25
+ 3.0625
+
+```
+
+### Reflect
+
+With the `Reflect` boundary condition knots repeat indefinitely, following the
+pattern shown below (Offset terms are not shown for brevity).
+
+```
+k[1], k[2], ..., k[end-1], k[end], k[end+1], ... k[2], k[1], k[2], ...
+```
+
+```jldoctest reflect-demo; setup = :(using Interpolations)
+julia> x = [1.0, 1.5, 1.75, 2.0];
+
+julia> etp = LinearInterpolation(x, x.^2, extrapolation_bc=Reflect());
+
+julia> kiter = knots(etp);
+
+julia> k = Iterators.take(kiter, 6) |> collect
+6-element Array{Float64,1}:
+ 1.0
+ 1.5
+ 1.75
+ 2.0
+ 2.25
+ 2.5
+
+```
+
+Evaluating the extrapolation at `etp.(k)` confirms that the extrapolated knots
+correspond to the correct inbound knots.
+
+```jldoctest reflect-demo
+julia> etp.(k)
+6-element Array{Float64,1}:
+ 1.0
+ 2.25
+ 3.0625
+ 4.0
+ 3.0625
+ 2.25
+
+```
+
+### Multiple Dimensions
+
+As with an `AbstractInterpolation`, iterating over knots for a
+multi-dimensional extrapolation also supported.
+
+```jldoctest; setup = :(using Interpolations)
+julia> x = [1.0, 1.5, 1.75, 2.0];
+
+julia> etp = LinearInterpolation((x, x), x.*x');
+
+julia> knots(etp) |> collect
+4×4 Array{Tuple{Float64,Float64},2}:
+ (1.0, 1.0)   (1.0, 1.5)   (1.0, 1.75)   (1.0, 2.0)
+ (1.5, 1.0)   (1.5, 1.5)   (1.5, 1.75)   (1.5, 2.0)
+ (1.75, 1.0)  (1.75, 1.5)  (1.75, 1.75)  (1.75, 2.0)
+ (2.0, 1.0)   (2.0, 1.5)   (2.0, 1.75)   (2.0, 2.0)
+
+```
+
+Because some boundary conditions generate an infinite sequences of knots,
+iteration over knots can end up "stuck" iterating along a single axis:
+
+```jldoctest; setup = :(using Interpolations)
+julia> x = [1.0, 1.5, 1.75, 2.0];
+
+julia> etp = LinearInterpolation((x, x), x.*x', extrapolation_bc=(Periodic(), Throw()));
+
+julia> knots(etp) |> k -> Iterators.take(k, 6) |> collect
+6-element Array{Tuple{Float64,Float64},1}:
+ (1.0, 1.0)
+ (1.5, 1.0)
+ (1.75, 1.0)
+ (2.0, 1.0)
+ (2.5, 1.0)
+ (2.75, 1.0)
+
+```
+
+Rearranging the axes so non-repeating knots are first can address this issue:
+
+```jldoctest; setup = :(using Interpolations)
+julia> x = [1.0, 1.5, 1.75, 2.0];
+
+julia> etp = LinearInterpolation((x, x), x.*x', extrapolation_bc=(Throw(), Periodic()));
+
+julia> knots(etp) |> k -> Iterators.take(k, 6) |> collect
+6-element Array{Tuple{Float64,Float64},1}:
+ (1.0, 1.0)
+ (1.5, 1.0)
+ (1.75, 1.0)
+ (2.0, 1.0)
+ (1.0, 1.5)
+ (1.5, 1.5)
+
+```
+
+### Directional Boundary Conditions
+
+If the boundary conditions are directional, the forward boundary condition is
+used to determine if the iterator will generate an infinite sequence of knots.
+
+For example the following extrapolation `etp`, will throw an error for values
+less than `1.0`, but will use `Periodic` extrapolation for values above `2.0`. As a
+result, the iterator will generate an infinite sequence of knots starting at `1.0`.
+
+```jldoctest iterate-directional-unbounded; setup = :(using Interpolations)
+julia> x = [1.0, 1.2, 1.3, 2.0];
+
+julia> etp = LinearInterpolation(x, x.^2, extrapolation_bc=((Throw(), Periodic()),));
+
+julia> kiter = knots(etp);
+
+julia> kiter |> k -> Iterators.take(k, 5) |> collect
+5-element Array{Float64,1}:
+ 1.0
+ 1.2
+ 1.3
+ 2.0
+ 2.2
+
+```
+
+We can also check if the iterator has a length using: `Base.IteratorSize`
+
+```jldoctest iterate-directional-unbounded
+julia> Base.IteratorSize(kiter)
+Base.IsInfinite()
+
+```
+
+Swapping the boundary conditions, results in a finite sequence of knots from
+`1.0` to `2.0`.
+
+```jldoctest iterate-directional-bounded; setup = :(using Interpolations)
+julia> x = [1.0, 1.2, 1.3, 2.0];
+
+julia> etp = LinearInterpolation(x, x.^2, extrapolation_bc=((Periodic(), Throw()),));
+
+julia> kiter = knots(etp);
+
+julia> collect(kiter)
+4-element Array{Float64,1}:
+ 1.0
+ 1.2
+ 1.3
+ 2.0
+
+```
+
+As expected the iterator now has a defined length:
+
+```jldoctest iterate-directional-bounded; setup = :(using Interpolations)
+julia> Base.IteratorSize(kiter)
+Base.HasLength()
+
+julia> length(kiter)
+4
+
+julia> size(kiter)
+(4,)
+
+```
diff --git a/src/Interpolations.jl b/src/Interpolations.jl
index b69fec22..beaff2fe 100644
--- a/src/Interpolations.jl
+++ b/src/Interpolations.jl
@@ -24,7 +24,9 @@ export
     Throw,
 
     LinearInterpolation,
-    CubicSplineInterpolation
+    CubicSplineInterpolation,
+
+    knots
 
     # see the following files for further exports:
     # b-splines/b-splines.jl
@@ -35,8 +37,9 @@ export
 using LinearAlgebra, SparseArrays
 using StaticArrays, WoodburyMatrices, Ratios, AxisAlgorithms, OffsetArrays
 
-using Base: @propagate_inbounds
-import Base: convert, size, axes, promote_rule, ndims, eltype, checkbounds, axes1
+using Base: @propagate_inbounds, HasEltype, HasLength, IsInfinite, SizeUnknown
+import Base: convert, size, axes, promote_rule, ndims, eltype, checkbounds, axes1,
+    iterate, length, IteratorEltype, IteratorSize
 
 abstract type Flag end
 abstract type InterpolationType <: Flag end
@@ -474,5 +477,6 @@ include("io.jl")
 include("convenience-constructors.jl")
 include("deprecations.jl")
 include("lanczos/lanczos.jl")
+include("iterate.jl")
 
 end # module
diff --git a/src/iterate.jl b/src/iterate.jl
new file mode 100644
index 00000000..972cf077
--- /dev/null
+++ b/src/iterate.jl
@@ -0,0 +1,188 @@
+# Similar to ExtrapDimSpec but for only a single dimension
+const ExtrapSpec = Union{BoundaryCondition,Tuple{BoundaryCondition,BoundaryCondition}}
+
+# Type Alias to get Boundary Condition or forward boundary conditions if
+# directional
+const FwdExtrapSpec{FwdBC} = Union{FwdBC, Tuple{BoundaryCondition, FwdBC}}
+
+"""
+    KnotIterator{T,ET}(k::AbstractArray{T}, bc::ET)
+
+Defines an iterator over the knots in `k` based on the boundary conditions `bc`.
+
+# Fields
+- `knots::Vector{T}` The interpolated knots of the axis to iterate over
+- `bc::ET` The Boundary Condition for the axis
+- `nknots::Int` The number of interpolated knots (ie. `length(knots)`)
+
+`ET` is `Union{BoundaryCondition,Tuple{BoundaryCondition,BoundaryCondition}}`
+
+# Iterator Interface
+The following methods defining Julia's iterator interface have been defined
+
+`IteratorSize(::Type{KnotIterator})` -> Will return one of the following
+- `Base.IsInfinite` if the iteration will produces an infinite sequence of knots
+- `Base.HasLength` if iteration will produce a finite sequence of knots
+- `Base.SizeUnknown` if we can't decided from only the type information
+
+`length` and `size` -> Are defined if IteratorSize is HasLength, otherwise will
+raise a MethodError.
+
+`IteratorEltype` will always return `HasEltype`, as we always track the data
+types of the knots
+
+`eltype` will return the data type of the knots
+
+`iterate` Defines iteration over the knots starting from the first one and
+moving in the forward direction along the axis.
+
+# Knots for Multi-dimensional Interpolants
+Iteration over the knots of a multi-dimensional interpolant is done by wrapping
+multiple KnotIterator within `Iterators.product`.
+
+"""
+struct KnotIterator{T,ET <: ExtrapSpec}
+    knots::Vector{T}
+    bc::ET
+    nknots::Int
+    KnotIterator{T,ET}(k::AbstractArray{T}, bc::ET) where {T,ET} = new(k, bc, length(k))
+end
+
+# Dispatch on outputs of k = getknots and bc = etpflag to create one
+# KnotIterator per interpolation axis
+function KnotIterator(k::NTuple{N,AbstractArray}, bc::NTuple{N, ExtrapSpec}) where {N}
+    # One ExtrapSpec per Axis
+    map(KnotIterator, k, bc)
+end
+function KnotIterator(k::NTuple{N,AbstractArray}, bc::BoundaryCondition) where {N}
+    # All Axes use the same BoundaryCondition
+    map(x -> KnotIterator(x, bc), k)
+end
+function KnotIterator(k::AbstractArray{T}, bc::ET) where {T,ET <: ExtrapDimSpec}
+    # Prior dispatches will end up here: Axis knots: k and boundary condition bc
+    KnotIterator{T,ET}(k, bc)
+end
+
+const RepeatKnots = Union{Periodic,Reflect}
+IteratorSize(::Type{KnotIterator}) = SizeUnknown()
+IteratorSize(::Type{KnotIterator{T}}) where {T} = SizeUnknown()
+
+IteratorSize(::Type{KnotIterator{T,ET}}) where {T,ET} = _knot_iter_size(ET)
+_knot_iter_size(::Type{<:BoundaryCondition}) = HasLength()
+_knot_iter_size(::Type{<:RepeatKnots}) = IsInfinite()
+_knot_iter_size(::Type{Tuple{RevBC, FwdBC}}) where {RevBC,FwdBC} = _knot_iter_size(FwdBC)
+
+length(iter::KnotIterator) = _knot_length(iter, IteratorSize(iter))
+_knot_length(iter::KnotIterator, ::HasLength) = iter.nknots
+size(iter::KnotIterator) = (length(iter),)
+
+IteratorEltype(::Type{KnotIterator{T,ET}}) where {T,ET} = HasEltype()
+eltype(::Type{KnotIterator{T,ET}}) where {T,ET} = T
+
+"""
+    knots(itp::AbstractInterpolation)
+    knots(etp::AbstractExtrapolation)
+
+Returns an iterator over knot locations for an AbstractInterpolation or
+AbstractExtrapolation.
+
+Iterator will yield scalar values for interpolations over a single dimension,
+and tuples of coordinates for higher dimension interpolations. Iteration over
+higher dimensions is taken as the product of knots along each dimension.
+
+ie. Iterator.product(knots on first dim, knots on 2nd dim,...)
+
+Extrapolations with Periodic or Reflect boundary conditions, will produce an
+infinite sequence of knots.
+
+# Example
+```jldoctest
+julia> using Interpolations;
+
+julia> etp = LinearInterpolation([1.0, 1.2, 2.3, 3.0], rand(4); extrapolation_bc=Periodic());
+
+julia> Iterators.take(knots(etp), 5) |> collect
+5-element Array{Float64,1}:
+ 1.0
+ 1.2
+ 2.3
+ 3.0
+ 3.2
+
+```
+"""
+function knots(itp::AbstractInterpolation)
+    # Construct separate KnotIterator for each dimension, and combine them using
+    # Iterators.product
+    k = getknots(itp)
+    bc = Throw()
+    iter = KnotIterator(k, bc)
+    length(iter) == 1 ? iter[1] : Iterators.product(iter...)
+end
+
+function knots(etp::AbstractExtrapolation)
+    # Construct separate KnotIterator for each dimension, and combine them using
+    # Iterators.product
+    k = getknots(etp)
+    bc = etpflag(etp)
+    iter = KnotIterator(k, bc)
+    length(iter) == 1 ? iter[1] : Iterators.product(iter...)
+end
+
+# For non-repeating ET's iterate through once
+iterate(iter::KnotIterator) = iterate(iter, 1)
+iterate(iter::KnotIterator, idx::Integer) = idx <= iter.nknots ? (iter.knots[idx], idx+1) : nothing
+
+# For repeating knots state is the knot index + offset value
+function iterate(iter::KnotIterator{T,ET}) where {T,ET <: FwdExtrapSpec{RepeatKnots}}
+    iterate(iter, (1, zero(T)))
+end
+
+# Periodic: Iterate over knots, updating the offset each cycle
+function iterate(iter::KnotIterator{T,ET}, state::Tuple) where {T, ET <: FwdExtrapSpec{Periodic}}
+    state === nothing && return nothing
+    curidx, offset = state[1], state[2]
+
+    # Increment offset after cycling over all the knots
+    if mod(curidx, iter.nknots-1) != 0
+        nextstate = (curidx+1, offset)
+    else
+        knotrange = iter.knots[end] - iter.knots[1]
+        nextstate = (curidx+1, offset+knotrange)
+    end
+
+    # Get the current knot
+    knot = iter.knots[periodic(curidx, 1, iter.nknots)] + offset
+    return knot, nextstate
+end
+
+# Reflect: Iterate over knots, updating the offset after a forward and backwards
+# cycle
+function iterate(iter::KnotIterator{T, ET}, state) where {T, ET <: FwdExtrapSpec{Reflect}}
+    state === nothing && return nothing
+    curidx, offset = state[1], state[2]
+
+    # Increment offset after a forward + backwards pass over the knots
+    cycleidx = mod(curidx, 2*iter.nknots-1)
+    if  cycleidx != 0
+        nextstate = (curidx+1, offset)
+    else
+        knotrange = iter.knots[end] - iter.knots[1]
+        nextstate = (curidx+1, offset+2*knotrange)
+    end
+
+    # Map global index onto knots, and get that knot
+    idx = reflect(curidx, 1, iter.nknots)
+    knot = iter.knots[idx]
+
+    # Add offset to map local knot to global position
+    if 0 < cycleidx <= iter.nknots
+        # Forward pass
+        knot = knot + offset
+    else
+        # backwards pass
+        knot = offset + 2*iter.knots[end] - knot
+    end
+
+    return knot, nextstate
+end
diff --git a/test/iterate.jl b/test/iterate.jl
new file mode 100644
index 00000000..bde0d411
--- /dev/null
+++ b/test/iterate.jl
@@ -0,0 +1,247 @@
+using Test
+using Interpolations
+
+# Unit Tests for Base.IteratorEltype and Base.IteratorSize methods (Type Info Only)
+@testset "iterate - interface" begin
+    import Interpolations.KnotIterator
+    # Always have an eltype as we explicitly track via T
+    @test Base.IteratorEltype(KnotIterator) == Base.HasEltype()
+    @test Base.IteratorEltype(KnotIterator{Int}) == Base.HasEltype()
+
+    # If missing ET type parameter -> SizeUnknown, as could be HasLength or
+    # IsInfinite
+    @test Base.IteratorSize(KnotIterator) == Base.SizeUnknown()
+    @test Base.IteratorSize(KnotIterator{Int}) == Base.SizeUnknown()
+    @test Base.IteratorSize(KnotIterator{Int,Flat}) == Base.HasLength()
+
+    # If ET is Directional -> Size based on Fwd Direction
+    @test Base.IteratorSize(KnotIterator{Int,Tuple{Flat,Periodic}}) == Base.IsInfinite()
+    @test Base.IteratorSize(KnotIterator{Int,Tuple{Periodic,Flat}}) == Base.HasLength()
+end
+
+# eltype units tests of KnotIterator
+@testset "iterate - eltype" for T ∈ [ Int, Float64, Any ]
+    x = convert(Vector{T}, collect(1:5))
+    itp = LinearInterpolation(x, x.^2)
+    kiter = knots(itp)
+    @test Base.IteratorEltype(typeof(kiter)) == Base.HasEltype()
+    @test typeof(kiter) <: Interpolations.KnotIterator{T}
+    @test eltype(kiter) == T
+end
+
+# size/length units tests for KnotIterator
+@testset "iterate - size - directional"  begin
+    ExtrapBC = [Throw(), Flat(), Periodic(), Reflect()]
+    @testset "RevBC: $RevBC, FwdBC: $FwdBC" for FwdBC ∈ ExtrapBC, RevBC ∈ ExtrapBC
+        etp = LinearInterpolation(1:10, rand(10),
+            extrapolation_bc=((RevBC, FwdBC),)
+        )
+        kiter = knots(etp)
+        @test typeof(kiter) <: Interpolations.KnotIterator
+        ktype = KnotIterator{Int, Tuple{typeof(RevBC), typeof(FwdBC)}}
+        @test typeof(kiter) == ktype
+
+        if typeof(FwdBC) <: Union{Periodic, Reflect}
+            @test Base.IteratorSize(kiter) == Base.IsInfinite()
+            @test_throws MethodError length(kiter)
+            @test_throws MethodError size(kiter)
+        else
+            @test Base.IteratorSize(kiter) == Base.HasLength()
+            @test length(kiter) == 10
+            @test size(kiter) == (10,)
+        end
+    end
+end
+
+# Unit Tests for 1D Interpolations
+@testset "iterate - 1d interpolation" begin
+    itp = interpolate(1:5, BSpline(Linear()))
+    kiter = knots(itp)
+
+    # Single axis -> KnotIterator
+    @test typeof(kiter) <: Interpolations.KnotIterator
+
+    # Check size, length and eltype
+    @test Base.IteratorSize(kiter) === Base.HasLength()
+    @test length(kiter) == 5
+    @test size(kiter) == (5,)
+    @test Base.IteratorEltype(kiter) == Base.HasEltype()
+    @test eltype(kiter) <: Int
+
+    # Check contents
+    @test collect(kiter) == [1, 2, 3, 4, 5]
+
+    # Sample for loop
+    kout = []
+    for (k,) ∈ kiter
+        push!(kout, k)
+    end
+    @test kout == collect(kiter)
+end
+
+# Unit Test for 2D Interpolations
+@testset "iterate - 2d interpolation" begin
+    itp = interpolate(rand(3,3), BSpline(Constant()))
+    kiter = knots(itp)
+
+    # 2D iterator -> ProductIterator wrapping KnotIterators
+    @test typeof(kiter) <: Iterators.ProductIterator
+    @test eltype(kiter.iterators) <: Interpolations.KnotIterator
+
+    # Check size, length and eltype
+    @test Base.IteratorSize(kiter) == Base.HasShape{2}()
+    @test length(kiter) == 9
+    @test size(kiter) == (3, 3)
+    @test Base.IteratorEltype(kiter) == Base.HasEltype()
+    @test eltype(kiter) <: Tuple{Int, Int}
+
+    # Collect knots and check again
+    k = collect(kiter)
+    @test length(k) == 9
+    @test size(k) == (3,3)
+    @test eltype(k) <: Tuple{Int, Int}
+    @test k == Iterators.product(1:3, 1:3) |> collect
+
+    # Sample for loop
+    kout = []
+    for (kx, ky) ∈ kiter
+        push!(kout, (kx, ky))
+    end
+    @test kout == Iterators.product(1:3, 1:3) |> collect |> vec
+
+end
+
+# Unit Tests for 1D extrapolations
+@testset "1D - iterate - $itp" for itp ∈ [ BSpline(Constant()) ]
+    x = 1.0:5.0
+    itp = interpolate(x.^2, itp)
+    @test typeof(itp) <: AbstractInterpolation
+    kiter = knots(itp)
+
+    # Checkout construction before proceeding
+    @test typeof(kiter) <: Interpolations.KnotIterator
+    @test eltype(kiter) <: Int
+    @test length(kiter) == 5
+    @test size(kiter) == (5,)
+    @test collect(kiter) == x
+
+    # Non-repeating Knots -> Should iterate over once and be done
+    @testset "extrapolation - $etp" for etp ∈ [ Throw(), Flat(), Line()]
+        extrp = extrapolate(itp, etp)
+        kiter = knots(extrp)
+        @test typeof(kiter) <: Interpolations.KnotIterator
+        @test eltype(kiter) <: Int
+        @test length(kiter) == 5
+        @test size(kiter) == (5,)
+        @test collect(kiter) == x
+    end
+
+    # Repeating Knots -> Should Iterate indefinitely
+    @testset "extrapolation - Periodic" begin
+        extrp = extrapolate(itp, Periodic())
+        k = knots(extrp)
+        @test typeof(k) <: Interpolations.KnotIterator
+        @test_throws MethodError length(k)
+        @test_throws MethodError size(k)
+        k10 = Iterators.take(k, 10) |> collect
+        @test k10 ≈ collect(1:10)
+        @test extrp.(k10) ≈ vcat(x[1:end-1].^2, x[1:end-1].^2, x[1:2].^2)
+    end
+    @testset "extrapolation - Reflect" begin
+        extrp = extrapolate(itp, Reflect())
+        k = knots(extrp)
+        @test typeof(k) <: Interpolations.KnotIterator
+
+        # Check length and size throw Method Errors (As Undefined)
+        @test Base.IteratorSize(k) == Base.IsInfinite()
+        @test_throws MethodError length(k)
+        @test_throws MethodError size(k)
+
+        k10 = Iterators.take(k, 10) |> collect
+        @test k10 ≈ collect(1:10)
+        @test extrp.(k10) ≈ vcat(x.^2, reverse(x[1:end - 1].^2), x[2].^2)
+    end
+end
+
+# Unit tests for iteration on an uneven grid - Periodic
+# Knots should repeat indefinitely with the first and last knots being co-located
+@testset "iterate - uneven - Periodic" begin
+    x = [1.0, 1.3, 2.4, 3.2, 4.0]
+    etp = LinearInterpolation(x, x.^2, extrapolation_bc=Periodic())
+    kiter = knots(etp)
+
+    @test typeof(kiter) <: Interpolations.KnotIterator
+    @test Base.IteratorSize(kiter) == Base.IsInfinite()
+    @test_throws MethodError length(kiter)
+    @test_throws MethodError size(kiter)
+
+    k = Iterators.take(kiter, 10) |> collect
+    @test k == [1.0, 1.3, 2.4, 3.2, 4.0, 4.3, 5.4, 6.2, 7.0, 7.3]
+    @test etp.(k) ≈ vcat(x[1:end-1], x[1:end-1], x[1:2]).^2
+end
+
+# Unit tests for iteration on an uneven grid - Reflect
+@testset "iterate - uneven - Reflect" begin
+    x = [1.0, 1.3, 2.4, 3.2, 4.0]
+    etp = LinearInterpolation(x, x.^2, extrapolation_bc=Reflect())
+    kiter = knots(etp)
+
+    @test typeof(kiter) <: Interpolations.KnotIterator
+    @test Base.IteratorSize(kiter) == Base.IsInfinite()
+    @test_throws MethodError length(kiter)
+    @test_throws MethodError size(kiter)
+
+    k = Iterators.take(kiter, 10) |> collect
+    @test k == [1.0, 1.3, 2.4, 3.2, 4.0, 4.8, 5.6, 6.7, 7.0, 7.3]
+    @test etp.(k) ≈ vcat(x, reverse(x[1:end - 1]), x[2]).^2
+end
+
+# Unit tests for 2D iteration with directional boundary conditions that are
+# bounded (ie. knots do not repeat indefinitely)
+@testset "2D - iteration - bounded - $bc" for bc ∈ [Line(), (Throw(), Line())]
+    etp = LinearInterpolation(([1, 2, 3], [1, 2, 3]), rand(3, 3);
+        extrapolation_bc=(Throw(), bc)
+    )
+    kiter = knots(etp)
+    @test typeof(kiter) <: Iterators.ProductIterator
+    @test eltype(kiter.iterators) <: Interpolations.KnotIterator
+
+    @test Base.IteratorSize(kiter) == Base.HasShape{2}()
+    @test length(kiter) == 9
+    @test size(kiter) == (3, 3)
+
+    @test Base.IteratorEltype(kiter) == Base.HasEltype()
+    @test eltype(kiter) == Tuple{Int,Int}
+
+    k = Iterators.take(kiter, 5) |> collect
+    @test length(k) == 5
+    @test typeof(k) <: Vector{Tuple{Int,Int}}
+    kexp = Iterators.product(1:3, 1:3) |> x -> Iterators.take(x, 5) |> collect
+    @test k == kexp
+end
+
+# Unit tests for 2D iteration with directional boundary conditions that are
+# unbounded (ie. knots do repeat indefinitely)
+@testset "2D - iteration - Unbounded - $bc" for bc ∈ [Periodic(), (Throw(), Periodic())]
+    etp = LinearInterpolation(([1, 2, 3], [1, 2, 3]), rand(3, 3);
+        extrapolation_bc=(Line(), bc)
+    )
+    kiter = knots(etp)
+    @test typeof(kiter) <: Iterators.ProductIterator
+    @test eltype(kiter.iterators) <: Interpolations.KnotIterator
+
+    # Check length and size throw ArgumentErrors (via Iterators.product)
+    @test Base.IteratorSize(kiter) == Base.IsInfinite()
+    @test_throws ArgumentError length(kiter)
+    @test_throws ArgumentError size(kiter)
+
+    @test Base.IteratorEltype(kiter) == Base.HasEltype()
+    @test eltype(kiter) == Tuple{Int,Int}
+
+    k = Iterators.take(kiter, 20) |> collect
+    @test length(k) == 20
+    @test typeof(k) <: Vector{Tuple{Int,Int}}
+    kexp = Iterators.product(1:3, Iterators.countfrom(1)) |>
+        x -> Iterators.take(x, 20) |> collect
+    @test k == kexp
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index e34ff1fb..dbc3d8df 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -49,5 +49,6 @@ const isci = get(ENV, "CI", "") in ("true", "True")
     include("io.jl")
     include("convenience-constructors.jl")
     include("readme-examples.jl")
+    include("iterate.jl")
 
 end