an instance of the geographic vorticity test case

[winash12]

I am trying to get a geographic vorticity test case going and I get the following error. But first the example code

`
import xarray as xr
import numpy as np
import metpy.calc as mpcalc
from metpy.units import units
import pyproj
import sys
def main():
setUpData()

def setUpData():

crs = pyproj.CRS(4326)
latitude = np.linspace(-90.,90.,721,endpoint=False)*units.degree
longitude = np.linspace(0.,360. ,1440,endpoint=False)*units.degree
u = np.empty((721,1440))
v = np.empty((721,1440))
longitude,latitude = np.meshgrid(longitude.m_as('degree'),latitude.m_as('degree'))
u = - np.sin(latitude/2)*units('m/s')
v = - np.cos(longitude)*units('m/s')
u = xr.DataArray(u,coords={'longitude':(['longitude','latitude'],
                                        longitude),
                           'latitude':(['longitude','latitude'],
                                       latitude)})

v = xr.DataArray(v,coords={'longitude':(['longitude','latitude'],longitude),
                           'latitude':(['longitude','latitude'],latitude)})

vort = mpcalc.vorticity(u,v,longitude,latitude,crs)

main()`

I get this error -

File "/home/winash12/testpython/test_rot_inv.py", line 28, in setUpData vort = mpcalc.vorticity(u,v,longitude,latitude,crs) File "/usr/local/lib/python3.9/dist-packages/metpy/calc/tools.py", line 1032, in wrapper bound_args = newsig.bind(*args, **kwargs) File "/usr/lib/python3.9/inspect.py", line 3062, in bind return self._bind(args, kwargs) File "/usr/lib/python3.9/inspect.py", line 2988, in _bind raise TypeError( TypeError: too many positional arguments

[DanielAdriaansen]

What version of MetPy are you using?

This looks like the problem:
too many positional arguments

Looking at how you're calling the function, you're providing 5 positional arguments but the function only takes u and v as positional arguments, and longitude, latitude, and crs need to be provided as keyword arguments I think like this:
vort = mpcalc.vorticity(u,v,longitude=longitude,latitude=latitude,crs=crs)

Note this is for MetPy v1.6 here: https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.vorticity.html

[winash12]

Yes you are correct - https://github.com/Unidata/MetPy/blob/main/tests/calc/test_kinematics.py#L70 and I am using 1.6. If I use GFS data with just the positional arguments with u and v it's fine. Here I am using mock up data.

My version of the mock is taken from https://github.com/Unidata/MetPy/blob/main/conftest.py#L169

Now I get this error.

vort = mpcalc.vorticity(u,v,longitude=longitude,latitude=latitude,crs=crs) File "/usr/local/lib/python3.9/dist-packages/metpy/calc/tools.py", line 1096, in wrapper grid_deltas = grid_prototype.metpy.grid_deltas File "/usr/local/lib/python3.9/dist-packages/metpy/xarray.py", line 544, in grid_deltas raise AttributeError( AttributeError: Grid deltas cannot be calculated since horizontal dimension coordinates cannot be found.

[DanielAdriaansen]

First, I would review this section of the documentation:
https://unidata.github.io/MetPy/latest/tutorials/xarray_tutorial.html#coordinates-and-coordinate-reference-systems

This is the relevant part:

If your dataset doesn’t have a CF-conforming grid mapping variable, you can manually specify the CRS using the .assign_crs() method:

So you need to use assign_crs(), and do that using one of the CF-compliant grid mappings listed here:
https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#appendix-grid-mappings

Based on your mock up, it seems like you just have a global latitude/longitude grid is that right? If so, you can do:

u = u.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)
v = v.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)

This assigns a metpy_crs coordinate variable to your Xarray objects which MetPy can interpret for calculations (like vorticity). Depending how you're loading the GFS data, enough projection information is probably present for MetPy to work with it without this step but without seeing the code you're using and how you're accessing the GFS I can't be sure.

But try using assign_crs() and see if it works. And please note, I am assuming the latitude_longitude grid and Earth radius of 6371229 m, so if it is different please choose the grid mapping that you need to use.

[DanielAdriaansen]

Actually, I note that your mock latitude/longitude are 2D, so the grid mapping of latitude_longitude is incorrect here (that assumes 1D lat/lon) so either you need to change to 1D lat/lon if it is a lat/lon grid or choose another grid mapping for your data.

[winash12]

