Python：使用GDAL绑定在内存中执行gdalwarp

Question

Python：使用GDAL绑定在内存中执行gdalwarp

pythongisgdal

12

我目前在R中有一个处理链，它下载MODIS数据，然后从系统调用gdalwarp重新投影特定的子数据集（例如NDVI）到WGS1984。然后将生成的GeoTiffs收集到HDF5文件中进行进一步处理。

现在我正在将处理链移植到python，想知道是否有一种方法可以使用gdal模块的功能跳过将GeoTiffs写入磁盘的步骤。

明确一下问题：

我是否能够仅使用gdal模块的python绑定执行gdalwarp而不写入磁盘？

我进行了一些研究，最接近我的问题的答案是这些帖子：

第一种方法需要一个模板，所以不是我要找的。

第二种方法看起来更有前途，它使用了AutoCreateWarpedVRT函数，这似乎是我想要的。虽然与答案中的示例相反，我的结果与参考结果不匹配（与任何误差阈值无关）。

在直接调用gdalwarp的以前实现中，除了目标参考系统外，我还指定了目标分辨率。因此，我认为这可能是有所区别的 - 但我无法在python中通过gdal绑定设置它。

以下是我的尝试（抱歉，没有MODIS数据无法重现）：

import gdal
import osr

ds = gdal.Open('/data/MOD13A2.A2016305.h18v07.005.2016322013359.hdf')

t_srs = osr.SpatialReference()
t_srs.ImportFromEPSG(4326)

src_ds = gdal.Open(ds.GetSubDatasets()[0][0], gdal.GA_ReadOnly)

dst_wkt =t_srs.ExportToWkt()
error_threshold = 0.125
resampling=gdal.GRA_NearestNeighbour

tmp_ds = gdal.AutoCreateWarpedVRT( src_ds,
                                  None, # src_wkt : left to default value --> will use the one from source
                                   dst_wkt,
                                   resampling,
                                   error_threshold)

# create tiff

dst_ds = gdal.GetDriverByName('GTiff').CreateCopy('warp_test.tif', tmp_ds)
dst_ds = None

这是供参考的：

gdalwarp -ot Int16 -tr 0.00892857142857143 0.00892857142857143 -t_srs EPSG:4326 "HDF4_EOS:EOS_GRID:MOD13A2.A2016305.h18v07.005.2016322013359.hdf:MODIS_Grid_16DAY_1km_VI:1 km 16 days NDVI" MOD13A2.A2016305.h18v07.005.2016322013359_MODIS_Grid_16DAY_1km_VI_1_km_16_days_NDVI.tif

比较：

i1 = gdal.Open('warp_test.tif')
i2 = gdal.Open('MOD13A2.A2016305.h18v07.005.2016322013359_MODIS_Grid_16DAY_1km_VI_1_km_16_days_NDVI.tif')

# test

print(i1.RasterXSize,i1.RasterYSize)
1267 1191

#reference

print(i2.RasterXSize,i2.RasterYSize)
1192 1120

i1.GetRasterBand(1).Checksum() == i2.GetRasterBand(1).Checksum()
False

所以您可以看到，使用gdal.AutoCreateWarpedVRT函数会导致数据集具有不同的尺寸和分辨率。

- Val

你的示例代码中的 sr、dst_wkt 和 tmp_ds2 是什么？它们无法被重现，因此很难理解问题是什么。 - Rutger Kassies

1

抱歉@Rutger，这段代码是从多个部分复制的，我应该更好地整理旧变量定义。我更新了我的问题。 - Val

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Rutger Kassies · Accepted Answer

如果你想模仿你的“参考”调用gdalwarp，你可以使用：

import gdal

ds = gdal.Warp('warp_test.tif', infile, dstSRS='EPSG:4326',
               outputType=gdal.GDT_Int16, xRes=0.00892857142857143, yRes=0.00892857142857143)
ds = None

如果您不想将输出写入磁盘文件，可以转而使用内存中的VRT文件，例如：

ds = gdal.Warp('', infile, dstSRS='EPSG:4326', format='VRT',
               outputType=gdal.GDT_Int16, xRes=0.00892857142857143, yRes=0.00892857142857143)

当然，您可以将任何格式的数据进行变形，但对于除了 VRT 以外的文件，变形后的结果实际上会存储在内存中。