Project

General

Profile

expr - New options with release 2.0.5

Added by Ralf Mueller over 2 years ago

Fellow expr-essionists!

In release 2.0.5 CDO there was an addition to the expr family of operators regarding the mathematical functions with input values:

I won't cover the first 4 functions because their purpose is pretty obvious. But hypot and atan2 are more interesting

  1. hypot(x,y) is a shortcut for √(x 2 + y 2). The alternative operator combination is
    cdo -sqrt -add -mul -selvar,x <ifile> -selvar,x <ifile> -mul -selvar,y <ifile> -selvar,y <ifile> <ofile>
    which is quite lengthy compared to hypot(x,y). So it's good to keep in mind, that expr offers an alternative here - an even shorter one than the former solution sqrt(u*u+v*v).
  2. atan2(x,y) calculates the inverse tangent of x/y but with an improved error handling (like what happens in case x=0 or y=0). Please check the linked documentation above for more information on that.
Now here is the task: Let's assume you have northward and eastward ocean (or wind) velocity vectors u and v and you would like to know their direction in a more convenient coordinate system like for navigation with a compass on a ship, an airplane or a simple map. In this system the main directions are
  • north: 0 degrees
  • east: 90 degrees
  • south: 180 degrees
  • west: 270 degrees

In fact it is a rotated polar coordinated system. Unfortunately this is not what atan2 returns natively (it instead follows the mathematical definitions of polar coordinates): It (a) returns angles in radian and (b) it returns +pi/2 for a vector, which points to north. So we need to tweak the results a bit to make atan2 usable.

  • transformation of radian to degrees can be done with the deg function in expr: deg(atan2(u,v))
  • shifting the results by 90 degrees is easy: We just subtract the result from 90: angle=90.0-deg(atan2(u,v))

A first CDO attempt looks like

cdo expr,"arg=90.0-deg(atan2(v,u));vel=hypot(u,v)" veloc.nc velarg.nc

The only issue left is that there are no negative angles in the navigation coordinate system. But this can be ensured after the atan2 computation with a ternary operator: arg=(arg<0)?arg+360. Which reads like: if the angle arg is smaller than zero, add 360 and leave it untouched otherwise. The final CDO call now is

cdo expr,"arg=90.0-deg(atan2(v,u));arg=(arg<0)?arg+360:arg;vel=hypot(u,v)" veloc.nc velarg.nc
In case you want to add the new variables arg and vel to the existing variables u and v in veloc.nc to keep them attached to one another, you can use the aexpr operator instead of expr.

AND NOW for the WIND - yes, I know: In meteorology wind directions are handled differently. Instead of naming directions wrt. where the air flows, the naming follows the direction from where the air comes. So vectors are the same, but naming is different: west wind causes an eastward ocean current ... don't ask

Another simplification is using the equality:

90.0-deg(atan2(v,u)) = deg(atan2(u,v))
which turns the above call into
cdo expr,"arg=deg(atan2(u,v));vel=hypot(u,v)" veloc.nc velarg.nc

Happy navigating with CDO!