CODEWITHSUNDEEP
pythongeospatialgdalosgeo
topcat
topcat

Reputation: 586

Unable to replicate gdal output

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.

Upvotes: 3

Views: 420

Answers (1)

Criminal_Affair_At_SO
Criminal_Affair_At_SO

Reputation: 3443

  • First option, get GDAL 3.4 where this problem is solved, GRIBs get automagically transformed from 0-360 to -180-180 when being converted from GRIB to GeoTIFF
  • Second option, use geosub available from NPM (npm -g install geosub) to download NOAA's GRIBs if this is what you are using, it can do this for you
  • Third option, use gdalwarp --config CENTER_LONG 0 which has been there since the early days

(Disclaimer: I am the author of the GRIB 0-360 translation in GDAL and the geosub package)

Upvotes: 1

Related Questions