I have yet again trouble 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 projectnew_posting : iterable The desired new pixel spacing"""from osgeo import gdal, osrgeo_t = DEM.GetGeoTransform()x_size = DEM.RasterXSize * _np.int(_np.abs(geo_t[1]/new_posting[0])) # Raster xsizey_size = DEM.RasterYSize * _np.int(_np.abs(geo_t[5]/new_posting[1]))# Raster ysize #Compute new geotransformnew_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 destthen 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?
thank you for your input
أكثر...
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 projectnew_posting : iterable The desired new pixel spacing"""from osgeo import gdal, osrgeo_t = DEM.GetGeoTransform()x_size = DEM.RasterXSize * _np.int(_np.abs(geo_t[1]/new_posting[0])) # Raster xsizey_size = DEM.RasterYSize * _np.int(_np.abs(geo_t[5]/new_posting[1]))# Raster ysize #Compute new geotransformnew_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 destthen 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?
thank you for your input
أكثر...