NWS NDFD Hawaii Grids converting variable(step, values) to variable(step, x, y)

[edrewitz]

Hi everyone!

I know some of you saw me post this question on bluesky social a couple of weeks ago, though I figured I'd post it here too in case anyone who might not have seen it on bluesky might have some ideas.

I am trying to plot the data from the NWS NDFD Hawaii grids on a map in cartopy. The problem is the data is structured as ds['parameter'][step, values] rather than the usual ds['parameter'][step, x, y] (see screenshot). In order for the plot to work I need everything in (step, x, y).

I've tried both using ds.metpy.parse_cf() as well as removing it and trying to assign a custom crs using the ds.metpy.assign_crs() method. Unfortunately, I am still seeing the dataset as (step, values).

I am curious if anyone here has any tips on how to get the Hawaii data into the usual (step, x, y). I've only seen this issue with Hawaii. Alaska and CONUS work fine. Hawaii uses a mercator projection, Alaska a stereographic and CONUS uses Lambert Conformal.

Thanks for all of your help!

[DanielAdriaansen]

@edrewitz I finally had a chance to take a look, and think I have a solution.

First, the documentation here:
https://vlab.noaa.gov/web/mdl/ndfd-spatial-reference-system#guam-hawaii-puerto-rico-grids-northern-pacific-and-oceanic-grids

Gave me some clues that the 2D data were somehow collapsed into 1D (the values dimension). The trick here was actually figuring out how. Based on that documentation, I thought the data were stored left-->right (west-->east) for each row, starting at the bottom row of the grid and moving upward to the top of the grid. That was partially right. The odd thing here, is that every other row has the data reversed, so the data are actually collapsed in a zig-zag pattern such that the first row is the first 321 points of the 72225 length array, but then the second row of the grid is the second set of 321 points (points 322-642), but those points are reversed (i.e., stored right-->left or east-->west). I am not sure if this is intentional or not, but alas I think that's the way it is.

I could not find an elegant way to do this in Xarray directly, so accept my NumPy based solution below:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import sys
import xarray as xr

# Get the NDFD file as the command line argument
ndfd = sys.argv[1]

# Open using cfgrib/xarray
ds = xr.open_dataset(ndfd,engine='cfgrib',backend_kwargs={'indexpath':''})

# Number of rows (y/j dimension)
# "Number of Points on the parallel"
nrow = 225

# Number of columns (x/i dimension)
# Number of Points on the Meridian"
ncol = 321

# Select the first step (there are three forecast times/steps in the file I tested)
ds = ds.isel(step=0)

# Pull out the variable, latitude, and longitude as individual numpy objects
var = ds['tmax'].values
lat = ds['latitude'].values
lon = ds['longitude'].values

# Create some test arrays
var2d = np.empty([nrow,ncol])
lat2d = np.empty([nrow,ncol])
lon2d = np.empty([nrow,ncol])

# Loop over each of the full arrays and insert data from the large vector where it should be in 2D
for i in range(0,nrow):
  start = i*ncol
  end = start+ncol
  # For even rows, the data are stored west-->east
  if i%2==0:
    var2d[i,:] = var[start:end]
  # For odd rows, the data are stored east-->west
  else:
    var2d[i,:] = np.flip(var[start:end],axis=0)

  # Lat/lon data are the same for each grid cell in this projection
  lat2d[i,:] = lat[start:end]
  lon2d[i,:] = lon[start:end]

# For the 1D lat/lon coordinates, it's the first row of longitudes and
# the first column of latitudes
lon1d = lon2d[0,:]
lat1d = lat2d[:,0]

# Convert to Celsius for convenience
outvar = var2d-273.15

# Quick plot with Cartopy directly
fig = plt.figure(figsize=(22,15))
ax = plt.axes(projection=ccrs.Mercator(central_longitude=202,min_latitude=10,max_latitude=25.0))
ax.pcolormesh(lon1d,lat1d,outvar,transform=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5)
ax.add_feature(cfeature.STATES, linewidth=0.5)
ax.add_feature(cfeature.BORDERS, linewidth=0.5)
plt.show()

Give this a shot and see if you get something reasonable. I've attached a screengrab of what I get when I run this code.