Reproject file to line up with geographic features on a map

I have been trying to reproject data to match the "geographic features" on a map (and also shapefiles) in QGIS or on the web. I’m ingesting raw GRIB data (available here). I’m using python to open the grib and get the values. I’m the using gdal to try to save to file.

Here is the python script:

import numpy as np
import pygrib
from osgeo import osr, gdal

grbs = pygrib.open('data/hrrr.t23z.wrfprsf00.grib2')
lats, lons = grbs[607].latlons()
surfacePres = grbs[607].values

flipped_data = np.flipud(surfacePres)

nRows = surfacePres.shape[0]
nCols = surfacePres.shape[1]

xres = (lats.max() - lats.min()) / nRows
yres = (lons.max() - lons.min()) / nCols

driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create('sfcpres.tif', nCols, nRows, 1, gdal.GDT_Float32)
dst_ds.SetGeoTransform([lons.min(), yres, 0, lats.max(), 0, -xres])

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS("WGS84")
dst_ds.SetProjection(srs.ExportToWkt())

band = dst_ds.GetRasterBand(1)
band.SetNoDataValue(9999)
band.WriteArray(flipped_data)

I have printed nRows, nCols, xres, and yres.

nRows 1059
nCols 1799
xres 0.029723824652320985
yres 0.04067720231374391

When I go to visualize my tif on a map (opening in QGIS), the file is not matching up with the base maps. You can clearly see this near the great lakes, and if you look closely on the west coast, the "darker colors" (purple, etc) are over water, when they should actually be over land.

When I go to open the original grib2 file (attached above – the "file"), QGIS uses proj4 in order to overlay it on my map. This depiction is what I am expecting.

Finally, this last image shows the difference between the two. The gray image is what I am currently generating, and the colored image is what I am trying to get.

from osgeo import osr, gdal

dataset = gdal.Open('data/hrrr.t23z.wrfprsf00.grib2')

band607 = gdal.Translate('band607.tif', dataset, bandList=[607])

warp_srs = osr.SpatialReference()
warp_srs.ImportFromEPSG(4326)

warped_b607 = gdal.Warp('warped_b607.tif', band607, dstSRS=warp_srs, dstNodata=9999)

dataset = None
band607 = None
warped_b607 = None