gdalwarp causing shift in pixels

Question

Can somebody help me with gdalwarp issue? I am having an input raster and would like to clip it with polygons. After running the below code, while overlaying output raster on input raster, I see there is a slight shift in pixels of output raster.

Can somebody help me how to match the output and input raster so that the pixel borders matches exactly?

import ogr
import subprocess
inraster = 'C:/Users/claudio/workspace/test/inraster.tif'
print(inraster)
inshape = 'C:/Users/claudio/workspace/test/country.shp'
ds = ogr.Open(inshape)
lyr = ds.GetLayer(0)
lyr.ResetReading()
ft = lyr.GetNextFeature()
while ft:
country_name = ft.GetFieldAsString('admin')
outraster = inraster.replace('.tif', '_%s.tif' % country_name.replace(' ', '_'))    
subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, '-crop_to_cutline', '-cwhere', "'admin'='%s'" % country_name])
ft = lyr.GetNextFeature()  
ds = None

Answer 1

Keeping the pixels aligned with the original raster's is kind of outside the scope of why gdalwarp exists. It would seem that is more in the territory of gdal_translate, but gdal_translate doesn't have an option to clip to a cutline. Regardless, if you can calculate the optimal bounds of the image ahead of time or in this script (aligned to pixel dimensions of course), then these can be passed in to gdalwarp's te option. te specifies the target extent.

import osr 
import ogr 
import gdal
import subprocess

inraster = 'C:/Users/claudio/workspace/test/inraster.tif'
inshape = 'C:/Users/claudio/workspace/test/country.shp'


# open raster and get its georeferencing information
dsr = gdal.Open(inraster, gdal.GA_ReadOnly)
gt = dsr.GetGeoTransform()
srr = osr.SpatialReference()
srr.ImportFromWkt(dsr.GetProjection())

# open vector data and get its spatial ref
dsv = ogr.Open(inshape)
lyr = dsv.GetLayer(0)
srv = lyr.GetSpatialRef()

# make object that can transorm coordinates
ctrans = osr.CoordinateTransformation(srv, srr)

lyr.ResetReading()
ft = lyr.GetNextFeature()
while ft: 
    # read the geometry and transform it into the raster's SRS
    geom = ft.GetGeometryRef()
    geom.Transform(ctrans)
    # get bounding box for the transformed feature
    minx, maxx, miny, maxy = geom.GetEnvelope()

    # compute the pixel-aligned bounding box (larger than the feature's bbox)
    left = minx - (minx - gt[0]) % gt[1]
    right = maxx + (gt[1] - ((maxx - gt[0]) % gt[1]))
    bottom = miny + (gt[5] - ((miny - gt[3]) % gt[5]))
    top = maxy - (maxy - gt[3]) % gt[5]

    country_name = ft.GetFieldAsString('admin')
    outraster = inraster.replace('.tif', '_%s.tif' % country_name.replace(' ', '_'))        
    subprocess.call([
        'gdalwarp', inraster, outraster,
        '-cutline', inshape, '-cwhere', "'admin'='%s'" % country_name,
        # target resolution
        '-tr', str(abs(gt[1])), str(abs(gt[5])),
        # target extent
        '-te', str(left), str(bottom), str(right), str(top)
    ]) 
    ft = lyr.GetNextFeature()  
ds = None

Answer 2

The standard option for gdalwarp is nearest neighbor interpolation. Nearest neighbor keeps the pixel values untouched, but result in slight shifts in the pixel position. It should be a shift below half a pixel. If it is more some other problem might occur.

You can use the -r resampling method to choose another interpolation. But all will alter your image in some way when the extent or resolution of your image changes. For example bilinear will interpolate on the basis of raster values, which will result in the slight shift of colors but keeps your geometry untouched. Check out the resampling options and what they do. Then you can choose the one that fits your needs the best.

Answer 3

I have encountered the same shift issue. This shift is not limited to GDAL/QGIS. It is also present in ArcMap, MapInfo and Global Mapper (maybe gdal is used in the backend...). A work around is to convert the GeoTiff to an ERS, fix the easting and northing in the .ers file (open with a text editor) and convert the ERS back to a GeoTiff. The conversion between GTiff and ERS with GDAL isn't too slow (about 2-3 min for a 60GB tif). However I am privileged to very good hardware.