TransWikia.com

Unable to replicate GDAL output

Geographic Information Systems Asked on November 1, 2021

I have a set of GRIB files that are fliped (longitude spans from 0 to 365), and I am using gdal to first transform the data to GeoTIFF, and then warp the gridded data to a standard WGS84 longitude (-180 to 180). So far, I have been using a combination of gdal_translate and gdalwarp from the command line, and using parallel to go through all the files fast. These are the functions in my bash script:

gdal_multiband_transform(){
    FILEPATH=$1
    SAVEPATH=$2

    NUM_BANDS=$(gdalinfo $FILEPATH | grep 'Band' | wc -l)

    if [[ $NUM_BANDS -eq 1 ]]
    then
        echo "Extracting 1 band from $FILEPATH"
        gdal_translate -of GTiff -b 1 $FILEPATH $SAVEPATH
    else
        echo "Extracting 2 bands from $FILEPATH"
        gdal_translate -of GTiff -b 1 -b 2 $FILEPATH $SAVEPATH
    fi
}

warp_raster(){

    echo "Rewarp all rasters in $PATH_TO_GRIB"
    find $PATH_TO_GRIB -type f -name '*.tif' | parallel -j 5 -- gdalwarp -t_srs WGS84 {} {.}_warp.tif 
        -wo SOURCE_EXTRA=1000 --config CENTER_LONG 0 -overwrite
}

warp_raster

Now, I wanted to replicate this same behavior in Python using the osgeo library. I skipped the translation part since I realize osgeo.gdal can warp the GRIB file directly instead of having to convert/translate to a GeoTIFF format. For that I used the following Python code:

from osgeo import gdal

OPTS = gdal.WarpOptions(dstSRS='WGS84',
                        warpOptions=['SOURCE_EXTRA=1000'],
                        options=['CENTER_LONG 0'])

       
try:
   ds = gdal.Open(filename)
except RuntimeError:
   ds = gdal.Open(str(filename))

if os.path.getsize(filename):
   ds_transform = gdal.Warp(file_temp_path,
                             ds,
                             options=OPTS)
   # is this a hack?
   ds_transform = None
else:
   print(f'{filename} is an empty file. No GDAL transform')

Where I define the same option from my bash script using the gdal.WarpOptions. The outcome is visually the same; the code achieves the main goal: warp the longitude between -180 and 180. But, when I take local statistics the differences are wild. Just a mean of the whole gridded data has a difference of 4 celsius (is surface temperature data). There’s any GDAL option I am missing in osgeo that is generating these differences? I do not want to use a bash script since I am looking for an only-Python implementation.

Also posted in SO

One Answer

I was struggling with the same problem in trying to get information out of NOAA's GFS model, and finally figured out that the issue in this situation is that CENTER_LONG 0 is being passed in to the wrong set of options.

The command-line --config arguments are global GDAL configuration settings, not warp-specific options. Therefore, to replicate this functionality in Python, you need to call:

gdal.SetConfigOption('CENTER_LONG', '0')

After that, the gdal.Warp() call will work correctly.

Answered by Arktronic on November 1, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP