（Incorrect calculation）Incorrect Meridional Scale (dy) Calculation in parse_grid_arguments for Partial differential calculation

[har20230609]

What went wrong?

When using the projection <Geographic 2D CRS: +proj=latlon +type=crs> with the datum World Geodetic System 1984, the calculation of vorticity from the wind field uv fails to provide accurate results when dx and dy are not explicitly specified. The same issue may occur with other projections as well.Specifically, the parse_grid_arguments decorator in the tool.py file under the calc directory returns an incorrect dy to the vorticity function. The dy is prematurely assigned as meridional_scale, resulting in significant errors in the final output.

Steps to Reproduce
Set the projection to <Geographic 2D CRS: +proj=latlon +type=crs>.
Attempt to compute vorticity using the wind field uv data.
Do not explicitly provide dx and dy.
Check the output from the vorticity function for accuracy.

Actual Result
The vorticity calculation yields inaccurate results due to the incorrect dy value. The param dy is prematurely assigned as meridional_scale

Code Trace
The issue originates in the parse_grid_arguments section related to metpy.grid_deltas, as shown below:

grid_deltas = grid_prototype.metpy.grid_deltas
bound_args.arguments['dx'] = grid_deltas['dx']
bound_args.arguments['dy'] = grid_deltas['dy']

Further tracing leads to the nominal_lat_lon_grid_deltas section:

forward_az, _, dy = geod.inv(lon_meridian_diff, lat[:-1], lon_meridian_diff, lat[1:], radians=False)

Operating System

Windows

Version

1.5.0

Python Version

3.8

Code to Reproduce

import metpy.calc as mpcalc
from metpy.cbook import example_data

# load example data
ds = example_data()
u,v = ds.uwind,ds.vwind
# Calculate the vertical vorticity of the flow
a = mpcalc.vorticity(u, v)

Errors, Traceback, and Logs

[dcamron]

Hello, thanks for opening your first issue with us. Please provide these significant errors in final value and what values and methods you are comparing against, and please elaborate on what errors you believe to be tracing through grid_deltas. Factors like meridional_scale are considered separately from deltas in our projection processing, and I don't trace any errors in dy through your provided example.

In case it helps, there are thorough discussions on these implementations and statistical comparison of their outputs to other sources in #893 and #2743. Please follow up with any additional questions plus the clarifications above so that we can understand your issue and help. Thanks!

[har20230609]

Hello, thanks for opening your first issue with us. Please provide these significant errors in final value and what values and methods you are comparing against, and please elaborate on what errors you believe to be tracing through grid_deltas. Factors like meridional_scale are considered separately from deltas in our projection processing, and I don't trace any errors in dy through your provided example.

In case it helps, there are thorough discussions on these implementations and statistical comparison of their outputs to other sources in #893 and #2743. Please follow up with any additional questions plus the clarifications above so that we can understand your issue and help. Thanks!

thanks for your reply
I understand the calculation process discussed in #893. I think that in the coordinate system talked above , the actual projection process should involve projecting spherical coordinates onto ellipsoidal coordinates. This is the basis for my def get_delta_factor below, which references Chapter 5 of the book at https://doi.org/10.1007/978-3-642-41245-5.

From the code below , we can see that MetPy encounters errors when calculating delta in this simple coordinate system. Additionally, when considering the ellipsoidal coordinate system, 𝑑𝑥 should be the semi-minor axis 𝑏 multiplied by dλ instead of the semi-major axis 𝑎 multiplied by dλ. Therefore, to be precise, there is also an issue with dx in MetPy. However, the numerical impact of this error is relatively small, which may make it difficult to detect.

Here is the code that demonstrates the problem

import numpy as np
from metpy.cbook import example_data

def get_delta_factor(lon, lat):
    a = 6378137.0
    e2 = 6.69437999014132e-3
    b = a * np.sqrt(1 - e2)

    dfi = np.diff(lat) / 180 * np.pi
    dlamda = np.diff(lon) / 180 * np.pi

    dy = a * dfi
    dx = b * dlamda

    m_scale = (1 - e2 * (np.sin(np.radians(lat)) ** 2)) ** (3 / 2) / (1 - e2)
    p_scale = np.sqrt(1 - e2 * (np.sin(np.radians(lat)) ** 2)) / np.cos(np.radians(lat))

    return dx, dy, p_scale, m_scale
if __name__ == "__main__":
    ds = example_data()
    u,v,lon,lat = ds.uwind,ds.vwind,ds.lon,ds.lat
    # calculate the param
    # 1.Here's what I think the calculations are
    dx, dy, p_scale, m_scale = get_delta_factor(lon, lat)
    # 2.Here's the delta you can trace in calculating vorticity from metpy
    delta_metpy = u.metpy.grid_deltas
    # 3.And You will notice that the dy used in metpy has been scaled
    dy_scaled = dy/m_scale[:-1]
    print("dy_scaled:"+" "*5 +"delta_metpy['dy']:")
    for value1,value2 in zip(dy_scaled.data,delta_metpy['dy'].magnitude):
        print(f"{value1:.5f},  {value2:.5f}")
    # Since the length of m_scale is one more than the total number of dy, the first n-1 values are simply divided by dy.
    # This may slightly affect the value of dy_metpy.
    # However, we can still observe that  dy_scaled is exactly the dalta_metpy['y'], the result after dy has been appropriately adjusted.
    # Moreover, the values of p_scale and m_scale that I calculated are consistent with the 
    # parallel_scale and meridional_scale obtained from MetPy, though I haven't included those here.

looking forward to your response. Thanks for your time!

[dcamron]

Perfect, this is a super helpful follow-up! We are prioritizing v1.7 releasing this week, but will take a look at this as soon as we can. Thanks.