Gradients for extrapolation objects are wrong when one of the variables is extrapolated

`Interpolations.gradient` returns the wrong gradients (when comparing against finite difference) when extrapolating. This seems like a bug to me, but please explain if this behavior is expected.

See example below.

```
using Interpolations
using Test

@testset "all" begin

@testset "cubic spline" begin

f1 = (x,y) -> sin(x) * 3*cos(y)

x = range(-π/2, stop=π/2, length=11)
y = range(-π/2, stop=π/2, length=11)

z = [f1.(xx, yy) for yy in y, xx in x]

zfit = Interpolations.interpolate(z, Interpolations.BSpline(Interpolations.Cubic(Interpolations.Line(Interpolations.OnGrid()))))
zfit_scaled = Interpolations.scale(zfit, x, y)
zfits_scaled_extrap = Interpolations.extrapolate(zfit_scaled, Interpolations.Line())

ϵ = 1e-6

xcheck = π * 1.0
ycheck = π * 1.0

∂f∂x_fd = (zfits_scaled_extrap(xcheck + ϵ, ycheck) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_scaled_extrap(xcheck, ycheck + ϵ) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_scaled_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4 atol=1e-12
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4 atol=1e-12

xcheck = π * 0.4
ycheck = π * 0.4

∂f∂x_fd = (zfits_scaled_extrap(xcheck + ϵ, ycheck) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_scaled_extrap(xcheck, ycheck + ϵ) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_scaled_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4 atol=1e-12
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4 atol=1e-12

# extrapolate x axis
xcheck = π * 1.0
ycheck = π * 0.4

∂f∂x_fd = (zfits_scaled_extrap(xcheck + ϵ, ycheck) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_scaled_extrap(xcheck, ycheck + ϵ) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_scaled_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4 atol=1e-10
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4 atol=1e-10 # fails

# extrapolate y axis
xcheck = π * 0.4
ycheck = π * 1.0

∂f∂x_fd = (zfits_scaled_extrap(xcheck + ϵ, ycheck) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_scaled_extrap(xcheck, ycheck + ϵ) - zfits_scaled_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_scaled_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4 atol=1e-12 # fails
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4 atol=1e-12

end

@testset "linear grid" begin

f1 = (x,y) -> sin(x) * 3*cos(y)

x = collect(range(-π/2, stop=π/2, length=11))
y = collect(range(-π/2, stop=π/2, length=11))

# make them not equidistant
x[3] += 0.1
y[5] += 0.1

z = [f1.(xx, yy) for yy in y, xx in x]

knots = (x, y)

zfit = Interpolations.interpolate(knots, z, Interpolations.Gridded(Interpolations.Linear()))
zfits_extrap = Interpolations.extrapolate(zfit, Interpolations.Line())

ϵ = 1e-4

xcheck = π * 1.0
ycheck = π * 1.0

∂f∂x_fd = (zfits_extrap(xcheck + ϵ, ycheck) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_extrap(xcheck, ycheck + ϵ) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4 atol=1e-12
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4 atol=1e-12

xcheck = π * 0.45
ycheck = π * 0.45

∂f∂x_fd = (zfits_extrap(xcheck + ϵ, ycheck) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_extrap(xcheck, ycheck + ϵ) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4

xcheck = π * 1.0
ycheck = π * 0.45

∂f∂x_fd = (zfits_extrap(xcheck + ϵ, ycheck) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_extrap(xcheck, ycheck + ϵ) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4 # fails

# extrapolate
xcheck = π * 0.45
ycheck = π * 1.0

∂f∂x_fd = (zfits_extrap(xcheck + ϵ, ycheck) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂y_fd = (zfits_extrap(xcheck, ycheck + ϵ) - zfits_extrap(xcheck, ycheck)) / ϵ
∂f∂x, ∂f∂y = Interpolations.gradient(zfits_extrap, xcheck, ycheck)
@test ∂f∂x_fd ≈ ∂f∂x rtol=1e-4 # fails
@test ∂f∂y_fd ≈ ∂f∂y rtol=1e-4

end

end
```

This returns

```
cubic spline: Test Failed at /app/gradient_interpolations.jl:47
Expression: ≈(∂f∂y_fd, ∂f∂y, rtol = 0.0001, atol = 1.0e-10)
Evaluated: -1.3405915728625928 ≈ 1.82219278625514e-16 (rtol=0.0001, atol=1.0e-10)
cubic spline: Test Failed at /app/gradient_interpolations.jl:56
Expression: ≈(∂f∂x_fd, ∂f∂x, rtol = 0.0001, atol = 1.0e-12)
Evaluated: -3.2627963864051424 ≈ -2.853013321706128 (rtol=0.0001, atol=1.0e-12)
linear grid: Test Failed at /app/gradient_interpolations.jl:106
Expression: ≈(∂f∂y_fd, ∂f∂y, rtol = 0.0001)
Evaluated: -0.7221353894237836 ≈ 2.861851643454053e-17 (rtol=0.0001)
linear grid: Test Failed at /app/gradient_interpolations.jl:115
Expression: ≈(∂f∂x_fd, ∂f∂x, rtol = 0.0001)
Evaluated: -3.6730303186716107 ≈ -2.950894929250398 (rtol=0.0001)
```

Showing that these gradients are substantially off when compared to the finite difference result.

A similar thing is observed when plotting the result.

```
f1 = (x,y) -> sin(x) * 3*cos(y)
x = range(-π/2, stop=π/2, length=11)
y = range(-π/2, stop=π/2, length=11)
z = [f1.(xx, yy) for yy in y, xx in x]

zfit = Interpolations.interpolate(z, Interpolations.BSpline(Interpolations.Cubic(Interpolations.Line(Interpolations.OnGrid()))))
zfit_scaled = Interpolations.scale(zfit, x, y)
zfits_scaled_extrap = Interpolations.extrapolate(zfit_scaled, Interpolations.Line())

yplot = range(-π, stop=π, length=101)
plot(yplot, zfits_scaled_extrap.(π, yplot), label="z")
plot(yplot, [Interpolations.gradient(zfits_scaled_extrap, π, yy)[2] for yy in yplot], label="∂z∂y")
```

which returns