python plot of wrf output does not align with cartopy feature lines

Question

I am trying to plot precipitation from the CONUS404 project that was produced with wrf. The map projection of the wrf output is Lambert Conformal. However, when I make the plot in cartopy and use the state boundaries feature, the data does not appear in the correct locations. For example, in the plot the dark blue band is supposed to line up with St. Louis, MO, but is plotted to the north.

My code:

import numpy as np
import xarray as xr
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import geocat.viz.util as gv
import cmaps

### Import and process data 
conus404 = xr.open_dataset('wrf2d_d01_2022-07-26.nc')
var = conus404['PREC_ACC_NC'].sum(dim='Time').to_numpy()
xlat, xlong = conus404['XLAT'].to_numpy(), conus404['XLONG'].to_numpy()

### Figure parameters
mpl.colormaps.unregister('cmp_haxby_r')
colormap = cmaps.cmp_haxby_r
LatMax,LatMin,LonMax,LonMin = 57, 18, -125, -70
data_proj = ccrs.PlateCarree()
map_proj = ccrs.LambertConformal(central_longitude=conus404.STAND_LON,central_latitude=conus404.MOAD_CEN_LAT,
                                 standard_parallels=(conus404.TRUELAT1, conus404.TRUELAT2))

### Intialize figure
fig, ax = plt.subplots(figsize=(8,8),dpi=300,subplot_kw=dict(projection=map_proj))
ax.set_extent([LonMin,LonMax,LatMin,LatMax],data_proj)
ax.add_feature(cfeature.COASTLINE,linewidths=0.25)
ax.add_feature(cfeature.STATES,linewidths=0.15)

### Ticks and gridlines
gl = ax.gridlines(crs=data_proj,draw_labels=True,alpha=0.5,lw=0.1,color='k',linestyle='--',
                  x_inline=False,y_inline=False)
gl.top_labels, gl.bottom_labels, gl.left_labels, gl.right_labels = False, True, True, False
gl.ylocator, gl.xlocator = mpl.ticker.MultipleLocator(10.), mpl.ticker.MultipleLocator(10.)
gl.ylabel_style = {'size':10,'rotation':0,'ha':'right','va':'center'}
gl.xlabel_style = {'size':10,'rotation':0,'ha':'center','va':'top'}
gl.ypadding, gl.xpadding = 5, 5

### Plot variable
method = 'pcolormesh'
if method == 'pcolormesh': # Much faster but doesn't allow user to specify number of color intervals
    cntr = ax.pcolormesh(xlong,xlat,var,cmap=colormap,transform=data_proj,vmin=0,vmax=60) 
if method == 'contourf': # Much slower but allows user to specify number of color intervals
    cntr = ax.contourf(xlong,xlat,var,cmap=colormap,transform=data_proj,levels=np.arange(0,61,1),extend='both')

### Color Bar
bar = fig.colorbar(cntr,ax=ax,orientation='vertical',ticks=mpl.ticker.MultipleLocator(20.),
                   shrink=0.6,pad=0.02,spacing='proportional')
bar.ax.yaxis.set_minor_locator(mpl.ticker.MultipleLocator(1))
bar.ax.tick_params(labelsize=10)
bar.set_label('mm',rotation=0,loc='center',labelpad=0.1)

The weird thing is that I do not have this issue when I run the same plotting code for a different wrf variable (pressure: 'P') as the data plots exactly where it should.

Any idea why precipitation seems to be out of place?