Unidata / MetPy

MetPy is a collection of tools in Python for reading, visualizing and performing calculations with weather data.
https://unidata.github.io/MetPy/
BSD 3-Clause "New" or "Revised" License
1.25k stars 414 forks source link

xarray accessor assign_y_x and HRRR grib data #2940

Open keltonhalbert opened 1 year ago

keltonhalbert commented 1 year ago

What went wrong?

I spent some time looking through open issues and stack overflow to make sure I am doing this correctly, but it appears MetPy 1.4 with Python 3.10 is broken with the assign_y_x accessor. Previously this worked in adding the X and Y coordinates in meters to HRRR grib data, but running the example provided no longer works. Downgrading to MetPy 1.3 didn't fix the issue either, but it could be something tied up in upgrading from Python 3.8 to Python 3.10.

I also tried setting the tolerance argument (100 meters, 1000 meters, 3000 meters, and 6000 meters), and still am not having success. Any suggestions or ideas to recover this functionality?

Operating System

Linux

Version

1.4

Python Version

3.10

Code to Reproduce

import xarray as xr
import fsspec
import metpy
import pint

def main():

    print(xr.__version__)
    print(metpy.__version__)
    url = "simplecache::gs://high-resolution-rapid-refresh/hrrr.20230303/conus/hrrr.t00z.wrfsfcf00.grib2"
    ds = xr.open_dataset(fsspec.open_local(url, simplecache={'cache-storage': '/tmp'}), engine="cfgrib",
                         backend_kwargs={'filter_by_keys': {'typeOfLevel': 'surface', 'shortName': 'cape' }})
    ds = ds.metpy.assign_crs(
        grid_mapping_name = "lambert_conformal_conic",
        latitude_of_projection_origin=ds.cape.GRIB_LaDInDegrees,
        longitude_of_projection_origin=ds.cape.GRIB_LoVInDegrees,
        standard_parallel=(ds.cape.GRIB_Latin1InDegrees, ds.cape.GRIB_Latin2InDegrees),
        earth_radius=6371229.)
    print(ds.latitude.values)
    print(ds.longitude.values)
    print(ds)

    ds = ds.metpy.assign_y_x(tolerance=pint.Quantity(6000.0, "m"))
    print(ds)

if __name__ == "__main__":
    main()

Errors, Traceback, and Logs

2022.12.0
1.3.0
[[21.138123   21.14511004 21.1520901  ... 21.1545089  21.14753125
  21.14054663]
 [21.16299459 21.1699845  21.17696744 ... 21.17938723 21.1724067
  21.16541921]
 [21.18786863 21.19486142 21.20184723 ... 21.20426802 21.19728462
  21.19029425]
 ...
 [47.78955926 47.799849   47.81012868 ... 47.81369093 47.80341474
  47.79312849]
 [47.81409316 47.82438621 47.8346692  ... 47.83823259 47.8279531
  47.81766354]
 [47.8386235  47.84891986 47.85920615 ... 47.86277069 47.85248789
  47.84219502]]
[[237.280472   237.30713868 237.3338097  ... 287.6569408  287.68361332
  287.71028151]
 [237.27297501 237.29964881 237.32632695 ... 287.66442108 287.69110073
  287.71777603]
 [237.26547368 237.2921546  237.31883986 ... 287.6719057  287.69859247
  287.7252749 ]
 ...
 [225.9351904  225.97171577 226.00825329 ... 298.97907406 299.0156158
  299.05214538]
 [225.91986142 225.95639874 225.99294822 ... 298.99437498 299.03092868
  299.06747022]
 [225.90452027 225.94106954 225.97763099 ... 299.00968806 299.04625373
  299.08280723]]
<xarray.Dataset>
Dimensions:     (y: 1059, x: 1799)
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
    surface     float64 ...
    latitude    (y, x) float64 21.14 21.15 21.15 21.16 ... 47.86 47.85 47.84
    longitude   (y, x) float64 237.3 237.3 237.3 237.4 ... 299.0 299.0 299.1
    valid_time  datetime64[ns] ...
    metpy_crs   object Projection: lambert_conformal_conic
