Das Ergebnis von gdalwarp wird mit gdal Python-Bindungen repliziert

20

Ich versuche, die GDAL-Pythonbindungen neu zu projizieren / zu resampeln, erhalte jedoch geringfügig andere Ergebnisse als mit dem Befehlszeilendienstprogramm gdalwarp.

Siehe Update unten für ein kürzeres Beispiel

Dieses Skript veranschaulicht den Python-Ansatz:

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

Vergleichen Sie mit Ausgabe von:

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

Sie unterscheiden sich in der Größe (um 2 Zeilen und 1 Spalte) sowie mit einigen unterschiedlichen Pixelwerten in der Nähe der Kanten.

Siehe transparente Überlagerung von test_epsg3573.tif und test_epsg3573_gdalwarp.tif weiter unten. Wenn die Bilder identisch wären, gäbe es nur schwarze und weiße Pixel, kein Grau.

QGIS-Overlay von test_epsg3573.tif und test_epsg3573_gdalwarp.tif

Getestet mit Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1

Update :

Hier ist ein viel kürzeres Beispiel. Dies scheint nicht auf ein Upsampling zurückzuführen zu sein, da das Folgende ebenfalls zu widersprüchlichen Ergebnissen führtgdalwarp

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

Und dies ist der GDALWARP-Aufruf, von dem ich erwarte, dass er derselbe ist, aber nicht:

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

Das folgende Bild zeigt jedes resultierende Binärbild, das mit 50% Transparenz überlagert ist. Die hellgrauen Pixel sind Inkonsistenzen zwischen den beiden Ergebnissen.

Inkonsistenz in QGIS dargestellt

python gdal gdalwarp numpy Bruce Wallin
quelle
1
Hast du es versucht gdal.AutoCreateWarpedVRT(source_file, source_srs_wkt, dest_srs_wkt)?
user2856
Danke Luke, kannte diese Funktion nicht. Versuchte gerade, aber einige Pixel sind immer noch unterschiedlich zwischen den beiden. Das heißt, die Geotransformationen und Formen der Raster sind identisch (wenn sie nicht aufgesampelt werden), aber einige Pixel scheinen anders aufgesampelt zu werden. Dies zeigt zumindest, dass das Problem auch dann noch besteht, wenn kein Upsampling durchgeführt wird.
Bruce Wallin

Antworten:

16

Ich erhalte die gleichen Ergebnisse wie gdalwarpaus , gdal.AutoCreateWarpedVRTwenn ich die Fehlerschwelle auf 0,125 gesetzt den Standard entsprechen (-et) in gdalwarp . Alternativ können Sie -et 0.0in Ihrem Anruf festlegen gdalwarp, dass der Standardwert in übereinstimmt gdal.AutoCreateWarpedVRT.

Beispiel

Erstellen Sie eine Referenz zum Vergleichen mit:

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

Führen Sie die Projektion in Python aus (basierend auf Code aus der Funktion "warp_27 ()" in der GDAL- Autotest-Suite ):

# 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
quelle