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
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
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP