GDAL – Troubleshooting GeoTIFF Creation Using gdal_translate

gdalgdal-translategdal2tiles

I am converting weather radar data (grib2) to raster tiles for Mapbox. The resulting tiles are mis-aligned from the original grib2 data and I am not sure why. I have a 4-step process using GDAL in the command line to convert the grib2 to raster tiles, which I will go into great detail below. The weather data provider provides an online viewer for their data, so I was able to make an example by matching the zoom and capturing the difference :

enter image description here

I am providing the grib2 file for anyone that needs it on Google Drive.

Edit : This is the colorization file needed for testing

Here is my process to convert this data to tiles, and I suspect the problem is with gdalwarp somehow, but I am not sure. EDIT – Based on comments, It now seems gdal_translate is the likely culprit.

  1. Use gdal_translate to convert the grib2 file to a GeoTiff

gdal_translate -b 1 MRMS_ReflectivityAtLowestAltitude_00.50_20211212-172838.grib2 Output.tif

  1. Use gdalwarp to reproject the Output.tif to 3857 for Mapbox.

gdalwarp -t_srs EPSG:3857 Output.tif Output-projected.tif

  1. I use gdaldem to colorize the file (color is different, this is intentional).

gdaldem color-relief -alpha Output-projected.tif colorization_file.txt Output-colorized.tif

  1. Last, I use gdal2tiles to tile the image for Mapbox.

python gdal2tiles.py --profile=mercator --processes=8 -x --zoom=0-13 Output-colorized.tif outputDir

I have ran gdalinfo on the output from each step in the process so that anyone may see in detail :

As a bonus, if anyone wanted to visually check this exact data location it is located in Florida at -81.42, 28.85. This is in case someone wants to check the data using their own software.

Edit : For clarity, I am using GDAL 3.4.0, released 2021/11/04

With all of the information above, I can only guess at what I am doing wrong. After hours of fiddling with the script I am still pretty clueless. I have an idea that it has to do with gdalwarp because this is the initial point where I had attempted to change the projection. If anyone needs more information, I will gladly edit this post to include it.

Best Answer

I downloaded the grib2 file to do the process step-by-step and identify in which point does this error occur. In a first glance, I also thought it might be an issue in the first step (transforming grib2 to tif), but I then found that the apparent shift happens in the GDAL Warp step and is simply due to a nearest neighbor interpolation.

First step:

Such as the author, I also used the snippet available in my website for saving a grib2 file to into Python.

(I would appreciate if the question author could edit the question and add that he has been trying to use this script for step #1 as well, thanks!)

To be exact, this is the script:

from osgeo import gdal
arq1 = gdal.Open (‘gribfile.grib2’)
GT_entrada = arq1.GetGeoTransform()
print(GT_entrada)
save_arq = gdal.Translate('tif_file.tif’,arq1)

I then opened both the grib2 and the tif files in QGIS. Same pixel size, same projection, same extent, same everything. I do not think this is where the shift happens.

Second step:

I didn’t do the GDAL Warp processing on CLI, I chose to use the Warp (Reproject) tool on QGIS, which uses the same algorithms. If you fill only the fields of the input raster and the output projection, it does the same calculation as shown in the question. The default interpolation type is Nearest Neighbour (or Neighbor), which will “copy” the value of the closest pixel on the original raster to the pixel on the new one. Because other values, such as the pixel size, and the extent, were not provided to the algorithm, the pixel size changed in this operation, as well as the number of pixels in each direction. This change off-centered the locations of pixels, in respect to the original raster. Combining this change with the interpolation used, it may appear, in first sight, that the pixels have “shifted”. I have another example of this “occurrence” on my blog. If I force the extent to change and the pixel size to stay the same, while performing Nearest Neighbor, I see a shift in the pixels. To make sure it was that, I also did the same operation (GDAL Warp), but with bilinear interpolation. In the resulting raster, an outer region filled with the interpolated values became visible. This reinforces the other findings.

I did not go further into the steps because I genuinely think that solving this will solve the whole issue.

I could not think of a complete solution, but I have some ideas of workarounds (but I am not sure these will work for your case and application, please use your discretion when applying them).

Suggested workarounds:

  • Using an interpolation other than Nearest Neighbor. Look into interpolation options and choose a different one.
  • Choose a fixed pixel size when interpolating with GDAL Warp. You can choose this based on the location you mainly want the map to be as correct as possible.
  • Choose a different output projection or use the original projection.
Related Question