How to convert X/Y model data to lat/long in MetPy?

[rolypolyman]

I am having some difficulty getting THREDDS HRRR data to import because it appears to come in X,Y coordinates instead of lat/long. I am wondering how I am supposed to transform it to lat/long so it can be plotted on a map. If there is a code sample somewhere I can look at that.

Thanks.

Pseudocode:

query.variables('Pressure_reduced_to_MSL_msl')
data = xr.open_dataset(NetCDF4DataStore(data))
slp = data['Pressure_reduced_to_MSL_msl']  # HRRR
slp = slp.sel(height_above_ground3 = 2 * units.meters, method='nearest').isel(time=time_index)
lon = **?????**
lat = **?????**
lon_2d, lat_2d = np.meshgrid(lon, lat)
crs = ccrs.LambertConformal(central_longitude=central_long, central_latitude=45)
data_crs = ccrs.PlateCarree()

<xarray.DataArray 'Pressure_reduced_to_MSL_msl' (time: 1, y: 555, x: 742)>
[411810 values with dtype=float32]
Coordinates:
    reftime  (time) datetime64[ns] 2024-03-18
  * x        (x) float32 -1.387e+03 -1.384e+03 -1.382e+03 ... 490.1 492.7 495.2
  * y        (y) float32 231.5 234.0 236.5 ... 1.633e+03 1.636e+03 1.638e+03
  * time     (time) datetime64[ns] 2024-03-18T06:00:00
Attributes:
    long_name:                       Pressure reduced to MSL @ Mean sea level
    units:                           Pa
    Grib_Variable_Id:                VAR_7-0--1-1_L101
    Grib2_Parameter:                 [0 3 1]
    Grib2_Parameter_Discipline:      Meteorological products
    Grib2_Parameter_Category:        Mass
    Grib2_Parameter_Name:            Pressure reduced to MSL
    Grib2_Level_Type:                101
    Grib2_Level_Desc:                Mean sea level
    Grib2_Generating_Process_Type:   Forecast
    Grib2_Statistical_Process_Type:  UnknownStatType--1
    grid_mapping:                    LambertConformal_1377X2145-38p23N-95p44W-2

[dopplershift]

What kind of plot are you trying to do? In my experience, it's quite rare that you should need lat/lon to do plotting, but you should be able to plot using the x/y so long as you pass a proper projection (which MetPy can generate from the metadata).

[rolypolyman]

Just a simple cartopy map of sea level pressure. All the examples I've viewed so far (MetPy Mondays etc) work with lat/long, and the map is in lat/long, but here THREDDS is giving X/Y so I'm not really understanding how to get from A to B.

[dopplershift]

Here's an example of how to plot without needing to do any manual coordinate manipulation:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy
import numpy as np
from siphon.catalog import TDSCatalog
import xarray as xr
from xarray.backends import NetCDF4DataStore

# Find the latest HRRR run
latest_hrrr = TDSCatalog('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/HRRR/CONUS_2p5km/catalog.xml').latest

# Download a subset of the data and open with xarray
ncss = latest_hrrr.subset()
query = ncss.query()
query.lonlat_box(-120, -60, 25, 50)
query.variables('Pressure_reduced_to_MSL_msl')
nc = ncss.get_data(query)
ds = xr.open_dataset(NetCDF4DataStore(nc))

# Pull out the pressure field at the first time, parsing the
# CF metadata for projection information
mslp = ds.metpy.parse_cf('Pressure_reduced_to_MSL_msl')[0]

# Convert pressure to mbar/hPa for easy contouring
mslp = mslp.metpy.convert_units('mbar')

# Set up a plot on a Azimuthal Equal Area/Plate Carree map
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Make use of metadata to find the x/y coordinates and the native data projection
x = mslp.metpy.x
y = mslp.metpy.y
data_proj = mslp.metpy.cartopy_crs

# Plot contours of every 4 mbar, passing transform so Cartopy knows how to properly
# georeference the x/y coordinates.
ax.contour(x, y, mslp, levels=np.arange(900, 1100, 4), transform=data_proj)

# Add coastlines and state borders
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.STATES)

which gives:

[rolypolyman]

Thank you, perfect, now it makes a little more sense. I appreciate the help.