python: perform gdalwarp in memory with gdal bindings

Question

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.

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.

To be clear, the question is:

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:

• 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.

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

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.

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

And this is for the reference:

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

The comparison:

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

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

Answer 1

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

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)

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.