Initial code for wind field reconstruction from vorticity and divergence

[winash12]

Description Of Changes

Checklist

• Complete draft of the code for the issue - Reconstructing the wind field from vorticity and divergence over finite domains #3207
• Tests are being added

[DWesl]

So I spoke to Sebastian Schemm - the original author of the paper and this is what he had to say and I quote "Two ways to solve this:

1. The paper says «the spherical grid of the employed limited-area model is mapped onto a flat Cartesian grid using a polar stereographic projection.». For that I use a standard routine.

For small enough areas, a map projection would work. Conformal projections (stereographic (polar or otherwise), Lambert Conformal, or Mercator (regular, transverse, or oblique)) will ensure the winds are in the right directions after mapping back to the sphere.

2. The Green’s function simply requires the distance between two grid points, in the paper called r, and the zonal distance x-x’ (or meridional), on the sphere. You compute those distances on the sphere using the great circle distance.

The metric terms in derivatives on the sphere might cause problems, but that affects primarily the tangential distances, not the radial ones, so it might work out.

[winash12]

So I spoke to Sebastian Schemm - the original author of the paper and this is what he had to say and I quote "Two ways to solve this:

1. The paper says «the spherical grid of the employed limited-area model is mapped onto a flat Cartesian grid using a polar stereographic projection.». For that I use a standard routine.

For small enough areas, a map projection would work. Conformal projections (stereographic (polar or otherwise), Lambert Conformal, or Mercator (regular, transverse, or oblique)) will ensure the winds are in the right directions after mapping back to the sphere.

2. The Green’s function simply requires the distance between two grid points, in the paper called r, and the zonal distance x-x’ (or meridional), on the sphere. You compute those distances on the sphere using the great circle distance.

The metric terms in derivatives on the sphere might cause problems, but that affects primarily the tangential distances, not the radial ones, so it might work out.

@DWesl So how would we move forward ? I think theoretically both approaches can be implemented.

With regards to the projection would the function interface need to accept a particular projection ? Then the vorticity and the distance would also need to be calculated with that projection.

Also each calculation of the r^2 in a particular projection would need to include the map scale factors.

In this notebook - https://github.com/evans36/miscellany/blob/main/Vorticity%20and%20Divergence%20Inversion.ipynb

xdiff = (i-x1)*dx[y1,x1].magnitude ydiff = (j-y1)*dy[y1,x1].magnitude

thei-x1 is one of those map scale factors.

So the map scale factors for the distance calculation would need to be part of the interface ?

I have invited @evans36 to this discussion.

[sgdecker]

You'll want to make changes so all the tests pass to get a formal review of this pull request. These are mainly style issues the linter is complaining about. The "Build Docs" checks are failing because there is no docstring in your example calculation.

More substantively, I notice you import xarray, which I believe no other part of the MetPy codebase does. I'm not sure what the Unidata maintainers will think of that.

[DWesl]

So I spoke to Sebastian Schemm - the original author of the paper and this is what he had to say and I quote "Two ways to solve this:

1. The paper says «the spherical grid of the employed limited-area model is mapped onto a flat Cartesian grid using a polar stereographic projection.». For that I use a standard routine.

For small enough areas, a map projection would work. Conformal projections (stereographic (polar or otherwise), Lambert Conformal, or Mercator (regular, transverse, or oblique)) will ensure the winds are in the right directions after mapping back to the sphere.

2. The Green’s function simply requires the distance between two grid points, in the paper called r, and the zonal distance x-x’ (or meridional), on the sphere. You compute those distances on the sphere using the great circle distance.

The metric terms in derivatives on the sphere might cause problems, but that affects primarily the tangential distances, not the radial ones, so it might work out.

@DWesl So how would we move forward ? I think theoretically both approaches can be implemented.

The version you have is already designed to take projected data, just note in the documentation that that is the case, and use one of the metpy or pint wrappers to ensure that the relevant arguments are passed as meters.

With regards to the projection would the function interface need to accept a particular projection ? Then the vorticity and the distance would also need to be calculated with that projection.

Which projection is best depends on a lot of things and would be best left to the user to decide.

Also each calculation of the r^2 in a particular projection would need to include the map scale factors.

In this notebook - https://github.com/evans36/miscellany/blob/main/Vorticity%20and%20Divergence%20Inversion.ipynb

xdiff = (i-x1)*dx[y1,x1].magnitude ydiff = (j-y1)*dy[y1,x1].magnitude

thei-x1 is one of those map scale factors.

I suspect i - x1 and j - y1 are distances in index space, and dy and dx are the map scale factors (or depending on your definition, dx/dx.mean() would be the map scale factors)

