How to transform .grib file into a GeoTIFF with correct projection using Python, GDAL, ArcPy

Question

I am trying to transform a .grib file into a GeoTIFF to be used in a GIS (ArcGIS to be particular), but am having trouble getting the image to project properly. I have been able to create a GeoTIFF, using GDAL in Python, that shows the data but is not showing up in the correct location when brought into ArcGIS. The resulting image is below.

The data I am working with can be downloaded from: https://gimms.gsfc.nasa.gov/SMOS/SMAP/L05/

I am trying to project the data into WGS84 Web Mercator (Auxiliary Sphere), EPSG: 3857

Note: I have tried bringing in the data via ArcMap by creating a Raster Mosaic which should be able to work with .grib data, but I didn't have any luck.

Update: I have also tried using the Project Raster tool, but ArcGIS does not like the default projection that comes from the .grib file and gives an error.

The code I'm using:

import gdal

src_filename = r"C:\att\project\arcshare\public\disaster_response\nrt_products\smap\20150402_20150404_anom1.grib"
dst_filename = r"C:\att\project\arcshare\public\disaster_response\nrt_products\smap\smap_py_test1.tif"

#Open existing dataset
src_ds = gdal.Open(src_filename)

#Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )

#Output to new format
dst_file = driver.CreateCopy( dst_filename, src_ds, 0 )

#Properly close the datasets to flush to disk
dst_ds = None
src_ds = None

I am not very well versed in using GDAL or GDAL in Python, so any help or tips would be greatly appreciated.

I recommend to use cfrgib https://stackoverflow.com/questions/27782574/extract-data-from-grib-weather-model/57338757#57338757. There you will have access to the latitude and longitude coordinates. With `from pyproj import Proj, transform ` you can define two projection and transform the coordinates. — dl.meteo, Aug 04 '19 at 13:10

score 1 · Answer 1 · answered Jul 22 '20 at 17:46

Try using gdal.Translate (in Python) or gdal_translate (from command line). Here are two examples of how I have used each approach in the past:

Option 1: Python approach

from osgeo import gdal

# Open existing dataset
src_ds = gdal.Open(src_filename)

# Ensure number of bands in GeoTiff will be same as in GRIB file. 
bands = [] # Set up array for gdal.Translate(). 
if src_ds is not None:
    bandNum = src_ds.RasterCount # Get band count
for i in range(bandNum+1): # Update array based on band count
    if (i==0): #gdal starts band counts at 1, not 0 like the Python for loop does.
        pass
    else:
        bands.append(i)

# Open output format driver
out_form= "GTiff"

# Output to new format using gdal.Translate. See https://gdal.org/python/ for osgeo.gdal.Translate options.
dst_ds = gdal.Translate(dst_filename, src_ds, format=out_form, bandList=bands)

# Properly close the datasets to flush to disk
dst_ds = None
src_ds = None

Option 2: Command line gdal_translate (called from Python) approach

import os

# Open output format driver, see gdal_translate --formats for list
out_form = "GTiff"

# Pull out specific band of interest
band=3

# Convert from GRIB to GeoTIFF using system gdal_translate
src_ds = src_filename
dst_ds = dst_filename
os.system("gdal_translate -b %s -of %s %s %s" %(str(band), out_form, src_ds, dst_ds))

I've had trouble in the past creating a multi-band GeoTiff using option 2, so I recommend using option 1 when possible.

score -1 · Answer 2 · answered Aug 04 '19 at 13:16

Something like this should transform your native coordinates into your desired projection. This is not tested, yet. (Could by latitude instead of latitudes).

from cfgrib import xarray_store
from pyproj import Proj, transform

grib_data = xarray_store.open_dataset('your_grib_file.grib')
lat = grib_data.latitudes.value
lon = grib_data.longitudes.value
lon_transformed, lat_transformed = transform (Proj(init='init_projection'), 
 Proj(init='target_projection', lon, lat)

How to transform .grib file into a GeoTIFF with correct projection using Python, GDAL, ArcPy

2 Answers2