title: python: perform gdalwarp in memory with gdal bindings
body: <p>I currently have a processing chain in <code>R</code> which downloads MODIS data and then calls <code>gdalwarp</code> from the system to reproject a specific subdataset (e.g. NDVI) into WGS1984. The resulting <code>GeoTiffs</code> are then collected into an <code>HDF5</code> file for further processing. </p>
  
  <p>Now I'm moving the processing chain to <code>python</code>, and I was wondering if there's a way to skip the step of writing <code>GeoTiffs</code> to disk with the functionalities of the <code>gdal</code> module.</p>
  
  <p>To be clear, the question is: </p>
  
  <p>Can i perform <code>gdalwarp</code> with using strictly the python bindings from the <code>gdal</code> module and without writing to disk? </p>
  
  <p>I've been researching a bit and the closest answers to my questions were these posts:</p>
  
  <ul>
  <li><a href="https://stackoverflow.com/questions/10454316/how-to-project-and-resample-a-grid-to-match-another-grid-with-gdal-python">How to project and resample a grid to match another grid with GDAL python?
  </a></li>
  <li><a href="https://gis.stackexchange.com/questions/139906/replicating-result-of-gdalwarp-using-gdal-python-bindings">Replicating result of gdalwarp using gdal Python bindings</a></li>
  </ul>
  
  <p>The first method requires a template, so not really what I'm looking for.</p>
  
  <p>The second method looks more promising, it's using the function <code>AutoCreateWarpedVRT</code> 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).</p>
  
  <p>In my previous implementation which calls <code>gdalwarp</code> 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.</p>
  
  <p>Here's what I tried (sorry, not reproducible without the MODIS data):</p>
  
  <pre><code>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 --&gt; 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
  </code></pre>
  
  <p>And this is for the reference:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>The comparison:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>So you can see, using the <code>gdal.AutoCreateWarpedVRT</code> function results in a dataset with different dimensions and resolution. </p>