Commit f6d55b42 authored by Patrick Wagner's avatar Patrick Wagner

Resolve "Add rotated velocities"

parent e8277fbb
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -2,9 +2,153 @@ from itertools import chain, islice
import numpy as np
from functools import reduce
import xarray as xr
import xarray.ufuncs as uf
def section_leg_indices(ji0=None, ji1=None):
def heading(ll_pairs):
"""
Compute the heading of a section leg defined by a two latlon pairs
Parameter
---------
ll_pairs : list of tuples
List of tuples of length 2 that defines section leg.
`[(j0,i0),(j1,i1)]`
Returns
--------
heading: Float
"""
conv = np.pi/180 # for degree to radian conversion
xa = ll_pairs[0][1]*conv
xb = ll_pairs[1][1]*conv
ya = -uf.log(np.tan(np.pi/4.-conv*ll_pairs[0][0]/2))
yb = -uf.log(np.tan(np.pi/4.-conv*ll_pairs[1][0]/2))
xb_xa = (xb-xa) % (2*np.pi)
if (xb_xa >= np.pi):
xb_xa -= 2 * np.pi
# if (xb_xa <= -np.pi): # Not needed because xb_xa >= 0
# xb_xa += 2 * np.pi
heading = uf.arctan2(xb_xa, (yb-ya))
if (heading < 0):
heading += 2*np.pi
return heading
def coordinate_rotation(lat, lon, dim):
"""
Compute the local inclination of the coordinate system, i.e. the heading
of every u or v interface. Which interface is processed is controlled by
paramter dim. Note that for the y case the result is offset by -pi/2.
Parameter
---------
lat : xarray.DataArray
Array of shape (1,y,x) giving latitude values
lon : xarray.DataArray
Array of shape (1,y,x) giving longitude values
dim: string
String to specify which dimension ('x' or 'y')
is processed, i.e. inclination of v or u interface.
Returns
--------
heading: xarray.DataArray
Array of shape (1,y,x) giving local inclination in radian.
"""
conv = np.pi/180 # for degree to radian conversion
xa = lon*conv
ya = -uf.log(uf.tan(np.pi/4.-conv*lat/2))
if dim == 'y':
xb = lon.roll(y=-1, roll_coords=True)*conv
yb = -uf.log(uf.tan(np.pi/4.-conv*lat.roll(y=-1, roll_coords=True)/2))
if dim == 'x':
xb = lon.roll(x=-1, roll_coords=True)*conv
yb = -uf.log(uf.tan(np.pi/4.-conv*lat.roll(x=-1, roll_coords=True)/2))
xb_xa = (xb-xa) % (2*np.pi)
xb_xa = (xb_xa - 2*np.pi).where(xb_xa >= np.pi, xb_xa)
xb_xa = (xb_xa + 2*np.pi).where(xb_xa <= -np.pi, xb_xa)
heading = uf.arctan2(xb_xa, (yb-ya))
heading = uf.arctan2(xb_xa, (yb-ya))
if dim == 'y':
heading[0, -1, :] = heading[0, -2, :]
if dim == 'x':
heading -= np.pi/2
heading[0, :, -1] = heading[0, :, -2]
return heading # *(-1)
def rotate_velocities(u, v, theta):
"""
Rotate velocities by angel theta in clockwise direction.
u_rot = u * xr.ufuncs.cos(theta) - v * np.sin(theta)
v_rot = u * xr.ufuncs.sin(theta) + v * np.cos(theta)
Parameter
---------
v: xarray.Dataset
Array of shape N providing meridional velocity component
u: xarray.Dataset
Array of shape N providing zonal velocity component
theta: xarray.Dataset or scalar
Angle [radian] used to rotate velocities. At least one dimension
as to match the dimension of u and v or theta has to be a scalar
value
Returns
-------
v_rot: xarray.Dataset
Array of shape N providing rotated meridional velocity component
u_rot: xarray.Dataset
Array of shape N providing rotated zonal velocity component
-------
"""
u_rot = u * xr.ufuncs.cos(theta) - v * np.sin(theta)
v_rot = u * xr.ufuncs.sin(theta) + v * np.cos(theta)
return u_rot, v_rot
def section_indices(ji_pairs=None, mesh_hgr=None):
"""Forms a staircase of consectuive segments from
ij_pairs[0] to ij_pairs[-1]. Number of segments will be one
less than number of input points.
Parameters
----------
ji_pairs : list of tuples
List of tuples of lenth n that defines section nodes.
`[(j0,i0),(j1,i1), ..., (jn,in)]`
Return
------
ji_pairs : Iterator
Iterator of tuples of length n
"""
if mesh_hgr:
ji_pairs_for_sub_sections = map(section_leg_indices,
ji_pairs[:-1], ji_pairs[1:],
[mesh_hgr]*(len(ji_pairs)-1))
else:
ji_pairs_for_sub_sections = map(section_leg_indices,
ji_pairs[:-1], ji_pairs[1:])
ji_pairs_stitched_together = reduce(
lambda x, y: chain(x, islice(y, 1, None)),
ji_pairs_for_sub_sections)
return ji_pairs_stitched_together
def section_leg_indices(ji0=None, ji1=None, mesh_hgr=None):
"""Return sequence of (j,i) forming a staircase of from ji0 to ji1.
Parameters
......@@ -46,40 +190,22 @@ def section_leg_indices(ji0=None, ji1=None):
ii_leg = np.round(np.linspace(i0, i1, n_of_points))
jj_leg = np.append([j0], np.cumsum(np.diff(ii_leg) == 0)
* np.sign(j1-j0)+j0)
ii_leg = np.append(i0-dx, ii_leg)
i0 -= dx
ii_leg = np.append(i0, ii_leg)
jj_leg = np.append(j0, jj_leg)
# cast to list
jj_leg, ii_leg = list(jj_leg.astype("int")), list(ii_leg.astype("int"))
return zip(jj_leg, ii_leg)
def section_indices(ji_pairs=None):
"""Forms a staircase of consectuive segments from
ij_pairs[0] to ij_pairs[-1]. Number of segments will be one
less than number of input points.
Parameters
----------
ji_pairs : list of tuples
List of tuples of lenth n that defines section nodes.
`[(j0,i0),(j1,i1), ..., (jn,in)]`
Return
------
ji_pairs : Iterator
Iterator of tuples of length n
"""
ji_pairs_for_sub_sections = map(section_leg_indices,
ji_pairs[:-1], ji_pairs[1:])
ji_pairs_stitched_together = reduce(
lambda x, y: chain(x, islice(y, 1, None)),
ji_pairs_for_sub_sections)
return ji_pairs_stitched_together
# Heading of section legs
if mesh_hgr:
latlon = [(mesh_hgr.gphif.isel(y=j0, x=i0).values.item(),
mesh_hgr.glamf.isel(y=j0, x=i0).values.item()),
(mesh_hgr.gphif.isel(y=j1, x=i1).values.item(),
mesh_hgr.glamf.isel(y=j1, x=i1).values.item())]
a = heading(latlon)
a = [a]*len(jj_leg)
return zip(jj_leg, ii_leg, a)
else:
return zip(jj_leg, ii_leg)
def reduce_dataset(data, vars_to_keep, coords_to_keep):
......@@ -267,7 +393,9 @@ def select_section(ji, gridU, gridV, mesh_hgr=None, mesh_zgr=None, mask=None,
Parameters
----------
ji : Iterator of tuples.
Iterator of n tuples of ji-coordinates defining the secion.
Iterator of n tuples of ji-coordinates defining the section.
If tuples include angle of section leg as last lement this angle is
used to rotate velocites.
gridU : xarray.Dataset
Dataset of shape (t,z,y,x) holding variables at U-point.
Use dataset returned by broken_line.reduce_dataset()
......@@ -307,9 +435,15 @@ def select_section(ji, gridU, gridV, mesh_hgr=None, mesh_zgr=None, mask=None,
"""
jj, ii = list(zip(*ji))
ji = list(zip(*ji))
jj = ji[0]
ii = ji[1]
ii = xr.DataArray(list(ii), dims='c')
jj = xr.DataArray(list(jj), dims='c')
if len(ji) == 3:
a = ji[2]
a = xr.DataArray(list(a[:-1]), dims='c')
iidiff = ii.diff('c')
jjdiff = jj.diff('c')
gridU = xr.concat([gridU, gridU.shift(y=-1)], dim='s').\
......@@ -327,9 +461,37 @@ def select_section(ji, gridU, gridV, mesh_hgr=None, mesh_zgr=None, mask=None,
.where(jjdiff, 0)
gridV = gridV.isel(s=1).where(iidiff == 1, gridV.isel(s=0))\
.where(iidiff, 0)
#
if ((len(ji) == 3) and mesh_hgr):
# Compute and select inclination of coordinate system
a_coord_x =\
coordinate_rotation(mesh_hgr.gphif, mesh_hgr.glamf, 'x')
a_coord_y =\
coordinate_rotation(mesh_hgr.gphif, mesh_hgr.glamf, 'y')
a_coord_x =\
select_auxillary_variables(a_coord_x, jj, ii, iidiff, 'x')
a_coord_y =\
select_auxillary_variables(a_coord_y, jj, ii, jjdiff, 'y')
a_coord = a_coord_x + a_coord_y
# correct angles of section legs by local
# inclination of coordinate system
a = a - a_coord
# rotate velocties
u_rot_gridu, v_rot_gridu =\
rotate_velocities(gridU['u_normal'], gridU['u_along'], a)
u_rot_gridv, v_rot_gridv =\
rotate_velocities(gridV['u_along'], gridV['u_normal'], a)
gridU = gridU.assign({'u_rot_normal': u_rot_gridu,
'u_rot_along': v_rot_gridu})
gridV = gridV.assign({'u_rot_normal': u_rot_gridv,
'u_rot_along': v_rot_gridv})
# Define direction of positiv transport
gridU['u_normal'] = gridU['u_normal'] * jjdiff
gridV['u_normal'] = gridV['u_normal'] * iidiff * (-1)
if (("u_along" in gridU.data_vars) & ("u_along" in gridV.data_vars)):
gridU['u_along'] = gridU['u_along'] * jjdiff
gridV['u_along'] = gridV['u_along'] * iidiff
section = gridU + gridV
section = section.assign({'ii': ii[:-1]})
section = section.assign({'jj': jj[:-1]})
......
......@@ -353,11 +353,14 @@ def test_along_section_velocity_zonal_section():
selected to give the across section velocity. The respective
orthogonal component is interpolated to this point to give the
along section velocity.
Along section velocities are defined positiv in the direction of the
progression of the section
We test a zonal section.
The marked values are used to calculate the along section velocity
at the V-points. We move from right to left.
at the V-points. We move from right to left. Note the reversed sign of
U_along.
Zonal flow field:
U:
......@@ -376,7 +379,7 @@ def test_along_section_velocity_zonal_section():
V = 1
U_along = [(5+6+6+7)/4, (2+3+3+4)/4 (3+4+4+5)/4, (4+5+5+6)/4]
= [6, 5, 4, 3]
= [-6, -5, -4, -3]
"""
u = np.array(list(range(5))*5).reshape(1, 5, 5) +\
np.array((list(range(5)))).reshape(1, 5, 1)
......@@ -396,7 +399,7 @@ def test_along_section_velocity_zonal_section():
mask = xr.Dataset({'umask': (['t', 'z', 'y', 'x'], np.ones((1, 1, 5, 5))),
'vmask': (['t', 'z', 'y', 'x'], np.ones((1, 1, 5, 5)))})
expected_u_along = [6, 5, 4, 3]
expected_u_along = [-6, -5, -4, -3]
ji_pairs = [(2, 4), (2, 0)]
......@@ -467,3 +470,68 @@ def test_along_section_velocity_meridional_section():
result_u_along = section['u_along']
assert (expected_u_along == result_u_along).all(),\
"Calculation of along section velocity failed."
def test_cross_section_velocity():
"""
We rotate the velocities according to the orientation of each seaction leg
to give cross section velocities. Additionally we take the local
inclination of the coordinate system (theta; relative to true north)
into account. Consider the follwing ocean model output:
U = 1 (in i-direction)
V = 0 (in j-direction)
theta = pi/8 rad
We want to compute velocities across a section with an heading of
alpha=3/8*pi (relative to true north) going from southwest to northeast.
The cross section velocity is given by:
U_cross = U * cos(alpha-theta) = cos(1/4*pi)
U_along = U * sin(alpha-theta) = sin(1/4*pi)
"""
# Create rotated coordinate system
theta = np.pi/8
lon0 = np.array(list(range(5))*5).reshape(1, 5, 5)
lat0 = np.array(list(range(5))*5).reshape(1, 5, 5, order='F')
lon = lon0 * xr.ufuncs.cos(-theta) - lat0 * np.sin(-theta)
lat = lon0 * xr.ufuncs.sin(-theta) + lat0 * np.cos(-theta)
mesh_hgr = xr.Dataset({'gphiv': (['t', 'y', 'x'],
lat),
'gphiu': (['t', 'y', 'x'],
lat),
'gphif': (['t', 'y', 'x'],
lat),
'glamv': (['t', 'y', 'x'],
lon),
'glamu': (['t', 'y', 'x'],
lon),
'glamf': (['t', 'y', 'x'],
lon),
'e1v': (['t', 'y', 'x'], np.ones((1, 5, 5))),
'e2v': (['t', 'y', 'x'], np.ones((1, 5, 5))),
'e1u': (['t', 'y', 'x'], np.ones((1, 5, 5))),
'e2u': (['t', 'y', 'x'], np.ones((1, 5, 5)))})
mesh_zgr = xr.Dataset({'e3u': (['t', 'z', 'y', 'x'],
np.ones((1, 1, 5, 5))),
'e3t': (['t', 'z', 'y', 'x'],
np.ones((1, 1, 5, 5))),
'e3v': (['t', 'z', 'y', 'x'],
np.ones((1, 1, 5, 5)))})
mask = xr.Dataset({'umask': (['t', 'z', 'y', 'x'], np.ones((1, 1, 5, 5))),
'vmask': (['t', 'z', 'y', 'x'], np.ones((1, 1, 5, 5)))})
gridU = xr.Dataset({'vozocrtx': (['depthu', 'y', 'x'],
np.ones((1, 5, 5)))})
gridV = xr.Dataset({'vomecrty': (['depthv', 'y', 'x'],
np.zeros((1, 5, 5)))})
ji_pairs = [(0, 0), (4, 4)]
ji = bl.section_indices(ji_pairs, mesh_hgr)
gridU, gridV = bl.shift_grids(gridU, gridV, mesh_hgr, mesh_zgr, mask)
section = bl.select_section(ji, gridU, gridV, mesh_hgr)
assert ((section.u_rot_along.isel(c=slice(0, -1)) - np.sin(np.pi/4))
< 1e-3).all(),\
"Rotation of velocities failed"
assert ((section.u_rot_normal.isel(c=slice(0, -1)) - np.cos(np.pi/4))
< 1e-3).all(),\
"Rotation of velocities failed"
import numpy as np
import xarray as xr
from xorca_brokenline import shift_grids
import pytest
def test_shifted_grids_UV():
......@@ -346,7 +345,6 @@ def test_reduce_datasets():
"Not all variables in processed gridU have the same shape."
@pytest.mark.xfail
def test_volume_weighted_interpolation():
"""
Interpolation on U/V-points hast to take the different volumes of the grid
......
import xorca_brokenline as bl
import numpy as np
import xarray as xr
# import pytest
def test_heading():
"""
The heading [0,360) is defined positive in clockwise direction with zero
indicating true north
Starting at (0,0) and given the following targets:
(10,0),(0,10),(-10,0),(0,-10)
The headings are (0,90,180,270)
"""
expected_values = [0, 90, 180, 270]
ll_pairs = [[(0, 0), (10, 0)], [(0, 0), (0, 10)], [(0, 0), (-10, 0)],
[(0, 0), (0, -10)]]
result = [bl.heading(ll)/np.pi*180 for ll in ll_pairs]
assert (expected_values == result), "Heading yields wrong result"
def test_rotated_velocities():
"""We define rotated velocities at the points of the transport-velocities,
i.e. rotated velocities sit at the same points as u_normal.
We assume theta to be the angle between the section and the y-axis and
rotate the coordinate system clockwise
Thus the rotation matrix is defined as:
R[ theta - pi/2 ] = [ cos theta -sin theta ]
[ sin theta cos theta]
u_rot = u * cos(theta) - v * sin(theta)
v_rot = u * sin(theta) + v * cos(theta)
u_rot gives the velocity perpindicular to the section.
Assume U = V = 1
theta = [0.25*pi, 0.5*pi, 0.75*pi, pi, 1.5*pi,1.75*pi, 2*pi]
u_rot = [0, -1, -sqrt(2), -1, 0, sqrt(2), 1]
v_rot = [sqrt(2), 1, 0, -1, -sqrt(2), -1, 0, 1]
"""
u = 1
v = 1
expected_u_rot = [0, -1, -np.sqrt(2), -1, 1, np.sqrt(2), 1]
expected_v_rot = [np.sqrt(2), 1, 0, -1, -1, 0, 1]
theta = [0.25*np.pi, 0.5*np.pi, 0.75*np.pi, np.pi, 1.5*np.pi, 1.75*np.pi,
2*np.pi]
u_rot, v_rot = bl.rotate_velocities(u, v, theta)
assert (np.round(expected_u_rot, 8) == np.round(u_rot, 8)).all(),\
"Rotation of velocities failed"
assert (np.round(expected_v_rot, 8) == np.round(v_rot, 8)).all(),\
"Rotation of velocities failed"
def test_coord_rotation():
"""
In order to produce correct cross-section velocities the local inclination
of the coordinate system has to be taken into account.
We desing a coordinate system rotated by 45° and use coordinate_rotation()
to recover that value. coordinate_rotation() can calculate the inclination
based on the the u and v interfaces. We test both ways.
"""
# Create rotated coordinate system
theta = np.pi/4
lon0 = np.array(list(range(5))*4).reshape(1, 4, 5)
lat0 = np.array(list(range(4))*5).reshape(1, 4, 5, order='F')
lon = lon0 * xr.ufuncs.cos(theta) - lat0 * np.sin(theta)
lat = lon0 * xr.ufuncs.sin(theta) + lat0 * np.cos(theta)
mesh_hgr = xr.Dataset({'gphiv': (['t', 'y', 'x'],
lat),
'glamv': (['t', 'y', 'x'],
lon)})
coord_rot_x = bl.coordinate_rotation(mesh_hgr.gphiv, mesh_hgr.glamv, 'x')\
- theta
coord_rot_y = bl.coordinate_rotation(mesh_hgr.gphiv, mesh_hgr.glamv, 'y')\
- theta
print(coord_rot_x)
assert (coord_rot_x < 1e-2).all(),\
"Estimate ofcoordinate system inclination is not correct"
assert (coord_rot_y < 1e-2).all(),\
"Estimate ofcoordinate system inclination is not correct"
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment