Issue with vertcross function when using Fortran-order data from .npy files
I encountered an issue when using the vertcross function from the wrf-python package to perform vertical interpolation on data stored in a .npy file. Specifically, when I load the .npy file , the data is stored in Fortran-order (column-major) format by default.
The vertcross function seems to require C-order (row-major) data storage format to work correctly. When passing Fortran-order data, I received a ValueError: invalid shape for argument 1 (got (38,), expected (360,)), which suggests a misinterpretation of the data shape due to the storage format mismatch.
Steps to Reproduce:
-
Load vorticity data counted by wrfout (e.g., vorticity) from a
.npyfile:vorticity = np.load("vorticity.npy") -
Attempt to use vertcross to perform vertical interpolation:
vs_cross = vertcross(vorticity, z, ...) -
Encounter the error: ValueError: invalid shape for argument 1 (got (38,), expected (360,)) The vertcross function should work correctly when the data is in Fortran-order (which is the default when reading .npy files). If the function requires C-order, it would be helpful to include a clear message in the documentation or provide an automatic conversion inside the function.
solution:
Manually converting the data to C-order format using np.ascontiguousarray() resolves the issue:
vorticity = np.ascontiguousarray(vorticity)
vs_cross = vertcross(vorticity, z, ...)
It would be helpful if the vertcross function could automatically handle Fortran-order data or provide better error handling or documentation related to data format expectations.
Behavior noted here: https://github.com/NCAR/wrf-python/blob/eec6637df55b05bc15112d1fb3c5c3917f26872a/src/wrf/extension.py#L178
Switching both arrays to fortran-contiguous z-y-x produces an array of nans:
>>> from wrf import vertcross, CoordPair
>>> vorticity = np.zeros((61, 184, 251), dtype="f4", order="f")
>>> vertical = vorticity.copy()
>>> vertical[:, :, :] = np.arange(61)[:, np.newaxis, np.newaxis]
>>> center = CoordPair(x=120, y=90)
>>> vertcross(np.asfortranarray(vertical), np.asfortranarray(vertical), pivot_point=center, angle=45)
<xarray.DataArray 'field3d_cross' (dim_0: 100, dim_1: 259)> Size: 104kB
array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]], dtype=float32)
Dimensions without coordinates: dim_0, dim_1
Attributes:
orientation: (30.0, 0.0) to (213.0, 183.0) ; center=CoordPair(x=120, y...
missing_value: 9.969209968386869e+36
_FillValue: 9.969209968386869e+36
Do you want to make a PR promoting that comment to a note in the docstring and propagate that note up the call chain? The alternative would be to ensure the array is C-contiguous at the relevant point, but I'm not sure where that is.
Thanks for your response! I'll take some time to look deeper into the code to determine whether it's better to submit a PR that adds a note in the docstring or to modify the code to ensure the array is C-contiguous automatically. This might take some time, but I'll follow up once I have a better understanding.