`from metpy.units import check_units
import os
import xarray as xr
import boto3
from botocore import UNSIGNED
from botocore.config import Config

def main():
if (not os.path.isfile('gfs.t12z.pgrb2.0p25.f000')):

client = boto3.client('s3', config=Config(signature_version=UNSIGNED))
client.download_file('noaa-gfs-bdp-pds', 'gfs.20230809/12/atmos/gfs.t12z.pgrb2.0p25.f000', 'gfs.t12z.pgrb2.0p25.f000')

u850 = xr.open_dataset('gfs.t12z.pgrb2.0p25.f000', engine='cfgrib',backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa', 'shortName': 'u', 'level': 850}})
u850 = u850.u

main()`

This is how I invoke the GFS data to load the u850 array. I am not sure where the CRS information is stored in the GFS file.

Regarding your suggestion to add these 2 methods -

u = u.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)
v = v.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)

these should be invoked on the xarray objects after they have been instantiated like this

u = xr.DataArray(u,coords={'longitude':(['longitude','latitude'], longitude), 'latitude':(['longitude','latitude'], latitude)},dims=('longitude','latitude'))

But that will be problematic because the error I reported is thrown when the object is instantiated.

This test case works - https://github.com/Unidata/MetPy/blob/main/tests/calc/test_kinematics.py#L70

and uses 2d arrays. Just not sure what is the CRS they are using.

`

[DanielAdriaansen]

Thanks for sending your sample code.

After this line:
u850 = xr.open_dataset('gfs.t12z.pgrb2.0p25.f000', engine='cfgrib',backend_kwargs={'filter_by_keys':{'typeOfLevel': 'isobaricInhPa', 'shortName': 'u', 'level': 850}})

You can see the projection info in the attributes of the u variable:

>>> u850.u.attrs
{'GRIB_paramId': 131, 'GRIB_dataType': 'fc', 'GRIB_numberOfPoints': 1038240, 'GRIB_typeOfLevel': 'isobaricInhPa',
'GRIB_stepUnits': 1, 'GRIB_stepType': 'instant', 'GRIB_gridType': 'regular_ll', 'GRIB_NV': 0, 'GRIB_Nx': 1440, 'GRIB_Ny': 721,
'GRIB_cfName': 'eastward_wind', 'GRIB_cfVarName': 'u', 'GRIB_gridDefinitionDescription': 'Latitude/longitude. Also called
equidistant cylindrical, or Plate Carree', 'GRIB_iDirectionIncrementInDegrees': 0.25, 'GRIB_iScansNegatively': 0,
'GRIB_jDirectionIncrementInDegrees': 0.25, 'GRIB_jPointsAreConsecutive': 0, 'GRIB_jScansPositively': 0,
'GRIB_latitudeOfFirstGridPointInDegrees': 90.0, 'GRIB_latitudeOfLastGridPointInDegrees': -90.0,
'GRIB_longitudeOfFirstGridPointInDegrees': 0.0, 'GRIB_longitudeOfLastGridPointInDegrees': 359.75, 'GRIB_missingValue':
3.4028234663852886e+38, 'GRIB_name': 'U component of wind', 'GRIB_shortName': 'u', 'GRIB_units': 'm s**-1', 'long_name': 'U
component of wind', 'units': 'm s**-1', 'standard_name': 'eastward_wind'}

MetPy must be able to decipher this automatically, but I'm not totally sure. You can always do this:
u850 = u850.metpy.parse_cf()

Which then attaches the metpy_crs coordinate:

u850
<xarray.Dataset> Size: 4MB
Dimensions:        (latitude: 721, longitude: 1440)
Coordinates:
    time           datetime64[ns] 8B ...
    step           timedelta64[ns] 8B ...
    isobaricInhPa  float64 8B ...
  * latitude       (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude      (longitude) float64 12kB 0.0 0.25 0.5 ... 359.2 359.5 359.8
    valid_time     datetime64[ns] 8B ...
    metpy_crs      object 8B Projection: latitude_longitude
Data variables:
    u              (latitude, longitude) float32 4MB ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP
    history:                 2024-06-05T15:48 GRIB to CDM+CF via cfgrib-0.9.1...

If you're creating data from hand, you have to do all of the work to define the CRS/projection yourself.

But that will be problematic because the error I reported is thrown when the object is instantiated.

I disagree. The error you sent is coming from the line where you are trying to calculate vorticity, not the lines where you are instantiating the Xarray DataArrays:

vort = mpcalc.vorticity(u,v,longitude=longitude,latitude=latitude,crs=crs) File "/usr/local/lib/python3.9/dist
packages/metpy/calc/tools.py", line 1096, in wrapper grid_deltas = grid_prototype.metpy.grid_deltas File
"/usr/local/lib/python3.9/dist-packages/metpy/xarray.py", line 544, in grid_deltas raise AttributeError( AttributeError: Grid deltas
cannot be calculated since horizontal dimension coordinates cannot be found.

What are you trying to accomplish? I have lost track of that. It seems to work for your GFS data. What are you trying to do with your mock data?

[winash12]

I am submitting a PR that uses the vorticity function. Therefore I need to submit a unit test case with mock data that replicates the GFS data. The data can be constant but it needs to be of a sufficient spatial scale so that the voriticity makes physical sense.

` crs = pyproj.CRS(4326)
latitude = np.linspace(-90.,90.,721,endpoint=False)*units.degree
longitude = np.linspace(0.,360. ,1440,endpoint=False)*units.degree
u = np.empty((721,1440))
v = np.empty((721,1440))
longitude,latitude = np.meshgrid(longitude.m_as('degree'),latitude.m_as('degree'))