So the map scale factors for the distance calculation would need to be part of the interface ?

The very simplest way would be scalar dx and dy as part of the interface, forcing the vorticity and divergence to be on a regular rectangular grid.

The next simplest way would be having x and y coordinates as part of the interface, and using those for the distance calculations. This also lets the user choose a projection that minimizes the distortion in the region they care about without you having to deal with that.

That is, xdiff = (x - x1).magnitude; ydiff = (y -y1).magnitude, although I would suggest leaving the units there so the multiplication two lines later gets the result units right, or using a pint or metpy decorator to force the input units to be 1/s and m so the output will be m/s.

Next most complicated would be averaging the grid scales at each end of the distance calculation. This should give decent results as long as the domain is mostly on the same side of the equator.

xdiff = (i - i1) * 0.5 * (dx[i, j] + dx[i1, j1]); ydiff = (j - j1) * 0.5 * (dy[i, j] + dy[i1, j1])

After that, you get into pyproj geodesic calculations (solution two above), and I think at that point it would be faster to use the pyspharm getuv global spectral functions. This would be the only method that would handle the date line properly.

[winash12]

You'll want to make changes so all the tests pass to get a formal review of this pull request. These are mainly style issues the linter is complaining about. The "Build Docs" checks are failing because there is no docstring in your example calculation.

More substantively, I notice you import xarray, which I believe no other part of the MetPy codebase does. I'm not sure what the Unidata maintainers will think of that.

Thanks Prof. for your first review of this PR. Regarding the import of xarray it's not absolutely necessary. It is required to initialize these 2 data arrays

upsi = xa.zeros_like(umask) vpsi = xa.zeros_like(vmask)

Both the function interfaces have not been agreed upon by all. If the user can create these dataarrays and pass them in the calling interface that would be fine. That is why the review is important. None of the MetPy contributors have as yet looked at this PR.

[DWesl]

More substantively, I notice you import xarray, which I believe no other part of the MetPy codebase does. I'm not sure what the Unidata maintainers will think of that.

MetPy provides an accessor for XArray Datasets and DataArrays, so it gets used somewhere.
In metpy.calc, they have decorators that turn xarray.DataArrays into pint.Quantitys on entering the function and turn the result back into a DataArray on exiting the function. Most functions therefore seem written to take and return pint.Quantity. There is a pint decorator that would convert the function inputs into specified units and provide the units for the output, so the function could be written to take and return NumPy arrays.

Thanks Prof. for your first review of this PR. Regarding the import of xarray it's not absolutely necessary. It is required to initialize these 2 data arrays

upsi = xa.zeros_like(umask) vpsi = xa.zeros_like(vmask)

Would output_u = input_vorticity * dx give the correct size and datatype? If passed instances of pint.Quantity, it should also provide the correct units.

[winash12]

More substantively, I notice you import xarray, which I believe no other part of the MetPy codebase does. I'm not sure what the Unidata maintainers will think of that.

MetPy provides an accessor for XArray Datasets and DataArrays, so it gets used somewhere. In metpy.calc, they have decorators that turn xarray.DataArrays into pint.Quantitys on entering the function and turn the result back into a DataArray on exiting the function. Most functions therefore seem written to take and return pint.Quantity. There is a pint decorator that would convert the function inputs into specified units and provide the units for the output, so the function could be written to take and return NumPy arrays.

Thanks Prof. for your first review of this PR. Regarding the import of xarray it's not absolutely necessary. It is required to initialize these 2 data arrays
upsi = xa.zeros_like(umask) vpsi = xa.zeros_like(vmask)

Would output_u = input_vorticity * dx give the correct size and datatype? If passed instances of pint.Quantity, it should also provide the correct units.

@DWesl I am trying to fix the linting and flake8 errors this week so that I can be set up for a formal code review.

I am stuck with the flake8 and linting errors as for each error I fix and another 10 pop up. On the example code I checked in apparently cfgrib is not an engine I can use because eccodes is not used.

[DWesl]

I am stuck with the flake8 and linting errors as for each error I fix and another 10 pop up.

Yeah, that happens sometimes. If the problem is making the changes locally then waiting for GitHub to provide test results, you can look through the workflows in .github/workflows (the file you want probably has "lint" in the name and will end with either .yml or .yaml) to see what it installs and what commands it runs after, so you can install those packages and run those commands locally until they succeed.

If the problem is it reporting lots of little nitpicky things, you could copy the new functions to a new file, run an auto-formatter on it (black is popular and checks that the code produces the same AST after it's done as before it started, autopep8, yapf, and ruff are other options), then copy the functions from the resulting file back to where you put them originally. You will want to make sure the tests still pass before committing, but this should provide a decent first pass over the code. (Running the auto-formatters on the original file is likely to change lots of other functions you're not otherwise touching, and bigger changes take longer to review).