Dimensions without coordinates: y, x
Data variables:
    cape        (y, x) float32 ...
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:                 2023-03-03T22:33 GRIB to CDM+CF via cfgrib-0.9.1...
/home/wxbyte/installs/miniconda3/envs/zarr-tiler/lib/python3.10/site-packages/metpy/xarray.py:355: UserWarning: More than one time coordinate present for variable "cape".
  warnings.warn('More than one ' + axis + ' coordinate present for variable'
Traceback (most recent call last):
  File "/home/wxbyte/installs/zarr-tiler/demo_bug.py", line 28, in <module>
    main()
  File "/home/wxbyte/installs/zarr-tiler/demo_bug.py", line 23, in main
    ds = ds.metpy.assign_y_x(tolerance=pint.Quantity(6000.0, "m"))
  File "/home/wxbyte/installs/miniconda3/envs/zarr-tiler/lib/python3.10/site-packages/metpy/xarray.py", line 948, in assign_y_x
    y, x = _build_y_x(grid_prototype, tolerance)
  File "/home/wxbyte/installs/miniconda3/envs/zarr-tiler/lib/python3.10/site-packages/metpy/xarray.py", line 1147, in _build_y_x
    raise ValueError('Projected y and x coordinates cannot be collapsed to 1D within '
ValueError: Projected y and x coordinates cannot be collapsed to 1D within tolerance. Verify that your latitude and longitude coordinates correspond to your CRS coordinate.
keltonhalbert commented 1 year ago

It appears that with some modification, I can force it to work. I first have a statically defined proj string that I assign, and then have to manually specify force and tolerance values. If this is anticipated or expected behavior, feel free to close the issue. However, there is a lot of example code out there that uses the previous approach, so maybe still worth fixing?

import xarray as xr
from pyproj import CRS
import fsspec
import metpy
import pint

HRRR_PROJ_STR = "+proj=lcc +a=6371229 +b=6371229 +lon_0=262.5 +lat_0=38.5 +lat_1=38.5 +lat_2=38.5 +type=crs"

def main():

    print(xr.__version__)
    print(metpy.__version__)
    url = "simplecache::gs://high-resolution-rapid-refresh/hrrr.20230303/conus/hrrr.t00z.wrfsfcf00.grib2"
    ds = xr.open_dataset(fsspec.open_local(url, simplecache={'cache-storage': '/tmp'}), engine="cfgrib",
                         backend_kwargs={'filter_by_keys': {'typeOfLevel': 'surface', 'shortName': 'cape' }})
    print(ds)
    ds = ds.metpy.assign_crs(CRS(HRRR_PROJ_STR).to_cf())
    ds = ds.metpy.parse_cf().metpy.assign_y_x(force=False, tolerance=None)
    print(ds)

if __name__ == "__main__":
    main()
jthielen commented 1 year ago

I think this may be a typo/copy-paste error! Note that in your initial example code, you used an invalid key longitude_of_projection_origin rather than longitude_of_central_meridian (as used in the CF conventions for LCC). If you update this, your original code should run as-is, even with the default tolerance:

import xarray as xr
import fsspec
import metpy

url = "simplecache::gs://high-resolution-rapid-refresh/hrrr.20230303/conus/hrrr.t00z.wrfsfcf00.grib2"
ds = xr.open_dataset(fsspec.open_local(url, simplecache={'cache-storage': '/tmp'}), engine="cfgrib",
                     backend_kwargs={'filter_by_keys': {'typeOfLevel': 'surface', 'shortName': 'cape' }})
ds = ds.metpy.assign_crs(
    grid_mapping_name = "lambert_conformal_conic",
    latitude_of_projection_origin=ds.cape.GRIB_LaDInDegrees,
    longitude_of_central_meridian=ds.cape.GRIB_LoVInDegrees,
    standard_parallel=(ds.cape.GRIB_Latin1InDegrees, ds.cape.GRIB_Latin2InDegrees),
    earth_radius=6371229.)
print(ds.latitude.values)
print(ds.longitude.values)
print(ds)

ds = ds.metpy.assign_y_x()
print(ds)
[[21.138123   21.14511004 21.1520901  ... 21.1545089  21.14753125
  21.14054663]
 [21.16299459 21.1699845  21.17696744 ... 21.17938723 21.1724067
  21.16541921]
 [21.18786863 21.19486142 21.20184723 ... 21.20426802 21.19728462
  21.19029425]
 ...
 [47.78955926 47.799849   47.81012868 ... 47.81369093 47.80341474
  47.79312849]
 [47.81409316 47.82438621 47.8346692  ... 47.83823259 47.8279531
  47.81766354]
 [47.8386235  47.84891986 47.85920615 ... 47.86277069 47.85248789
  47.84219502]]
[[237.280472   237.30713868 237.3338097  ... 287.6569408  287.68361332
  287.71028151]
 [237.27297501 237.29964881 237.32632695 ... 287.66442108 287.69110073
  287.71777603]
 [237.26547368 237.2921546  237.31883986 ... 287.6719057  287.69859247
  287.7252749 ]
 ...
 [225.9351904  225.97171577 226.00825329 ... 298.97907406 299.0156158
  299.05214538]
 [225.91986142 225.95639874 225.99294822 ... 298.99437498 299.03092868
  299.06747022]
 [225.90452027 225.94106954 225.97763099 ... 299.00968806 299.04625373
  299.08280723]]
<xarray.Dataset>
Dimensions:     (y: 1059, x: 1799)
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
    surface     float64 ...
    latitude    (y, x) float64 21.14 21.15 21.15 21.16 ... 47.86 47.85 47.84
    longitude   (y, x) float64 237.3 237.3 237.3 237.4 ... 299.0 299.0 299.1
    valid_time  datetime64[ns] ...
    metpy_crs   object Projection: lambert_conformal_conic
Dimensions without coordinates: y, x
Data variables:
    cape        (y, x) float32 ...
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:                 2023-03-03T22:49 GRIB to CDM+CF via cfgrib-0.9.1...
<xarray.Dataset>
Dimensions:     (y: 1059, x: 1799)
Coordinates:
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
    surface     float64 ...
    latitude    (y, x) float64 21.14 21.15 21.15 21.16 ... 47.86 47.85 47.84
    longitude   (y, x) float64 237.3 237.3 237.3 237.4 ... 299.0 299.0 299.1
    valid_time  datetime64[ns] ...
    metpy_crs   object Projection: lambert_conformal_conic
  * y           (y) float64 -1.587e+06 -1.584e+06 ... 1.584e+06 1.587e+06
  * x           (x) float64 -2.698e+06 -2.695e+06 ... 2.693e+06 2.696e+06
Data variables:
    cape        (y, x) float32 ...
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:                 2023-03-03T22:49 GRIB to CDM+CF via cfgrib-0.9.1...

I am not sure what would have changed in regards to this error after updating to more recent versions, however. Perhaps PROJ/pyproj became more strictly conforming to CF conventions?