u = - np.sin(latitude/2)*units('m/s')
v = - np.cos(longitude)*units('m/s')
u = xr.DataArray(u,coords={'longitude':(['longitude','latitude'],
                                        longitude),
                           'latitude':(['longitude','latitude'],
                                       latitude)},dims=('longitude','latitude'))

v = xr.DataArray(v,coords={'longitude':(['longitude','latitude'],longitude),
                           'latitude':(['longitude','latitude'],latitude)},dims=('longitude','latitude'))


u = u.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)
v = v.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)    
vort = mpcalc.vorticity(u,v,longitude=longitude,latitude=latitude,crs=crs)

main()`

Yes you are right with this code I am getting this error

/home/winash12/testpython/test_rot_inv.py:33: UserWarning: Horizontal dimension numbers not found. Defaulting to (..., Y, X) order. vort = mpcalc.vorticity(u,v,longitude=longitude,latitude=latitude,crs=crs) Traceback (most recent call last): File "/home/winash12/testpython/test_rot_inv.py", line 34, in <module> main() File "/home/winash12/testpython/test_rot_inv.py", line 8, in main setUpData() File "/home/winash12/testpython/test_rot_inv.py", line 33, in setUpData vort = mpcalc.vorticity(u,v,longitude=longitude,latitude=latitude,crs=crs) File "/usr/local/lib/python3.9/dist-packages/metpy/calc/tools.py", line 1096, in wrapper grid_deltas = grid_prototype.metpy.grid_deltas File "/usr/local/lib/python3.9/dist-packages/metpy/xarray.py", line 533, in grid_deltas dx, dy = nominal_lat_lon_grid_deltas( File "/usr/local/lib/python3.9/dist-packages/metpy/xarray.py", line 1330, in wrapper result = func(*bound_args.args, **bound_args.kwargs) File "/usr/local/lib/python3.9/dist-packages/metpy/calc/tools.py", line 877, in nominal_lat_lon_grid_deltas raise ValueError( ValueError: Cannot calculate nominal grid spacing from longitude and latitude arguments that are not one dimensional.

[winash12]

@DanielAdriaansen I got it to work eventually - contrived maybe but does work -

`import xarray as xr
import numpy as np
import metpy.calc as mpcalc
from metpy.units import units
import pyproj
import sys
def main():
setUpData()

def setUpData():

lat = np.linspace(90.,-90.,721,endpoint=True)*units.degree
lon = np.linspace(0.,360. ,1440,endpoint=False)*units.degree
u = np.empty((721,1440))
v = np.empty((721,1440))
longitude = xr.DataArray(lon,
                   attrs={'standard_name': 'longitude', 'units': 'degrees_east'})
latitude = xr.DataArray(lat,
                   attrs={'standard_name': 'latitude', 'units': 'degrees_north'})
u1 = -2. * np.sin(lat/2)*units('m/s')
v1 = - np.cos(lon)*units('m/s')
u1,v1 = np.meshgrid(v1,u1)
u = xr.DataArray(u1 ,coords=(latitude,longitude),dims=('latitude','longitude'))

v = xr.DataArray(v1 ,coords=(latitude,longitude),dims=('latitude','longitude'))

u = u.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)
v = v.metpy.assign_crs(grid_mapping_name='latitude_longitude',earth_radius=6371229.0)    
vort = mpcalc.vorticity(u,v)
`