Why does the HRRR model have 2-dimensional grids for latitude and longitude?

[WToma]

After reading https://cfconventions.org/Data/cf-conventions/cf-conventions-1.1/build/ch05s02.html and several Stackoverflow answers I'm still unclear how to properly interpret the situation when the latitude and longitude coordinates are 2D themselves.

Following the example on the Herbie main page:

>>> H = Herbie('2021-01-01 12:00', model='hrrr', product='sfc', fxx=6)
>>> H.download(':500 mb')
>>> x = H.xarray('TMP:2 m')
>>> x.coords['latitude']
<xarray.DataArray 'latitude' (y: 1059, x: 1799)>
array([[21.138123, 21.14511 , 21.15209 , ..., 21.154509, 21.147531, 21.140547],
       [21.162995, 21.169984, 21.176967, ..., 21.179387, 21.172407, 21.165419],
       [21.187869, 21.194861, 21.201847, ..., 21.204268, 21.197285, 21.190294],
       ...,
       [47.789559, 47.799849, 47.810129, ..., 47.813691, 47.803415, 47.793128],
       [47.814093, 47.824386, 47.834669, ..., 47.838233, 47.827953, 47.817664],
       [47.838623, 47.84892 , 47.859206, ..., 47.862771, 47.852488, 47.842195]])
Coordinates:
    time               datetime64[ns] 2021-01-01T12:00:00
    step               timedelta64[ns] 06:00:00
    heightAboveGround  float64 2.0
    latitude           (y, x) float64 21.14 21.15 21.15 ... 47.86 47.85 47.84
    longitude          (y, x) float64 ...
    valid_time         datetime64[ns] 2021-01-01T18:00:00
Dimensions without coordinates: y, x
Attributes:
    units:          degrees_north
    standard_name:  latitude
    long_name:      latitude
# similar for longitude, but the values are different

It looks like both latitude and longitude describe the whole grid each just by themselves, and two grids are different.

Trying to make sense of this, I looked at pluck_points (it's easier for me to understand than nearest_points):

Herbie/hrrrb/tools.py

Lines 132 to 143 in f1ae91b

	# Find absolute difference between requested point and the grid coordinates.
	abslat = np.abs(ds.latitude - p_lat)
	abslon = np.abs(ds.longitude - p_lon)
	# Create grid of the maximum values of the two absolute grids
	c = np.maximum(abslon, abslat)
	# Find location where lat/lon minimum absolute value intersects
	if ds.latitude.dims == ('y', 'x'):
	y, x = np.where(c == np.min(c))
	elif ds.latitude.dims == ('x', 'y'):
	x, y = np.where(c == np.min(c))

The way I'm reading this is: "the location you're searching for may be best matched either by a point in latitide, or by a point in longitude, or by point wherein the lat is in latitude and the long is in longitude, or the other way". In other words, latitude and longitude are both full grids in and of themselves, we might as well just call them grid1 and grid2. If we entirely disregard one of them, we still have complete information. I'm extremely suspicious of this interpretation, but it's the only one that makes any sense in my head so far :D

So with that, the questions are

• what's the best way to think of these 2D latitude / longitude coordinates?
• is there a rule, or rule of thumb, about when should I use one or the other?
• ... or are they OK to mix and match, like in my interpretation above?

Thank you!

[blaylockbk]

Many limited-area NWP models (the HRRR included) do not have a Cartesian coordinate system. The HRRR model domain is on a lambert-conformal projection, and 2D coordinates for both lat and lon are required to describe the location of each model point on the earth. The alternative is to use the projection's coordinate system, but that is not as intuitive as lat/lon coordinates.

1. The best way to think of 2D lat/lon is to understand that they are on projected coordinates.
2. If you want to locate the data on a map, you need to know both 2D lat/lon coordinates.

For what it's worth, plucking out nearest-neighbor points for curvilinear grids is not trivial; xarray has been trying to add this capability for several years. My pluck_points method you were looking at is illustrated here: https://github.com/blaylockbk/pyBKB_v3/blob/master/demo/Nearest_lat-lon_Grid.ipynb. That might help you see how it finds the nearest point in a grid.

[blaylockbk]

That is exactly right!

Here is how the HRRR 2-m temperatures looks on an equirectangular (PlateCarree) map versus a Lambert-conformal map; same data, different map projections. Notice how the parallels curve in the second plot.

I highlighted rows 600-650 in grey (ds.isel(y=slice(600-650))

[blaylockbk]

Thanks for the good question.

[williamhobbs]

@blaylockbk, this might be a basic (and slightly off-topic) question, but how did you get the latitude and longitude values to show up on the plot axes in the examples above?

It seems that something like

ax.set_xticks(ds.longitude)
ax.set_yticks(ds.latitude)

would work, but latitude and longitude are 2D with HRRR.

[blaylockbk]

Does this answer your question? https://scitools.org.uk/cartopy/docs/latest/matplotlib/gridliner.html

[williamhobbs]

Yes! Thanks!

[maxdow]

Thanks for this helpful informations :)
Is there a way to convert this HRRR grid to a regular cartesian grid ? I can't find a way to do this but i'm sure it is possible with xarray ?

I try to achieve this kind of indexing ( example from Icon dataset )

[blaylockbk]

As far as I understand, the only way to convert the lat/lon to a cartesian grid is to transform the lat/lon points to the HRRR map projection coordinates (lambert conic conformal). You can look into cartopy and metpy for methods to do that.

I've done something like this before in Herbie here

Herbie/herbie/accessors.py

Lines 262 to 276 in a45da72

	# Apply the MetPy method `assign_y_x` to the dataset
	# https://unidata.github.io/MetPy/latest/api/generated/metpy.xarray.html?highlight=assign_y_x#metpy.xarray.MetPyDataArrayAccessor.assign_y_x
	ds = ds.metpy.assign_y_x()
	# Convert the requested [(lon,lat), (lon,lat)] points to map projection.
	# Accept a list of point tuples, or Shapely Points object.
	# We want to index the dataset at a single point.
	# We can do this by transforming a lat/lon point to the grid location
	crs = ds.metpy_crs.item().to_cartopy()
	transformed_data = crs.transform_points(
	ccrs.PlateCarree(), point_df.longitude, point_df.latitude
	)
	xs = transformed_data[:, 0]
	ys = transformed_data[:, 1]

[maxdow]

thanks for your feedback.
I tested this approach with pyproj

import cartopy.crs as ccrs

projection = ccrs.LambertConformal(
    central_longitude=262.5,
    central_latitude=38.5,
    standard_parallels=(38.5, 38.5),
    globe=ccrs.Globe(semimajor_axis=6371229, semiminor_axis=6371229),
)
transformer = pyproj.Transformer.from_crs(projection.to_wkt(),4326,always_xy=True)

xv, yv = np.meshgrid(ds.x, ds.y)
lon, lat = transformer.transform(xv, yv)

I will try to figure out how to change the dimensions on the grid according to this projection