I have a geographical matrix M of dimensions nLat=181 and nLon=221 , with lats = 30:0.25:90 and lons=-30:0.25:25.
I want to compute the gradient vector of M w.r.t to lats , lons.
I have tried ForwardDiff , Iterations , ScatteredIterations , Flux and my problem is that they all accept continuous functions as input.
I have tried the method interpolating , or this one, but to no avail.
Here is my current code
x = lats
y = lons
points = (x,y)
samples = M
itp = Interpolations.interpolate(points, samples, Gridded(Linear()))
itpGradient = Interpolations.gradient.(itp,x,y)
I get the following error
DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 181 and 221")
You could adapt the 1D solution in the post linked to 2D as follows:
using Interpolations
f(x,y) = x^2 + y^2 # function used just to generate data
x = y = -2:0.02:2
z = f.(x, y') # z is the numeric data at coordinates (x, y)
itp = interpolate((x,y), z, Gridded(Linear()));
grad = gradient.(Ref(itp), x, y')
If the data matrix z is noisy, we could fit smoothing 2D splines before estimating the gradients. An example using Dierckx:
using Interpolations, Dierckx
f(x,y) = x^2 + y^2
x, y = -2:0.02:2, -3:0.02:3
nx, ny = length(x), length(y)
z = f.(x, y') + rand(nx, ny) # noisy data matrix over grid defined by (x,y)
s0 = 0.01*nx*ny*hypot(extrema(z)...) # test different smoothing factors
spl = Spline2D(x, y, z; kx=3, ky=3, s=s0)
zsmooth = evalgrid(spl, x, y)
itp = interpolate((x,y), zsmooth, Gridded(Linear()));
grad = gradient.(Ref(itp), x, y')
