Resampling DEM using gdal in python?

Question

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

    Parameters
    ----------
    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)
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(DEM.GetProjection())
    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?

Answer 1

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
band.WriteArray(res)

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

Answer 2

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.