python:使用gdal绑定在内存中执行gdalwarp [英] python: perform gdalwarp in memory with gdal bindings

问题描述

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

I currently have a processing chain in R which downloads MODIS data and then calls gdalwarp from the system to reproject a specific subdataset (e.g. NDVI) into WGS1984. The resulting GeoTiffs are then collected into an HDF5 file for further processing.

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

Now I'm moving the processing chain to python, and I was wondering if there's a way to skip the step of writing GeoTiffs to disk with the functionalities of the gdal module.

要清楚，问题是:

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

Can i perform gdalwarp with using strictly the python bindings from the gdal module and without writing to disk?

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

I've been researching a bit and the closest answers to my questions were these posts:

• 如何投影并重新采样网格以与GDAL python匹配另一个网格?
• 使用gdal Python绑定复制gdalwarp的结果

• How to project and resample a grid to match another grid with GDAL python?
• Replicating result of gdalwarp using gdal Python bindings

第一种方法需要模板，所以不是我真正想要的.

The first method requires a template, so not really what I'm looking for.

第二种方法看起来更有希望，它使用的功能是AutoCreateWarpedVRT，这似乎正是我想要的.尽管与答案中的示例相反，但我的结果与参考值不匹配(独立于任何错误阈值).

The second method looks more promising, it's using the function AutoCreateWarpedVRT which seems to be quite what I want. Although, in contrary to the example in the answer, my result doesn't match the reference (independently of any error threshold).

在我先前的直接调用gdalwarp的实现中，除了目标参考系统之外，我还指定了目标分辨率.因此，我认为这可能会有所作为-但我无法在python的gdal绑定中进行设置.

In my previous implementation which calls gdalwarp directly, I've specified a target resolution in addition to the target reference system. So I assume that's something that could make the difference - but I haven't been able to set it within the gdal bindings in python.

这是我尝试过的(抱歉，没有MODIS数据无法重现):

Here's what I tried (sorry, not reproducible without the MODIS data):

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函数可以得到具有不同尺寸和分辨率的数据集.

So you can see, using the gdal.AutoCreateWarpedVRT function results in a dataset with different dimensions and resolution.

推荐答案

如果要模仿对gdalwarp的引用"调用，可以使用:

If you want to mimic your "reference" call to gdalwarp you can use:

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文件，例如:

If you dont want to output to a file on disk, you can warp to an in-memory VRT file, for example:

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

您当然可以将其扭曲到内存中的任何格式，但是对于VRT以外的文件，扭曲的结果实际上将存储在内存中.

You can of course warp to any format in memory, but for files other than VRT the warped result will actually be stored in-memory.

这篇关于python:使用gdal绑定在内存中执行gdalwarp的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！