Normally we do 1-degree of latitude constant (y-dimension) and equal to 110000 m. So your dy
depends on how many degrees per grid point you have (yresolution). However, dx
will vary according to the latitude.
dy=110000*yesolution;
Now, we may use centered finite differences to compute what you want. (Below is a Matlab code, but I believe generic enough to be reproducible in other languages).
for y=2:length(lat)-1 dx=abs(110000*cos(latx(y)*(2*pi/360))*xresolution); for x=2:length(lon)-1 div(x,y,1)= (u(x+1,y,1)-u(x-1,y,1))/(2*dx) + (v(x,y+1,1)-v(x,y-1,1)/(2*dy); end end
Note that you will have an empty frame around your divergence field when x=1, x=max(x), y=1 and y=max(y), once they do not have two neighboors to compute the differences. The same is observed when you do this with hdivg
function in Grads, for example.
By the way, this is based on and yields the same result as Grads (I have checked!) using cdiff to reproduce hdivg
.
Also note that the cosine function here works with radians, so if you compute cosines in degrees directly (i.e. cosd
function in Matlab), you should omit the term scaling by (2*pi/360)
, which is just a conversion.
Hope it helps!
You say you have a grid, and your point ($x$,$y$) is at the ith and jth gridpoint, where i indicates the left-right index on the grid and j indicates the north-south index on the grid. Therefore $u$ and $v$ can be expressed as $u(x_{i,j},y_{i,j})$ and $v(x_{i,j},y_{i,j})$.
$$\frac{du}{dx}\approx\frac{\Delta u}{\Delta x}=\frac{u(x+\Delta x,y)-u(x-\Delta x,y)}{\Delta x}=\frac{u(x_{i+1,j},y)-u(x_{i-1,j},y)}{x_{i+1,j}-x_{i-1,j}}=\frac{u_{i+1,j}-u_{i-1,j}}{x_{i+1,j}-x_{i-1,j}}$$
A similar process can be applied to the v component, just change the index, and the variables