CODEWITHSUNDEEP
pythontransformgdalgeotiffgrib
glayne
glayne

Reputation: 175

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

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.

What the data looks like when I put it into ArcGIS

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.

Upvotes: 1

Views: 5253

Answers (2)

WxCZar
WxCZar

Reputation: 31

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.

Upvotes: 1

dl.meteo
dl.meteo

Reputation: 1786

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)

Upvotes: -1

Related Questions