change isotropize to satisfy Parseval's thereom

[miniufo]

This PR tries to solve #219.

Also, we add 3D detrending codes, which may be useful when calculating 3D (2D space and 1D time) frequency-wavenumber spectrum.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[roxyboy]

Thanks for the PR @miniufo !
I think changing the mean to sum makes sense. Could you run black on your code for formatting and also check whether the tests (here) pass by running pytest locally? It seems that the automatic Github workflow for tests isn't working at the moment...

[miniufo]

Hi @roxyboy, following your suggestions, I've reformatted the modified codes using black. But the assertion of all-close for isotropize at here:

xrft/xrft/tests/test_xrft.py

Lines 953 to 959 in 9ac46f9

	def _test_iso(theta):
	ps = xrft.power_spectrum(theta, spacing_tol=spacing_tol, dim=dims)
	ps = np.sqrt(ps.freq_x**2 + ps.freq_y**2)
	ps_iso = xrft.isotropize(ps, fftdim, nfactor=nfactor, truncate=truncate)
	assert len(ps_iso.dims) == 1
	assert ps_iso.dims[0] == "freq_r"
	npt.assert_allclose(ps_iso, ps_iso.freq_r**2 * 2 * np.pi, atol=0.02)

cannot be passed. The slopes of the two arrays are the same but differ by a constant that may depend on the number of points within each radial band, as sum is used instead of mean. Do you have any suggestion on changing the test codes?

[roxyboy]

Hi @roxyboy, following your suggestions, I've reformatted the modified codes using black. But the assertion of all-close for isotropize at here:

xrft/xrft/tests/test_xrft.py

Lines 953 to 959 in 9ac46f9

	def _test_iso(theta):
	ps = xrft.power_spectrum(theta, spacing_tol=spacing_tol, dim=dims)
	ps = np.sqrt(ps.freq_x**2 + ps.freq_y**2)
	ps_iso = xrft.isotropize(ps, fftdim, nfactor=nfactor, truncate=truncate)
	assert len(ps_iso.dims) == 1
	assert ps_iso.dims[0] == "freq_r"
	npt.assert_allclose(ps_iso, ps_iso.freq_r**2 * 2 * np.pi, atol=0.02)

cannot be passed. The slopes of the two arrays are the same but differ by a constant that may depend on the number of points within each radial band, as sum is used instead of mean. Do you have any suggestion on changing the test codes?

Apologies for my slow reply, I was away for the past week. Thanks for the Black formatting!

I understand that the assertion isn't passing because the radial wavenumber is taken using a mean but the spectrum is isotropized using sum. Since your PR is about making sure the isotropic spectra satisfy the Parseval's theorem, instead of trying to use the wavenumber, we should probably remove line 955 and change the assertion to something like npt.assert_allclose( (ps_iso*ps_iso.freq_r_diff).sum(), ps.sum() ).

[miniufo]

I've changed the assertion to npt.assert_allclose( (ps_iso).sum(), ps.sum() ). Since sum is used azimuthally in isotropize, another radial sum is enough to get the total energy. Adding an attribute ps_iso.freq_r.spacing would be much better, as there are ps.freq_x.spacing and ps.freq_y.spacing. But I find ps_iso.freq_r.diff('freq_r') is not a constant.

I also changed atol to 1 in test_isotropic_ps_slope:

xrft/xrft/tests/test_xrft.py

Line 1012 in 9ac46f9

npt.assert_allclose(a, s, atol=0.1)

Not sure if it is fine. It cannot get passed when atol=0.1.

[roxyboy]

Does this mean that this PR results in the slopes differing by a factor of one...? If so, this is not good.

[roxyboy]

I also changed atol to 1 in test_isotropic_ps_slope:

xrft/xrft/tests/test_xrft.py

Line 1012 in 9ac46f9

npt.assert_allclose(a, s, atol=0.1)

Not sure if it is fine. It cannot get passed when atol=0.1.

Yes, this is what I was personally worried that could happen when switching the isotropization from mean to sum. In a discrete (and not continuous) system as we have for data analyses, the number of points within each radial band differs. A mean isn't sensitive to this but a sum is. The spectral slopes disagreeing by a factor of one is not going to be acceptable unfortunately...

[miniufo]

Hi @roxyboy, I just looked into the test, and found that the frequency range to be fitted is crucial.

xrft/xrft/tests/test_xrft.py

Line 1011 in 9ac46f9

y_fit, a, b = xrft.fit_loglog(iso_ps.freq_r.values[4:], iso_ps.values[4:])

Previously, indices from 4 to end are used, which include the highest frequencies that are steeper than expected. Now I switch the indices from 0 to -35, which overlap with the expected range, and this test passed even with a smaller `atol=0.06.

The results looks like:

The red curve is the fitted line over indices from 0 to -35. You may want to check if this is OK.

[roxyboy]

@miniufo Ok, this is great! Can you run the Black formatting once more and I can merge this :)

[miniufo]

Oops, I forgot to do this after any changes. Can you check it again?

[roxyboy]

@miniufo Sorry, it seems that some files are still missing the Black formatting...

[miniufo]

That's strange. I format all .py files. Possibly missing the .ipynb files. Please try again and see if it is fine now.

[roxyboy]

@miniufo I agree that this is kind of unexpected behavior but the GitHub log is saying that test_xrft-checkpoint.py is missing the formatting...

> Run black --check xrft
  black --check xrft
  shell: /usr/bin/bash -e {0}
  env:
    pythonLocation: /opt/hostedtoolcache/Python/3.8.18/x64
    LD_LIBRARY_PATH: /opt/hostedtoolcache/Python/3.8.18/x64/lib
would reformat xrft/tests/.ipynb_checkpoints/test_xrft-checkpoint.py

[roxyboy]

Thanks for your work @miniufo ! The checks are all passing now so I will merge this.

[miniufo]

That's great~

[miniufo]

Hi @roxyboy, I've recently found that in a 2D homogeneous turbulence, the zonal/meridional spectra are not the same as isotropic spectra. The following codes passed the test of Parseval's thereom:

import xrft
import numpy.testing as npt
import xarray as xr
import numpy as np

ds = xr.open_dataset('/mnt/Data3/2DTurbulence/8192/non/sim_bih_ray_flow_part09.nc', chunks={'time':1})
psi = ds.psi[0].load()  # 8192 * 8192 streamfunction

psihx = xrft.power_spectrum(psi, dim=['x'], scaling='density').mean('y')  # x-component
psihy = xrft.power_spectrum(psi, dim=['y'], scaling='density').mean('x')  # y-component
psihi = xrft.isotropic_power_spectrum(psi, dim=['x', 'y'], scaling='density',
                                      truncate=True)                      # r-component

npt.assert_almost_equal((psi**2 * dx).sum('x').mean('y'),   # variance in grid space
                         psihx.sum('freq_x'), decimal=5)    # variance in spectral space
npt.assert_almost_equal((psi**2 * dy).sum('y').mean('x'),   # variance in grid space
                         psihy.sum('freq_y'), decimal=5)    # variance in spectral space
npt.assert_almost_equal((psi**2 * dx * dy).sum(['x', 'y']), # total variance in grid space
                         psihi.sum('freq_r'), decimal=5)    # total variance in spectral space

but give the plot like:

which is a bit surprising to me that the radial spectra differs from x-/y-components by a constant. The constant could be dkr / (Lx * Ly). Here the important thing is the unit of the result (also determines the magnitude). I am expecting the radial spectra is of the same magnitude as the x-/y-components.

Also, I find the docstring for scaling='density' or scaling='spectrum' is a bit unclear. I feel the difference can be made more clearly by the units of the outputs. I guess one is [A]^2 m^3 while the other is [A]^2 m^2. Do you have any comments on this?

[roxyboy]

Hi @roxyboy, I've recently found that in a 2D homogeneous turbulence, the zonal/meridional spectra are not the same as isotropic spectra.

Yes, they are not the same due to the Jacobian that is needed for coordinate transformation from Cartesian to polar coordinates.

Also, I find the docstring for scaling='density' or scaling='spectrum' is a bit unclear. I feel the difference can be made more clearly by the units of the outputs. I guess one is [A]^2 m^3 while the other is [A]^2 m^2. Do you have any comments on this?

Regarding the density and spectrum docstrings, they are somewhat of a (confusing but commonly accepted) jargon used in the field of FFT and we decided to adopt them.

[miniufo]

Thank you for sharing these. @roxyboy. I was trying to evaluate the spectra from a model simulation agaist SWOT observation. The calculation based on the model outputs uses 2D isotropic spectrum while the one based on SWOT uses along-track 1D power spectrum. Then how to reconcile these if they are not at the same energy level?

[roxyboy]

I'm going to add back the option of taking the azimuthal mean in isotropizing the 2D spectra #225 .
The summation approach is going to introduce strong isotropizing artifacts when the aspect ratio of the Fourier transformed domain is highly uneven.