How to use the gdal.ReprojectImage function in GDAL

ndv = src_b.GetNoDataValue()
        dst_b = dst_ds.GetRasterBand(i)
        if ndv is not None:
            dst_b.SetNoDataValue(ndv)
            dst_b.Fill(ndv)

    if interpolation == 'near':
        gdal_mode = gdal.GRA_NearestNeighbour
    elif interpolation == 'mode':
        gdal_mode = gdal.GRA_Mode
    elif interpolation == 'average':
        gdal_mode = gdal.GRA_Average
    else:
        raise ValueError("Invalid interpolation mode %s" % interpolation)

    res = gdal.ReprojectImage(
        src_ds,
        dst_ds,
        src_ds.GetProjectionRef(),
        dst_ds.GetProjectionRef(),
        gdal_mode
    )
    assert res == 0, 'Error reprojecting, res=%d' % res
    dst_arr = dst_ds.ReadAsArray()

    # GDAL ReadAsArray returns (bands, height, width) but we expect
    # (height, width, bands)
    if len(dst_arr.shape) == 3:
        dst_arr = dst_arr.transpose(1, 2, 0)

    # TODO: This assumes that the no data value is the same for all bands
    if ndv is not None:

high = match_ds.RasterYSize
        
        # Output / destination
        try:
            # try create File driver based in path
            dst = gdal.GetDriverByName('GTiff').Create(output, wide, high, 1, gdalconst.GDT_Float32)
        except RuntimeError:
            raise GeoAlgorithmExecutionException("Could not generate output file")
            
        dst.SetGeoTransform( match_geotrans )
        dst.SetProjection( match_proj)
        dst.GetRasterBand(1).SetNoDataValue(nodata) # write old nodata value

        # Do the work
        if interp == 'Bilinear':
            gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Bilinear)
        elif interp == 'Cubic':
            gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Cubic)        
        elif interp == 'Cubicspline':
            gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_CubicSpline)            
        elif interp == 'Lanczos':
            gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_Lanczos)            
        elif interp == 'NearestNeighbour':
            gdal.ReprojectImage(src, dst, src_proj, match_proj, gdalconst.GRA_NearestNeighbour)
                    
        del dst # Flush

im = Image.fromarray(array, 'RGBA')  # Always four bands
            im1 = im.resize((tilesize, tilesize), Image.ANTIALIAS)
            if path.exists(tilefilename):
                im0 = Image.open(tilefilename)
                im1 = Image.composite(im1, im0, im1)
            im1.save(tilefilename, self.tiledriver)

        else:

            # Other algorithms are implemented by gdal.ReprojectImage().
            dsquery.SetGeoTransform((0.0, tilesize / float(querysize), 0.0,
                                     0.0, 0.0, tilesize / float(querysize)))
            dstile.SetGeoTransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0))
            
            
            res = gdal.ReprojectImage(dsquery, dstile, None, None,
                                      self.resampling)
            if res != 0:
                self.error("ReprojectImage() failed on %s, error %d" %
                           (tilefilename, res))

dataset = get_dataset(raster)
    datatype = dataset.GetRasterBand(1).DataType
    resized_ds = gdal.GetDriverByName('MEM').Create(
        '', dimensions[0], dimensions[1],  dataset.RasterCount, datatype)
    for i in range(1, resized_ds.RasterCount+1):
        nodatavalue = dataset.GetRasterBand(i).GetNoDataValue()
        resized_band = resized_ds.GetRasterBand(i)
        resized_arr = resized_band.ReadAsArray()
        if nodatavalue:
            resized_arr[resized_arr == 0] = nodatavalue
            resized_band.SetNoDataValue(nodatavalue)
        resized_band.WriteArray(resized_arr)
    resized_ds.SetGeoTransform(transform)
    resized_ds.SetProjection(projection)

    gdal.ReprojectImage(dataset, resized_ds)
    return resized_ds

modis = osr.SpatialReference()
    modis.ImportFromProj4(wkt_from)
    tx = osr.CoordinateTransformation(modis, wgs84)
    g = gdal.Open(dataset)
    geo_t = g.GetGeoTransform()
    print geo_t
    x_size = g.RasterXSize
    y_size = g.RasterYSize
    (ulx, uly, ulz) = tx.TransformPoint(geo_t[0], geo_t[3])
    (lrx, lry, lrz) = tx.TransformPoint(geo_t[0] + (geo_t[1]*x_size), geo_t[3]+ (geo_t[5]*y_size))
    mem_drv = gdal.GetDriverByName(file_format)
    dest = mem_drv.Create(output_file_name, int((lrx-ulx)/pixel_spacing), int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)
    new_geo = ([ulx, pixel_spacing, geo_t[2], uly, geo_t[4], -pixel_spacing])
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(wgs84.ExportToWkt())
    gdal.ReprojectImage(g, dest, modis.ExportToWkt(), wgs84.ExportToWkt(), gdal.GRA_Bilinear)
    print "reprojected"

