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/new_posting)) # Raster xsize y_size = DEM.RasterYSize * _np.int(_np.abs(geo_t/new_posting))# Raster ysize #Compute new geotransform new_geo = ( geo_t, new_posting, geo_t, \ geo_t, geo_t, new_posting ) 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?