Suite

Rééchantillonner DEM à l'aide de gdal en python ?

Rééchantillonner DEM à l'aide de gdal en python ?


En utilisant gdal, je veux rééchantillonner (dans ce cas, suréchantillonner) un DEM. La seule chose qui devrait changer est l'espacement des pixels, donc je veux interpoler entre les valeurs que j'ai actuellement. Pour cela j'ai la fonction actuelle :

def resample_DEM_and_image(DEM,new_posting): """ Cette fonction reprojette un gtiff vers une nouvelle publication Paramètres ---------- DEM : osgeo.gdal.Dataset Le geotiff à projeter new_posting : iterable Le nouvel espacement de pixels souhaité """ depuis 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 #Calculer la nouvelle géotransformation 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

alors j'avais l'intention de lire l'ensemble de données gdal que j'obtiens dans dest en utilisant dest.ReadAsArray(). Cependant, tout ce que j'obtiens, ce sont des valeurs de 0 dans un tableau qui a la forme correcte (c'est-à-dire n fois la taille du DEM initial, où n est le facteur de suréchantillonnage)

Qu'est-ce que je fais mal?


RasterCount doit-il être un paramètre dans mem_drv.Create ? quelque chose comme:

dest = mem_drv.Create(", x_size, y_size, 1, gdal.GDT_Float32)

EDIT: oh, je le vois ici - je n'ai pas vu le raster de destination comme paramètre de ReprojectImage au début - j'espère que le changement ci-dessus vous aidera. je pensais à l'origine à réinterpoler le tableau, puis à l'écrire dans la bande raster de destination -

band = dest.GetRasterBand(1) res = SomeReInterpolationFunction(**kwargs) # s'il s'agissait d'une fonction personnalisée pour interpoler un tableau band.WriteArray(res)

bien que cela ne semble pas être nécessaire dans ce cas.


Lenouvelle_publicationdoit inclure le signe de l'ancienne géotransformation, c'est-à-dire négatif pour la direction y et positif dans la direction x dans ce cas.


Voir la vidéo: Reproject, resample and clip raster data with GDAL in Python