Compute horizontal and vertical derivatives

[arctica]

Hello,

Given an array zdata (1200x1200), how can the vertical derivative be calculated given that dFdy and dFdx can be computed from:

// Compute dF/dx
dFdx = diff(zdata, 1, 2); // Difference along columns
dFdx = [dFdx, dFdx(:, $)]; // Maintain the same grid size by appending the last column

// Compute dF/dy
dFdy = diff(zdata, 1, 1); // Difference along rows
dFdy = [dFdy; dFdy($, :)]; // Maintain the same grid size by appending the last row

What I need to find is dF/dz. Using just diff(zdata) gives a row of data not a grid.

Lester

[QwillDune]

Hello,

If I understand correctly, you want to calculate gradients.

See this post

Calculating gradients

Hi I am wondering if Scilab has a built-in function similar to Python numpy.gradient, which takes an ‘edge-order’ … I can manually create a first-order gradient function in Scilab like this: function [grad]=dydx(y, x) k = length(x); grad = zeros(k); // initialize grad(1) = (y(2) - y(1))./(x(2) - x(1)); grad(k) = (y(k) - y(k-1))./(x(k) - x(k-1)); if k>2 then // for i = 2:k-1 do // grad(i) = ( (y(i+1) - y(i))./(x(i+1) - x(i)) + .. // …

BR
François

[cfuttrup]

I agree, the function I proposed in the article references by François should do the trick. I can propose an improvement to the function, checking that the input is validated a bit more.

[arctica]

Taking a slightly different approach, one can work the problem of computing dF/dx, dF/dy and dF/dz using FFT. This necessitates the need to pad the input data grid on all boundaries to avoid edge effects as a function of processing the data in the FFT domain. The final stage would be clipping the grid back to the original data dimensions.

The following works to generate the first derivatives in X, Y and Z dimensions, with X and Y derived from an azimuthal derivative function. This replicates the basic methodology of grdff in the Generic Mapping Tools (GMT) and thus provides a method of verifying that the results are pretty much the same with only tiny differences.

GMT syntax:
grdfft faa_test_1m.grd -Gfaa3_dx.grd -A90 -fg -V -N2400/2400 : dF/dx
grdfft faa_test_1m.grd -Gfaa3_dy.grd -A0 -fg -V -N2400/2400 : dF/dy
grdfft faa_test_1m.grd -Gfaa3_dz.grd -D -fg -V -N2400/2400 : dF/dz

Replicated in Scilab 2024.1.0

function [array_pad]=padding(array,ncols,nrows,pad)
    tmp = [array ; array($, :) .*. ones(pad, 1)];
    tmp = [tmp  tmp(:, $) .*. ones(1,pad)];
    tmp = [repmat(tmp(1, : ), pad, 1) ; tmp];
    array_pad = [repmat(tmp(:,1) , 1, pad) tmp]
endfunction

function [array_clip]=clip(array_pad,pad)
    array_clip = array_pad(pad+1:$-pad, pad+1:$-pad);
endfunction

function dxy_grid = directional_derivative_fft(grid, dx, dy, n, theta)
    // grid: Input 2D array representing the grid data
    // dx/dy: Spacing in cartesian distance (1 arc minute ~ 1850m)
    // theta: Angle in radians for the directional derivative

    // Get the dimensions of the grid
    [nx, ny] = size(grid);
    
    // Pad original grid to avoid FFT edge effects
    // new grid = 2*nx 2*ny should be fine to eliminate boundary effects
    pad = nx/2; // padding distance at each edge (user-defined)
    grid_padded = padding(grid,n*nx,n*ny,pad);

    // Perform 2D FFT on the grid
    fft_grid = fft2(grid_padded);
    [NX, NY] = size(grid_padded); // process padded grid

    // Create vectors of wave numbers for x and y directions
    kx = (2 * %pi) * [0:floor(NX/2), -ceil(NX/2)+1:-1] / (NX * dx);
    ky = (2 * %pi) * [0:floor(NY/2), -ceil(NY/2)+1:-1] / (NY * dy);

    // Create meshgrid of wave numbers
    [KY, KX] = meshgrid(ky, kx);

    // Calculate the Fourier-space derivatives
    fft_dx = %i * KX .* fft_grid;
    fft_dy = %i * KY .* fft_grid;

    // Compute the directional derivative using the angle theta
    fft_derivative = cos(theta) .* fft_dx + sin(theta) .* fft_dy;

    // Perform the inverse FFT to get the derivative in the spatial domain
    dxy_grid_complex = ifft(fft_derivative);

    // Extract the real part of the result
    dxy_grid = real(dxy_grid_complex);
    
    // Clip resulting grid to original limits
    dxy_grid = clip(dxy_grid,pad)

endfunction

function dz_grid = vertical_derivative_fft(grid, dx, dy, n)

    [nx, ny] = size(grid);
    pad = nx/2; // padding distance at each edge
    grid_padded = padding(grid,n*nx,n*ny,pad);

    fft_grid = fft2(grid_padded);
    [NX, NY] = size(grid_padded);

    kx = (2 * %pi) * [0:floor(NX/2), -ceil(NX/2)+1:-1] / (NX * dx);
    ky = (2 * %pi) * [0:floor(NY/2), -ceil(NY/2)+1:-1] / (NY * dy);
    [KY, KX] = meshgrid(ky, kx);
    K = sqrt(KX.^2 + KY.^2);
    K(1,1) = 1e-10; // Assign a small non-zero value

    fft_derivative = fft_grid .* K;
    dz_grid_complex = ifft(fft_derivative);
    dz_grid = real(dz_grid_complex);
    dz_grid = clip(dz_grid,pad);

endfunction

// Load a netCDF file
NC=h5open('faa_test_1m.grd','r');
latitude=h5read(NC,'/lat');
ny = length(latitude)
longitude=h5read(NC,'/lon');
nx = length(longitude) 
F_faa=h5read(NC,'/z');
h5close(NC)

// test
dy_grid = directional_derivative_fft(F_faa, 1850, 1850, 2, -%pi/2); // dF/dy
dx_grid = directional_derivative_fft(F_faa, 1850, 1850, 2, 0); // dF/dx
HGM = sqrt(dx_grid.^2 + dy_grid.^2); // Horizontal; Gradient Magnitude
dz_grid = vertical_derivative_fft(F_faa, 1850, 1850, 2); // dF/dz

tilt = atan(dz_grid./HGM); // Tilt derivative

// Clip data range to +/- 45 degrees (+/- %pi/4)
// Data range used to compute tilt depth
tilt_clip = tilt;
tilt_clip(tilt < -%pi/4 | tilt > %pi/4) = %nan;

Hope this is useful or perhaps can be improved upon. Note: the input netCDF data are in Geographic coordinates so a dx/dy scale of 1850m is applied to correspond to the linear distance equivalent to ~1 arc minute; if the data are already in linear coordinates (eg UTM etc) then dx and dy scales would be 1.

Lester