Spatial reference transformation for GRIB files using Python

Question

I want to use the Python programming language to create a transformation of the spatial reference system used by NOAA in the GRIB files created by their WaveWatchIII model (link to grib files).

Using the following code to import the files:

import gdal    

file = 'example.grib'
raw_data = gdal.Open(file, gdal.GA_ReadOnly)
message_count = raw_data.RasterCount
print(message_count) # files are onedimensional
print(raw_data.RasterXSize, raw_data.RasterYSize)#360 181    
message = raw_data.GetRasterBand(1) # single banded files

So in each grib file there is only one layer. I have a list of GPS coordinates for which I want to look up the corresponding data value from the grib file. I tried to transform the grib file to a different projection with no success. I have visualized the grib file using the following code:

a = message.ReadAsArray()
plt.figure()
plt.imshow(a, cmap='hot', interpolation='nearest', vmin=0, vmax=10)

This yields the following picture:

The following code is an attempt to change the projection:

import osr
# get projection from grib file
source = raw_data.GetSpatialRef()

# gps coordinate system
target = osr.SpatialReference()
target.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(source, target)
transform.TransformPoint(34, 59)

The GetSpatialRef() function returns the following:

'GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST]]

For the list of GPS coordinates I have I want to match it with the closest data point in the grib grid. The transformation in the code above simply returns (59.0, 34.0, 0.0)

Answer 1

You can use rioxarray:

Load in the data and convert the longitude coordinates

import rioxarray
rds = rioxarray.open_rasterio(
    "US058GOCN-GR1mdl.0110_0240_00000F0RL2020041000_0001_000000-000000scdy_wav_per"
)
# convert from 0-360 to (-180 to 180) to match EPSG:4326
rds = rds.assign_coords(x=(((rds.x + 180) % 360) - 180)).sortby('x')

The point is the same in this projection as it is in EPSG:4326

from pyproj import Transformer

transformer = Transformer.from_crs(rds.rio.crs, "EPSG:4326", always_xy=True)
transformer.transform(34, 59)

(34.0, 59.0)

You can pull out the value from the grid like so:

rds.sel(x=34, y=59, method='nearest').values

array([9999.])

Answer 2

I took the solution for warping the image from Frank Warmerdam (How to reproject raster from 0 360 to -180 180 with cutting 180 meridian) with the gdal commandline tool, to transform the raster into WGS84:

gdalwarp -t_srs WGS84 ~/0_360.tif 180.tif  -wo SOURCE_EXTRA=1000 \
           --config CENTER_LONG 0

This then allows be, to convert geographic coordinates into the row- and columns index of a numpy array to get the pixel value at a points position:

import gdal

file = r'c:\tmp\180.tif'
ds = gdal.Open(file, gdal.GA_ReadOnly)
count = ds.RasterCount
print(count) # files are onedimensional
print(ds.RasterXSize, ds.RasterYSize)#360 181
rb = ds.GetRasterBand(1) # single banded files

a = rb.ReadAsArray()

def from_pixel(px, py, rx, ry, dx,dy):

    x = dx*rx + px + rx/2
    y = dy*ry + py + ry/2
    return x,y

def from_point(px, py, rx, ry, x, y):

    # x = dx*rx + px + rx/2
    # y = dy*ry + py + ry/2
    dx = (x - px - rx/2) / rx
    dy = (y - py - ry/2) / ry
    return int(dx), int(dy)


px = ds.GetGeoTransform()[0]
py = ds.GetGeoTransform()[3]
rx = ds.GetGeoTransform()[1]
ry = ds.GetGeoTransform()[5]

print(px, py, rx, ry)
dx, dy = from_point(px,py,rx,ry,-43,45)
print(dx,dy,a[dy][dx])

The function from_pixel() is made after this post: Finding the pixel coordinates using GDAL/numpy, but I needed the invers function from_point to read the raster value. I choose the point (-43,54) lon/lat, which gives me raster index and value: 136, 45, 9.0297. I successfully checked the position/value in ArcGIS (QGIS will do as well) to be correct.