Hi,
I am trying to calculate the divergence of two vector fields (NO2 flux and wind fields (U & V vectors). Firstly, the NO2 flux is multiplied by the u and v wind vectors and generates the zonal and meridional fluxes. Secondly, the calculation of divergence must perform using zonal and meridional fluxes. As per knowledge, divergence is the summation of x and y derivatives.
I used grdmath DDX and DDY to calculate the x and y first derivatives individually for zonal and meridional fluxes. Here I want to use the second-order central finite difference method to calculate divergence. Do I need to reapply DDX and DDY again on the first derivatives to calculate the second derivatives? Or Do I directly use D2DX2 and D2DY2 for the second derivative?
Additionally, can I use grdgradient for the divergence? I know gradient is a vector field, whereas divergence is a scalar field, but I am asking for the sake of curiosity.
Good point. I just did a PR to clarify that all derivatives in grdmath are centred, So if you need D2DX2 you should use that directly and not repeat DDX. If this is a geographic grid (which it sounds like) then the DX shrinks with latitude.
And no on the grdgradient question. Perhaps we should add DIVERG or something to grdmath.
Thank you, @pwessel, for providing the clarification. This helps to resolve my doubts. It would be fascinating if grdmath could directly perform divergence with DIVERG command. The cherry on the cake!
Alternatively, you just create a grdmath macro called DIVERGENCE and set it to mean the sum of those two 2nd derivatives; see man page.
1 Like
@pwessel, Is it possible to use DDX and DDY on the XY grid instead of the geographical grid?
Of course, becomes Cartesian
@pwessel, Thank you for the clarification. I was asking because my XY grid is in a meter, which is basically distance. I doubt that DDX and DDY could directly consider this for derivatives.
Additionally, I know gmt mapproject can be used for the map transformation. I tried to convert the geographical grid into an XY grid in meters but couldn’t achieve the outcome.
gmt mapproject latlon.xyz -R-105/-101/30.8/33.2 -Jm0.025m > xy.xyz
My grid is in 0.025-degree x 0.025-degree spatial resolution.
Mmm, I think that maybe you should use grdproject for this conversion.
https://docs.generic-mapping-tools.org/dev/grdproject.html