add rtol for calculating grid spacing

[zfan001]

I had several running errors during running neighbourhood processing and regridding.

"ValueError: Coordinate longitude points are not equally spaced"

It is from calculate_grid_spacing function. It happened for several BOM-commonly-used grids.
Grid spacing is from np.diff(coord). I think that it is reasonable to add rtol=7.0e-5 which made all failure grids passed.

Testing:

• Ran unit and acceptance tests and they passed OK

[codecov]

Codecov Report

Merging #1545 (83bdfbb) into master (71fd210) will not change coverage.
The diff coverage is 100.00%.

@@           Coverage Diff           @@
##           master    #1545   +/-   ##
=======================================
  Coverage   97.93%   97.93%           
=======================================
  Files         106      106           
  Lines        9447     9447           
=======================================
  Hits         9252     9252           
  Misses        195      195           

Impacted Files	Coverage Δ
improver/utilities/spatial.py	98.66% <ø> (ø)
improver/regrid/grid.py	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 71fd210...83bdfbb. Read the comment docs.

[gavinevans]

Thanks for this PR, @zfan001 👍

I've suggested some modifications to this PR in the comments below. Primarily I've recommended to pass in the relative tolerance that you want from the CLI, rather than change the relative tolerance in the function calls. This is particularly the case for the centralised utilities, which are used in a number of different places and therefore it is hard to gauge the impact of the change to the relative threshold. There might also be a concern that this change could mask unintended changes to our grid spacing that would have previously caused an error. These changes will increase the size of the PR but will hopefully make this change more appropriate for different users.

[tjtg]

There are a few things going on here. I'll try and pull it apart into separate questions/issues, then give my opinion on solutions to those.

Some functionality in IMPROVER expects equally spaced or equal area grids. This is a useful simplifying assumption about the input which makes calculations simpler to implement and faster to run, as opposed to doing all calculations in a grid-generic way.

Equal spacing is determined by the calculate_grid_spacing function. This function is used in regrid, pysteps advection and optical flow.
Equal area is determined by the check_if_grid_is_equal_area function. This function is used in circular kernel, pysteps advection and spatial weighting.

Questions:

1. How strict should the check for an equally spaced grid be?
   • This is currently implemented by the rtol argument to calculate_grid_spacing, with a default value of 0.0 (no tolerance, grid spacings must be exactly identical even though they are floats)
2. What is meant by an equal area grid?
   • This is currently implemented by checking if the X and Y dimensions both have equal spacing in units of metres. I suspect that distance as converted to units of metres is measured in the 2D projection rather than an on-the-ground spherical or elliptical distance, but am not sure of this.
3. How strict should the check for an equal area grid be?
   • This is a related question to the first point above.

My opinions on those questions:

1. The check for equal spacing should allow for floating point rounding, plus a small amount of extra tolerance. See example from Add tolerance in grid even-spacing checking & calculation #1459 (comment) for more detail on why some tolerance is needed and an example. Spatial coordinates are often using float32, so float32 epsilon (~=1.19e-7) is the minimum plausible relative tolerance. I think a little more tolerance is OK and unlikely to cause problems, so perhaps1e-6 relative tolerance as a default for calculate_grid_spacing.
   • An exception to that: the use of calculate_grid_spacing in regrid should have a looser tolerance to allow for regridding of input data that has arrived from elsewhere with some mangling of coordinates (eg. poorly done conversion between formats). From memory, the nearest/bilinear regrid code does use the assumption that the input grid is equally spaced to index into the data, so the tolerance shouldn't be overly slack. The 7e-5 relative tolerance value from this PR could be OK here.
2. The conditions in check_if_grid_is_equal_area as currently implemented aren't a check for equal area. For example, data using a cylindrical equal area projection would not pass the check as the grid spacing in the latitude dimension is not equal. However, I think the check being made does make sense for the intended purpose (eg. advection is safe/straightforward in both X and Y spatial dimensions). I'd suggest renaming the function and docstring to describe what's actually happening - maybe something like "grid is spatially square"?
3. Since this check is built on calculate_grid_spacing I'd suggest a similar default tolerance eg. floating point tolerance plus a bit as described in (1) above.

I don't think exposing these tolerances in the CLI is needed if we can get agreement on sensible values. I suspect that some of the coordinates metadata in files @zfan001 has been using as input are problematic (as in, differences are larger than unavoidable float rounding differences).

Implementing that would mean:

• Close this PR
• Create a new PR that does the following:
  • Set rtol default in calculate_grid_spacing to something non-zero
  • Rename the check_if_grid_is_equal_area and provide a better description in its docstring
  • Where calculate_grid_spacing is used in regrid, set a larger tolerance

I'd suggest its worthwhile getting another person's opinion on this before proceeding as there are some complex questions about meanings and intention here. Perhaps @gavinevans or @MoseleyS on the MO side or @thbom001 on the BOM side?

[gavinevans]

Thanks for the detailed write-up, @tjtg. I've pointed @MoseleyS at this, as I'm unfamiliar with the details of some of the usages of the calculate_grid_spacing function, so couldn't provide a recommendation on what might be a sensible default for the relative tolerance. I've pointed this at @bjwheltor for guidance in terms of the renaming of the check_if_grid_is_equal_area function. If a sensible non-zero default for the calculate_grid_spacing function that meets @MoseleyS's requirements can be chosen then that would be the simplest way forward, as you've proposed.

[zfan001]

@tjtg @gavinevans Thank you for these very good comments and suggestions.

I add more explanation on my comment "There is an inherent problem in the grid-spacing rtol, which is related to dtype(float32) and grid resolution".

For example, generate an equal-space lat/lon grid : longitude: 160:170degree, and latitude: 35:45 degree. grid spacing: 1/3 degree. create grid cube using float32.
Longitude:
array([160. , 160.33333, 160.66667, 161. , 161.33333, 161.66667,
162. , 162.33333, 162.66667, 163. , 163.33333, 163.66667,
164. , 164.33333, 164.66667, 165. , 165.33333, 165.66667,
166. , 166.33333, 166.66667, 167. , 167.33333, 167.66667,
168. , 168.33333, 168.66667, 169. , 169.33333, 169.66667,
170. ], dtype=float32)
With this grid cube and rtol=3.0e-5, calculate_grid_spacing still generates failure.
If changing gap spacing from 1/3 to 1/30, calculate_grid_spacing will fail even with rtol =2.oe-4.

If the grid spacing has recurring decimals, the problem will probably turn up (e.g. BOM' MSAS grid ).
Another NWP grid (grib file) which failed calculate_grid_spacing is probably due to possible accuracy loss in the data format conversion somewhere.

I am happy to implement what Tom suggested next week, if there is no objection. I recommend to use default rtol=1.e-5.

[bjwheltor]

I thought I'd respond to a couple of the earlier comments, although I should make clear that I do not know the details of the code, so I'm responding the the information provided in those comments only.

I agree with the breakdown and assessment by @tjtg and with the specific proposals around implementing a (fairly tight) tolerance - although I'd want @MoseleyS to confirm this does not impact the nowcast - the re-naming / clarification of the check_if_grid_is_equal_area, although I'm not sure exactly how this should be renamed (it appears to do an approximation of what the names suggests, with various assumptions, so I'd want to check the code in detail before addressing this).

The last comment by @zfan001 slightly concerned me, but I may have misunderstood. This concerns latitude-longitude (Plate Carrée projection) grid for which is reasonable to check for equal grid spacing (to a fairly strict tolerance), but not for an equal area grid (which I assume is the intent of check_if_grid_is_equal_area ), as this type of grid is not an equal area grid and so should fail this check. I believe that this type of constraint is usually applied to preserve areal extends of precipitation, without the need to do more complex calculations of the actual grid box areas.

[MoseleyS]

We expected this error to come up some time, and I agree that the solution is to specify an rtol value. Ideally, we should specify the smallest rtol value that achieves what we want. For the Met Office grids, which are all defined as having intervals expressable as non-recurring base-10 values, a tolerance of zero is fine. I agree that for an interval of 1/3 degrees, this is not fine. However, a tolerance of 1e-5 is much too large. 1e-7 is more reasonable (see https://stackoverflow.com/questions/39134956/accuracy-of-float32 and remember that our values are in the range 1e0 to 1e+6).

There is a second issue entangled here, which really ought to be considered separately - that of when a grid is equal-area. As @tjtg pointed out, the method is specific to the only use-case we have so far, which is to detect a Lambert Azimuthal Equal Area grid with equal spacing in x and y coordinates, while excluding a regular Plate Carree grid, which is regular in x and y, but is not equal area. The assumption that a coordinate that can be expressed in metres is acceptable is probably not safe.
Our neighbourhood processing code is deliberately designed to reject grids that are not equal area, because this is assumed in the method we have chosen. For example, if you tried to use neighbourhood processing with a regular lat-lon grid spanning the whole globe, data near the poles would be stretched in y much more than in x. This can be solved by having a latitude-dependent neighbourhood kernel, but we haven't needed it yet, so it is not implemented.

The example of the 1/3 degree grid is very helpful in understanding the purpose of this PR. It ought to be included in a unit-test so that we don't lose sight of it in the future.

[zfan001]

@bjwheltor @MoseleyS Thank you very much for your comments.
The PR is for mainly related to "calculate_grid_spacing" function. The problem I mentioned applies for any projections (plate-Carree, Lambert Azimuthal Equal Area, Albers-equal-area). For example, as long as the recurring decimal spacing and float32 is used, grids with any projections will possibly fail " calculate_grid_spacing".

np.allclose is used to check equal spacing in " calculate_grid_spacing".
https://numpy.org/doc/stable/reference/generated/numpy.allclose.html

Its original default is rtol=1.e-5.
"numpy.allclose(a, b, rtol=1e-05, atol=1e-08, equal_nan=False)".

"If the following equation is element-wise True, then allclose returns True.
absolute(a - b) <= (atol + rtol * absolute(b))"

I recommended rtol=1.0e-5 because coord(float32) has 8 significant digits, but grid spacing will normally lose 3 significant digits in np.diff(). In my opionon, It should eliminate vast majority of problems without causing other problems.

I support that "check_if_grid_is_equal_area" issue is treated separately, though it calls "calculate_grid_spacing". It is seemingly a function only for equal-area-project grids.

To proceed this PR, I suggest to remove all of my changes int/improver/improver/utilities/spatial.py and only change default rtol to 1.0e-5 (or 1.0e-6 ?,1.0e-7 ?) in " calculate_grid_spacing".

[zfan001]

@tjtg @gavinevans @MoseleyS @bjwheltor I changed the code. Basically three are 3 changes in this PR.
(1) change default rtol=1.0e-6 in "calculate_grid_spacing".
(2) in new regridding, loosen rtol=7.0e-5.
(3) add two unit tests with a recurring decimal spacing.

Unit tests/acceptance tests passed OK.

Please review it and let me know if you want any change. Thanks.

[MoseleyS]

Thanks @zfan001
I was wrong about specifying a tolerance of 1e-8 as I had forgotten that some of the precision is taken up by the actual longitude of the domain (e.g. 152), so I now agree that a value around 1e-5 is appropriate.
There's still a couple of things to discuss, but I think we're nearly there.

[bjwheltor]

Thanks, @zfan001 and @MoseleyS for the further comments. I think I understand the concerns a bit better now and see why the tolerance change is required.

[zfan001]

@MoseleyS @gavinevans changes done as suggested. Thanks.

Concerns have been overcome and others have approved the result. Gavin is happy.

[zfan001]

@MoseleyS @tjtg @gavinevans Many thanks!