dataset = get_dataset(raster)
    datatype = dataset.GetRasterBand(1).DataType
    resized_ds = gdal.GetDriverByName('MEM').Create(
        '', dimensions[0], dimensions[1],  dataset.RasterCount, datatype)
    for i in range(1, resized_ds.RasterCount+1):
        nodatavalue = dataset.GetRasterBand(i).GetNoDataValue()
        resized_band = resized_ds.GetRasterBand(i)
        resized_arr = resized_band.ReadAsArray()
        if nodatavalue:
            resized_arr[resized_arr == 0] = nodatavalue
            resized_band.SetNoDataValue(nodatavalue)
        resized_band.WriteArray(resized_arr)
    resized_ds.SetGeoTransform(transform)
    resized_ds.SetProjection(projection)

    gdal.ReprojectImage(dataset, resized_ds)
    return resized_ds

print('Correct DEM for geoid')
        self._add_egm96(dem_tiff, egm_tiff, data_folder)

        if resample == 'regular_grid':
            # If regular grid is used, we convert using gdal.
            # First create a geotiff file, then resample geotiff file. We always use cubic interpolation.
            print('Resampling to new regular grid')
            dlat = lats[1] - lats[0]
            dlon = lons[1] - lons[0]
            dem_tiff_final = os.path.join(data_folder, 'dem.tiff')
            dem_data_final = self._create_georeference(latlim, lonlim, dlat, dlon, 'float32', dem_tiff_final)
            dem_data_final = None
            dem_data_final = gdal.Open(dem_tiff_final, gdal.GA_Update)
            dem_data = gdal.Open(dem_tiff, gdal.GA_Update)

            gdal.ReprojectImage(dem_data, dem_data_final, None, None, gdal.GRA_Cubic)
            dem_data_final = None
            dem_tiff = dem_tiff_final

        elif resample == 'irregular_grid':
            # Use a simple bilinear approach to find values for specific points.
            print('Resampling to new irregular grid')
            heights = self._simple_bilinear(lats, lons, dem_tiff)

            return heights

        if doris_input == True:
            # Create a binary output file
            command = 'gdal_translate -of MFF ' + dem_tiff + ' ' + dem_tiff[:-5] + '.raw'
            os.system(command)

            if not os.path.exists(os.path.dirname(out_file)):

tilesize / float(querysize),
                0.0,
                0.0,
                0.0,
                tilesize / float(querysize),
                ))
            dstile.SetGeoTransform((
                0.0,
                1.0,
                0.0,
                0.0,
                0.0,
                1.0,
                ))

            res = gdal.ReprojectImage(dsquery, dstile, None, None,
                    self.resampling)
            if res != 0:
                self.error('ReprojectImage() failed on %s, error %d'
                           % (tilefilename, res))

col=gland.RasterXSize
    rows=gland.RasterYSize

    # Create new raster
    mem_drv = gdal.GetDriverByName('MEM')
    dest1 = mem_drv.Create('', col, rows, 1, gdal.GDT_Float32)
    dest1.SetGeoTransform(geo_land)
    dest1.SetProjection(osng.ExportToWkt())

    # Perform the projection/resampling
    if method is 1:
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_NearestNeighbour)
    if method is 2:
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Bilinear)
    if method is 3:
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Lanczos)
    if method is 4:
        gdal.ReprojectImage(g, dest1, wgs84.ExportToWkt(), osng.ExportToWkt(), gdal.GRA_Average)
    return(dest1)

querysize = dsquery.RasterXSize
        tilesize = dstile.RasterXSize
        tilebands = dstile.RasterCount

        if self.options.resampling == 'average':
            for i in range(1,tilebands+1):
                res = gdal.RegenerateOverview( dsquery.GetRasterBand(i),
                    dstile.GetRasterBand(i), 'average' )
                if res != 0:
                    self.error("RegenerateOverview() failed")
        else:
            # Other algorithms are implemented by gdal.ReprojectImage().
            dsquery.SetGeoTransform( (0.0, tilesize / float(querysize), 0.0, 0.0, 0.0, tilesize / float(querysize)) )
            dstile.SetGeoTransform( (0.0, 1.0, 0.0, 0.0, 0.0, 1.0) )
            
            res = gdal.ReprojectImage(dsquery, dstile, None, None, self.resampling)    
            if res != 0:
                self.error("ReprojectImage() failed on %s, error %d" % (tilefilename, res))