Using gdal I want to resample (in this case, oversample) a DEM. The only thing that should change is the pixel spacing, so I want to interpolate in between the values that I currently have. For this I have the current function:

def resample_DEM_and_image(DEM,new_posting):
    This function reprojects a gtiff
    to a new posting

    DEM : osgeo.gdal.Dataset
    The geotiff to project
    new_posting : iterable
    The desired new pixel spacing
    from osgeo import gdal, osr
    geo_t = DEM.GetGeoTransform()

    x_size = DEM.RasterXSize * _np.int(_np.abs(geo_t[1]/new_posting[0])) # Raster xsize
    y_size = DEM.RasterYSize * _np.int(_np.abs(geo_t[5]/new_posting[1]))# Raster ysize

    #Compute new geotransform
    new_geo = ( geo_t[0], new_posting[0], geo_t[2], \
            geo_t[3], geo_t[4], new_posting[1] )
    mem_drv = gdal.GetDriverByName( 'MEM' )
    dest = mem_drv.Create('', x_size, y_size, DEM.RasterCount,1,gdal.GDT_Float32)
    res = gdal.ReprojectImage( DEM, dest, \
            DEM.GetProjection(), DEM.GetProjection(), \
            gdal.GRA_Bilinear )
    return dest

then I was planning to read the gdal dataset that I get in dest using dest.ReadAsArray(). However all I get is values of 0 in an array that does have the correct shape (i.e. n times the size of the initial DEM, where n is the oversampling factor)

what am I doing wrong?

2 Answers 2


should RasterCount be a parameter within mem_drv.Create? something like:

dest = mem_drv.Create('', x_size, y_size, 1, gdal.GDT_Float32)

EDIT: oh, i see it here - i didn't see the destination raster as a parameter to ReprojectImage at first - hopefully the above change will help. i was originally thinking of re-interpolating the array, then writing it to the dest raster band -

band = dest.GetRasterBand(1)
res = SomeReInterpolationFunction(**kwargs) # if this were a custom function to interpolate an array

though it doesn't appear to be needed in this case.

  • Thanks for your comment. I got the function from a colleague who is out of office at the time, which is why I'm asking here. Originally this line of code was actually : dest = mem_drv.Create('', x_size, y_size, DEM.RasterCount,numeric_dt_to_gdal_dt(DEM.GetRasterBand(1).DataType)). After I found the page you suggested I changed it to what I had and forgot to remove the RasterCount, but in any case, the result is the same with any of these 3 option...
    – uetli
    Oct 16, 2015 at 12:14
  • hmm.. i tried it on a geotiff i had lying around and it seemed to work - do you have a sample dataset? Oct 16, 2015 at 13:42
  • what are your new_posting values? Oct 17, 2015 at 11:37
  • 3
    i was initially using 0.5, 0.5 - though the original geotransform has a positive x and negative y. After changing to new_posting values of 0.5, -0.5 it seemed to work for me. Oct 17, 2015 at 12:14
  • great, that solved the problem. It totally makes sense, but didn't think about that. Thanks a lot
    – uetli
    Oct 19, 2015 at 6:17

The new_posting needs to include the sign of the old geotransform, i.e. negative for the y-direction and positive in the x-direction in this case.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.