On the example code I checked in apparently cfgrib is not an engine I can use because eccodes is not used.

I suspect that means eccodes is not installed by the existing infrastructure. You could add cfgrib/eccodes to the list of things they install, or you could convert the grib file to netCDF or CSV and use that in the example instead.

[winash12]

I am stuck with the flake8 and linting errors as for each error I fix and another 10 pop up.

Yeah, that happens sometimes. If the problem is making the changes locally then waiting for GitHub to provide test results, you can look through the workflows in .github/workflows (the file you want probably has "lint" in the name and will end with either .yml or .yaml) to see what it installs and what commands it runs after, so you can install those packages and run those commands locally until they succeed.

If the problem is it reporting lots of little nitpicky things, you could copy the new functions to a new file, run an auto-formatter on it (black is popular and checks that the code produces the same AST after it's done as before it started, autopep8, yapf, and ruff are other options), then copy the functions from the resulting file back to where you put them originally. You will want to make sure the tests still pass before committing, but this should provide a decent first pass over the code. (Running the auto-formatters on the original file is likely to change lots of other functions you're not otherwise touching, and bigger changes take longer to review).

On the example code I checked in apparently cfgrib is not an engine I can use because eccodes is not used.

I suspect that means eccodes is not installed by the existing infrastructure. You could add cfgrib/eccodes to the list of things they install, or you could convert the grib file to netCDF or CSV and use that in the example instead.

Yes the lint , ruff and flake8 errors are all very nitpicky. The thing is I have pep8 and pycodestyle, pylint (not ruff I could not install it) and when I run my changes through these 3 I get no errors at all. But when it goes through MetPy I get a gazillion errors - some of them very trivial and each one I fix brings about others.

yes I will convert my GRIB file to netCDF as it's not worth modifying what MetPy will need as part of it's installation.

[DWesl]

Interesting. My experience with pylint is it finds more things than flake8, which implies they have a few plugins. I don't remember offhand whether flake8 runs pycodestyle checks by default or only if you install a plugin, but I do remember the pep8 python package changed its name to pycodestyle a few years back to avoid confusion with PEP 8.

It looks like the linting requirements are pure python aside from ruff's dependency on Rust, which is rather less ubiquitous than the usual C, C++, or even Fortran dependencies, so installing everything between the first and last lines of the file should get you a flake8 installation that will match the one doing the checks.

For ruff, given the newness of Rust, I am somewhat surprised they do not have CI just make the changes on your PR or as a PR against the branch in your PR. You should still be able to put your functions in a new file, run black or yapf on them, then run flake8 on it until it passes (black is pure python and might have inspired ruff, I think yapf is also pure python), then copy the functions back where you got them and fix the last few flake8 issues with those files.

EDIT: On looking at the current errors, you are past most of those easy problems. The currently-flagged problems are extra spaces in and around docstrings, and writing docstring summaries in declarative rather than the imperative mood (removing the "s" at the end of the verb should fix this).

[DWesl]

A few suggestions for the lint and coverage problems, as requested on winash12#2 here. The comments might make more sense when viewed on the "Files Changed" tab of the PR.

[DWesl]

What problems are you experiencing with np.mgrid and np.meshgrid that you cannot use those? They look like they do something similar.

[DWesl]

Suggested change

	"""Tests the vectorized array indices for two different bounding boxes. """
	"""Tests the vectorized array indices for two different bounding boxes."""

[DWesl]

Suggested change

	""" Test the mask of a bounding box. """
	"""Test the mask of a bounding box."""

[DWesl]

Suggested change

	"""Tests the bounding box indices for 2 different cases. """
	"""Tests the bounding box indices for 2 different cases."""

[DWesl]

Suggested change

	"""Returns the vectorization indices for inner for loop in the
	wind field reconstruction method.
	"""Return the vector indices for the wind field reconstruction method.

It doesn't like the "s" on the verb, and thinks the summary is too long.

[DWesl]

Suggested change

	"""Returns the array indices of a bounding box."""
	"""Return the array indices of a bounding box."""

It wants imperative (telling the computer to do something) rather than declarative (saying what the computer does)

[DWesl]

This is complaining because you started the setup for the tests but didn't add the tests for these functions. Adding the tests or removing the import will make the warning go away. Adding the tests should increase the coverage Codecov reports for the functions.

[DWesl]

They'll want you to change these lines back before the merge; among other things, they cause the tests to fail with a SystemExit error rather than passing.