oceanspy
oceanspy copied to clipboard
hainegroup
Survey velocities don't match exactly original curvilinear velocities
@renskegelderloos The survey aligned velocities don't exactly match the original velocities for curvilinear grids. It might be that there's a small error due to the extra operation with sin and cos from model output. Let me know if you get results that don't look right.
This is where I noticed it: https://github.com/malmans2/oceanspy/blob/959a2488519008f50804325a0a69925304217fdf/oceanspy/tests/test_compute_functions.py#L406-L448
My understanding of this issue is that, this is a problem with the test rather than anything wrong with the code. The test is doing the following things:
- calculate the geographical_aligned velocity. Here is a test showing it is doing the correct thing.
U = od_in._ds['U']
V = od_in._ds['V']
U = od_in._grid.interp(U, "X", boundary="fill", fill_value=np.nan)
V = od_in._grid.interp(V, "Y", boundary="fill", fill_value=np.nan)
od_in = od_in.compute.geographical_aligned_velocities()
o1 = np.hypot(U,V)
o2 = np.hypot(od_in._ds['U_zonal'],od_in._ds['V_merid'])
np.allclose(o1,o2)# True
The result here is True, which means everything is working fine.
- Making a section, and then calculate the section aligned velocity from the data on the section.
s1 = np.hypot(od_surv._ds['U'],od_surv._ds['V']) # "original" data
s2 = np.hypot(od_surv._ds['U_zonal'],od_surv._ds['V_merid']) # zonal/meridional vel
s3 = np.sqrt(ds_out['tan_Vel']**2 + ds_out['ort_Vel']**2) # section aligned
np.allclose(s1,s3) ## False
That's more or less the original test. The current speed does not match. But s3 = s2, and if you do
#correspond to the case when rotate = False
od_surv._ds = od_surv._ds.drop(['U_zonal','V_merid'])
ds_out = survey_aligned_velocities(od_surv)
s4 = np.sqrt(ds_out['tan_Vel']**2 + ds_out['ort_Vel']**2)
np.allclose(s4,s1)## True
you also get s4 = s1. This suggest that the compute function is alright.
I can briefly explain what I think is going on here. The survey_aligned_vel function always prefer to use U_zonal. And subsample.survey_stations only interp the components of velocity, rather than their direction and strength, which actually make sense. say we have $$U_1^2+V_1^2 = u_1^2+v_1^2$$ $$U_2^2+V_2^2 = u_2^2+v_2^2$$ where u represent the original velocity on c-grid, U represent the transformed velocity. After xesmf interpolation: $$U = w_1U_1+w_2U_2$$ $$V = w_1V_1+w_2V_2$$ $$u = w_1u_1+w_2u_2$$ $$v = w_1v_1+w_2v_2$$ Then when you are interested in the current speed at the survey station, you got two different result $$s_1 = u^2+v^2 = w_1^2(u_1^2+v_1^2)+w_2^2(u_1^2+v_1^2)+2w_1w_2(u_1u_2+v_1v_2)$$ $$s_2 = U^2+V^2 = w_1^2(U_1^2+V_1^2)+w_2^2(U_1^2+V_1^2)+2w_1w_2(U_1U_2+V_1V_2)$$ The first two terms are the same, but the last term does not necessarily match.
I also don't think we need to change this current test, because it works. And you can argue that assert_almost_equal is a even better test.
So, in summary, you recommend changing the test to assert_almost_equal and closing?
the test is already assert_almost_equal, I suggest we don't change anything. It might be nice if we can refer to this issue in the comment.
Change
# TODO: not exact match. Check with Renske!
to
# see issue #70 for more info.
, so that less people will bother Renske in the future hhhh.