Réplication du résultat de gdalwarp à l'aide des liaisons gdal Python

20

J'essaie de re-projeter / rééchantillonner avec les liaisons Python GDAL, mais j'obtiens des résultats légèrement différents par rapport à ceux de l'utilitaire de ligne de commande gdalwarp.

Voir la mise à jour ci-dessous pour un exemple plus court

Ce script illustre l'approche Python:

from osgeo import osr, gdal
import numpy


def reproject_point(point, srs, target_srs):
    '''
    Reproject a pair of coordinates from one spatial reference system to
    another.
    '''
    transform = osr.CoordinateTransformation(srs, target_srs)
    (x, y, z) = transform.TransformPoint(*point)

    return (x, y)


def reproject_bbox(top_left, bottom_right, srs, dest_srs):
    x_min, y_max = top_left
    x_max, y_min = bottom_right
    corners = [
        (x_min, y_max),
        (x_max, y_max),
        (x_max, y_min),
        (x_min, y_min)]
    projected_corners = [reproject_point(crnr, srs, dest_srs)
                         for crnr in corners]

    dest_top_left = (min([crnr[0] for crnr in projected_corners]),
                     max([crnr[1] for crnr in projected_corners]))
    dest_bottom_right = (max([crnr[0] for crnr in projected_corners]),
                         min([crnr[1] for crnr in projected_corners]))

    return dest_top_left, dest_bottom_right


################################################################################
# Create synthetic data
gtiff_drv = gdal.GetDriverByName('GTiff')
w, h = 512, 512
raster = numpy.zeros((w, h), dtype=numpy.uint8)
raster[::w / 10, :] = 255
raster[:, ::h / 10] = 255
top_left = (-109764, 215677)
pixel_size = 45

src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(3413)

src_geotran = [top_left[0], pixel_size, 0,
               top_left[1], 0, -pixel_size]

rows, cols = raster.shape
src_ds = gtiff_drv.Create(
    'test_epsg3413.tif',
    cols, rows, 1,
    gdal.GDT_Byte)
src_ds.SetGeoTransform(src_geotran)
src_ds.SetProjection(src_srs.ExportToWkt())
src_ds.GetRasterBand(1).WriteArray(raster)


################################################################################
# Reproject to EPSG: 3573 and upsample to 7m
dest_pixel_size = 7

dest_srs = osr.SpatialReference()
dest_srs.ImportFromEPSG(3573)

# Calculate new bounds by re-projecting old corners
x_min, y_max = top_left
bottom_right = (x_min + cols * pixel_size,
                y_max - rows * pixel_size)
dest_top_left, dest_bottom_right = reproject_bbox(
    top_left, bottom_right,
    src_srs, dest_srs)

# Make dest dataset
x_min, y_max = dest_top_left
x_max, y_min = dest_bottom_right
new_rows = int((x_max - x_min) / float(dest_pixel_size))
new_cols = int((y_max - y_min) / float(dest_pixel_size))
dest_ds = gtiff_drv.Create(
    'test_epsg3573.tif',
    new_rows, new_cols, 1,
    gdal.GDT_Byte)
dest_geotran = (dest_top_left[0], dest_pixel_size, 0,
                dest_top_left[1], 0, -dest_pixel_size)
dest_ds.SetGeoTransform(dest_geotran)
dest_ds.SetProjection(dest_srs.ExportToWkt())

# Perform the projection/resampling
gdal.ReprojectImage(
    src_ds, dest_ds,
    src_srs.ExportToWkt(), dest_srs.ExportToWkt(),
    gdal.GRA_NearestNeighbour)

dest_data = dest_ds.GetRasterBand(1).ReadAsArray()

# Close datasets
src_ds = None
dest_ds = None

Comparer avec la sortie de:

gdalwarp -s_srs EPSG:3413 -t_srs EPSG:3573 -tr 7 7 -r near -of GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif

Ils diffèrent en taille (par 2 lignes et 1 colonne) ainsi que par certaines valeurs de pixels différentes près des bords.

Voir la superposition transparente de test_epsg3573.tif et test_epsg3573_gdalwarp.tif ci-dessous. Si les images étaient identiques, il n'y aurait que des pixels noir et blanc, pas de gris.

Superposition QGIS de test_epsg3573.tif et test_epsg3573_gdalwarp.tif

Testé avec Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1

Mise à jour :

Voici un exemple beaucoup plus court. Cela ne semble pas être causé par un suréchantillonnage, car les résultats suivants produisent également des résultats incompatibles avecgdalwarp

from osgeo import osr, gdal
import numpy


