python - 无法复制 gdal 输出

Question

我有一组翻转的 GRIB 文件(经度从 0 到 365)，我使用 gdal 首先将数据转换为 GeoTIFF，然后将网格化数据扭曲为标准 WGS84 经度(-180 到 180)。到目前为止，我一直在命令行中结合使用 gdal_translate 和 gdalwarp，并使用 parallel 快速浏览所有文件。这些是我的 bash 脚本中的函数:

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

现在，我想使用 osgeo 库在 Python 中复制同样的行为。我跳过了翻译部分，因为我意识到 osgeo.gdal 可以直接扭曲 GRIB 文件，而不必转换/翻译为 GeoTIFF 格式。为此，我使用了以下 Python 代码:

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')

我在 bash 脚本中使用 gdal.WarpOptions 定义相同的选项。结果在视觉上是一样的；代码实现了主要目标:在 -180 和 180 之间扭曲经度。但是，当我进行本地统计时，差异是巨大的。整个网格数据的平均值相差 4 摄氏度(是表面温度数据)。我在 osgeo 中缺少任何产生这些差异的 GDAL 选项吗？我不想使用 bash 脚本，因为我正在寻找纯 Python 实现。

Answer 1

• 第一个选项，获取解决此问题的 GDAL 3.4，当从 GRIB 转换为 GeoTIFF 时，GRIB 会自动从 0-360 转换为 -180-180
• 第二个选项，使用 NPM 提供的 geosub (npm -g install geosub) 下载 NOAA 的 GRIB，如果您正在使用它，它可以为您完成此操作
• 第三个选项，使用从早期就存在的gdalwarp --config CENTER_LONG 0

(免责声明:我是 GDAL 中 GRIB 0-360 翻译和 geosub 包的作者)