6

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.

15
  • Which projection is used for the source data (grib2 file ?). Is output.tif resulting from gdal_translate correctly projected ? Also, you may have a look on gis.stackexchange.com/questions/299954/…
    – simo
    Commented Dec 15, 2021 at 5:11
  • @simo I do not know the projection of the source file, but it might be because I don't know what I am looking at. I included a link to the gdalinfo for the source file. I see EPSG 9001, 9122, and 9122 listed. I am unsure which one is the actual projection, or if they all mean the same thing?
    – David
    Commented Dec 15, 2021 at 5:16
  • Try what happens if you edit the output from step 1 with gdal_edit -mo AREA_OR_POINT=Point Output.tif --config GTIFF_POINT_GEO_IGNORE YES and then repeat the next steps.
    – user30184
    Commented Dec 15, 2021 at 9:07
  • Please try that edit. I suspect that the issue comes from setting the anchor point either into the centre of a pixel or into the top left corner of the pixel.
    – user30184
    Commented Dec 15, 2021 at 16:59
  • 1
    @FaridCheraghi Thank you for noticing, I have edited the question with the data. You may also get it here : drive.google.com/file/d/14v1s2fvWKxCg_KlJ0pI8SHxRcCNB098K/…
    – David
    Commented Dec 18, 2021 at 8:39

3 Answers 3

2

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.
1
+350

Usually, data stored in a raster format is not meant to be visualized at a scale where the pixels are visible. I therefore do not consider this problem as a geometric accuracy issue but as a precision issue.

Because of the reprojection, your squares pixels in the original coordinate system will become rectangles in the Mercator projection ("vertically" stretched). The "vertical" size of your pixels will not be constant because it depends on the latitude, therefore you cannot precisely reproject your dataset with pixels of constant size. As far as I know, Mapbox requires sqaure pixels, so no matter what tool you use, the information content will be modified. Now, let's see how we can minimize the distortions of the information.

When the size of the pixel changes, you need to assign a new value to the projected pixels. As mentioned in the previous post, you could replace the "nearest neighbor" (default) algorythm with an algorithm that combines the values of the neighbours. I recommend -r bilinear or (for a smoother result) -r cubic. Note that you seem to have values of -999 that are in fact no data, so in order to avoid problems around the boundaries, define those values as nodata (using gdal_edit.py).

You can minimize the shift by forcing the pixel size based on the X-width of your pixels, using the option -tr 1111 1111 (to keep square pixels). This will not solve vertical alignment issues, but it will be ok horizontally. Note that if you work only on a small specific area, you can select a compromise between your X and Y values on this area.

Another workaround is to artificially increase the resolution of your raster, so that you can rebuild the "rectangles" to better match their footprint in the new projection (then using again the nearest neighbor resampling will preserve the original values of the raster). The drawback is that it will increase the size of your dataset : if you take 111 m for the size of your pixels, you multiply the (uncompressed) size of your raster by 100.

enter image description here

0

Because you are not sure of the source GRIB projection format, I will give you the advice to try specifying the projection for source file when using gdalwarp.

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

OR

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

See projection extent validity on https://epsg.io/ website :

6
  • Unfortunately, these do not work - the output projection is exactly the same as my problem example image.
    – David
    Commented Dec 15, 2021 at 7:03
  • from where did you get your GRIB files ?
    – simo
    Commented Dec 15, 2021 at 7:17
  • The GRIB2 files come from this directory with live data : mrms.ncep.noaa.gov/data/2D/MergedReflectivityAtLowestAltitude ... it is part of the MRMS (Multi-Radar-Multi-Sensor) project at nssl.noaa.gov/projects/mrms . The second link has information but I see nothing about projection.
    – David
    Commented Dec 15, 2021 at 7:32
  • I guess gdal should anyway understand .... Have you tried to convert the grib file from another converter such as mygeodata.cloud/converter/grib-to-geotiff or directly read GRIB file in a GIS ?
    – simo
    Commented Dec 15, 2021 at 8:37
  • 3
    In the WKT that gdalinfo shows EPSG:9001 means metre epsg.org/unit_9001/metre.html and EPSG:9122 means degree. They are codes for units, not for coordinate reference system so using them as -s_srs is wrong.
    – user30184
    Commented Dec 16, 2021 at 7:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.