1

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?

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

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 '15 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? – fluidmotion Oct 16 '15 at 13:42
  • what are your new_posting values? – fluidmotion Oct 17 '15 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. – fluidmotion Oct 17 '15 at 12:14
  • great, that solved the problem. It totally makes sense, but didn't think about that. Thanks a lot – uetli Oct 19 '15 at 6:17
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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