GDAL Python – How to Replicate GDALWarp Results with GDAL Python Bindings and Numpy

gdalgdalwarpnumpypython

I am trying to re-project/resample with the GDAL python bindings, but am getting slightly different results compared to those from the command line utility gdalwarp.

See update below for shorter example

This script illustrates the Python approach:

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

Compare with output of:

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

They differ in size (by 2 rows and 1 column) as well as with some differing pixel values near edges.

See transparent overlay of test_epsg3573.tif and test_epsg3573_gdalwarp.tif below. If images were identical there would only be black and white pixels, no grey.

QGIS overlay of test_epsg3573.tif and test_epsg3573_gdalwarp.tif

Tested with Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1

Update:

Here is a much shorter example. This seems to not be caused by upsampling as the following also produces results inconsistent with gdalwarp

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

And this is the gdalwarp call that I expect to be the same, yet is not:

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

The image below shows each resulting binary image overlayed at 50% transparency. The light grey pixels are inconsistencies between the two results.

Inconsistency illustrated in QGIS

Best Answer

I get the same results as gdalwarp from gdal.AutoCreateWarpedVRT if I set the error threshold to 0.125 to match the default (-et) in gdalwarp. Alternatively, you could set -et 0.0 in your call to gdalwarp to match the default in gdal.AutoCreateWarpedVRT.

Example

Create a reference to compare to:

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

Run the projection in Python (based on code from the "warp_27() function in the 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"
Related Question