时间戳

首页 > 解决方案 > 无法复制 gdal 输出

python - 无法复制 gdal 输出

问题描述

我有一组翻转的 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 的实现。

标签： pythongeospatialgdalosgeo

解决方案

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

（免责声明：我是GDAL和geosub包中GRIB 0-360翻译的作者）

推荐阅读