TransWikia.com

Reproject file to line up with geographic features on a map

Geographic Information Systems Asked by Craytor on July 5, 2021

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.
enter image description here

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.
enter image description here

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.
enter image description here

One Answer

Source grib file is projected in Lambert Conformal Conic projection, with 3000x3000 meters pixels. You are trying to copy the same values in the same rows and columns but in geographic coordinates. That method doesn't warp the information. If you look at the image, it still looks like a LCC projected image.

You can extract the 607 band and warp it:

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

Correct answer by Gabriel De Luca on July 5, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP