20

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

Improve this question
2
  • 1
    Have you tried gdal.AutoCreateWarpedVRT(source_file, source_srs_wkt, dest_srs_wkt)? Commented Mar 23, 2015 at 7:21
  • Thanks Luke, did not know this function. Tried just now, but some pixels are still different between the two. I.e., the geo-transforms and shapes of the rasters are identical (when not up-sampled), but some pixels seem to be resampled differently. This at least demonstrates the issue is still present even when not up-sampling. Commented Mar 23, 2015 at 15:20

1 Answer 1

Reset to default
17

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"
Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.