# Create synthetic data
gtiff_drv = gdal.GetDriverByName('GTiff')
w, h = 512, 512
raster = numpy.zeros((w, h), dtype=numpy.uint8)
raster[::w / 10, :] = 255
raster[:, ::h / 10] = 255
top_left = (-109764, 215677)
pixel_size = 45

src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(3413)

src_geotran = [top_left[0], pixel_size, 0,
               top_left[1], 0, -pixel_size]

rows, cols = raster.shape
src_ds = gtiff_drv.Create(
    'test_epsg3413.tif',
    cols, rows, 1,
    gdal.GDT_Byte)
src_ds.SetGeoTransform(src_geotran)
src_ds.SetProjection(src_srs.ExportToWkt())
src_ds.GetRasterBand(1).WriteArray(raster)

# Reproject to EPSG: 3573
dest_srs = osr.SpatialReference()
dest_srs.ImportFromEPSG(3573)

int_ds = gdal.AutoCreateWarpedVRT(src_ds, src_srs.ExportToWkt(), dest_srs.ExportToWkt())

# Make dest dataset
dest_ds = gtiff_drv.Create(
    'test_epsg3573_avrt.tif',
    int_ds.RasterXSize, int_ds.RasterYSize, 1,
    gdal.GDT_Byte)
dest_ds.SetGeoTransform(int_ds.GetGeoTransform())
dest_ds.SetProjection(int_ds.GetProjection())
dest_ds.GetRasterBand(1).WriteArray(int_ds.GetRasterBand(1).ReadAsArray())

# Close datasets
src_ds = None
dest_ds = None

Et c'est l'appel gdalwarp que je m'attends à être le même, mais ce n'est pas le cas:

gdalwarp -s_srs EPSG:3413 -t_srs EPSG:3573 -of GTiff test_epsg3413.tif test_epsg3573_gdalwarp.tif

L'image ci-dessous montre chaque image binaire résultante superposée à 50% de transparence. Les pixels gris clair sont des incohérences entre les deux résultats.

Incohérence illustrée dans QGIS

python gdal gdalwarp numpy Bruce Wallin
la source
1
As-tu essayé gdal.AutoCreateWarpedVRT(source_file, source_srs_wkt, dest_srs_wkt)?
user2856
Merci Luke, ne connaissait pas cette fonction. J'ai essayé tout à l'heure, mais certains pixels sont toujours différents entre les deux. C'est-à-dire que les géotransformations et les formes des rasters sont identiques (lorsqu'elles ne sont pas suréchantillonnées), mais certains pixels semblent être rééchantillonnés différemment. Cela démontre au moins que le problème est toujours présent même lorsqu'il n'est pas suréchantillonné.
Bruce Wallin

Réponses:

16

Je reçois les mêmes résultats que gdalwarpde gdal.AutoCreateWarpedVRTsi je mets le seuil d'erreur à 0,125 pour correspondre à la valeur par défaut (-et) dans gdalwarp . Vous pouvez également définir -et 0.0votre appel pour gdalwarpqu'il corresponde à la valeur par défaut de gdal.AutoCreateWarpedVRT.

Exemple

Créez une référence à comparer:

gdalwarp -t_srs EPSG:4326 byte.tif warp_ref.tif

Exécutez la projection en Python (basé sur le code de la fonction "warp_27 () dans la suite d'autotests GDAL ):

# Open source dataset
src_ds = gdal.Open('byte.tif')

# Define target SRS
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(4326)
dst_wkt = dst_srs.ExportToWkt()

error_threshold = 0.125  # error threshold --> use same value as in gdalwarp
resampling = gdal.GRA_NearestNeighbour

# Call AutoCreateWarpedVRT() to fetch default values for target raster dimensions and geotransform
tmp_ds = gdal.AutoCreateWarpedVRT( src_ds,
                                   None, # src_wkt : left to default value --> will use the one from source
                                   dst_wkt,
                                   resampling,
                                   error_threshold )

# Create the final warped raster
dst_ds = gdal.GetDriverByName('GTiff').CreateCopy('warp_test.tif', tmp_ds)
dst_ds = None

# Check that we have the same result as produced by 'gdalwarp -rb -t_srs EPSG:4326 ....'

ref_ds = gdal.Open('warp_ref.tif')
ref_cs = ref_ds.GetRasterBand(1).Checksum()

ds = gdal.Open('warp_test.tif')
cs = ds1.GetRasterBand(1).Checksum()

if cs == ref_cs:
    print 'success, they match'
else:
    print "fail, they don't match"
user2856
la source