python:使用gdal绑定在内存中执行gdalwarp

Question

我现在有一个处理链中R,这将下载数据MODIS,然后调用gdalwarp从系统重新投影的特定子数据集(例如NDVI)成WGS1984.GeoTiffs然后将结果收集到HDF5文件中以供进一步处理.

现在我移动处理链python,我想知道是否有跳过写入的步骤的方式GeoTiffs与的功能到磁盘gdal模块.

需要说明的是,问题是:

我可以执行gdalwarp严格使用gdal模块的python绑定而无需写入磁盘吗？

我一直在研究,最接近我的问题的答案是这些帖子:

• 如何投影和重新采样网格以与GDAL python匹配另一个网格？
• 使用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具有不同尺寸和分辨率的数据集中的函数结果.

Answer 1

如果您想模仿您的"参考"电话,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以外的文件,扭曲的结果实际上将存储在内存中.