gdalwarp with GCPs via GDAL Python bindings

gdal_warp.py

from pathlib import Path
from osgeo import gdal, osr

# Adapted from https://svn.osgeo.org/gdal/trunk/autotest/alg/warp.py
def warp_with_gcps(input_path, output_path, gcps, gcp_epsg=3301, output_epsg=3301):
    # Open the source dataset and add GCPs to it
    src_ds = gdal.OpenShared(str(input_path), gdal.GA_ReadOnly)
    gcp_srs = osr.SpatialReference()
    gcp_srs.ImportFromEPSG(gcp_epsg)
    gcp_crs_wkt = gcp_srs.ExportToWkt()
    src_ds.SetGCPs(gcps, gcp_crs_wkt)

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

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

    # 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)
    dst_xsize = tmp_ds.RasterXSize
    dst_ysize = tmp_ds.RasterYSize
    dst_gt = tmp_ds.GetGeoTransform()
    tmp_ds = None

    # Now create the true target dataset
    dst_path = str(Path(output_path).with_suffix(".tif"))
    dst_ds = gdal.GetDriverByName('GTiff').Create(dst_path, dst_xsize, dst_ysize, src_ds.RasterCount)
    dst_ds.SetProjection(dst_wkt)
    dst_ds.SetGeoTransform(dst_gt)
    dst_ds.GetRasterBand(1).SetNoDataValue(0)

    # And run the reprojection
    gdal.ReprojectImage(src_ds,
                        dst_ds,
                        None,  # src_wkt : left to default value --> will use the one from source
                        None,  # dst_wkt : left to default value --> will use the one from destination
                        resampling,
                        0,  # WarpMemoryLimit : left to default value
                        error_threshold,
                        None,  # Progress callback : could be left to None or unspecified for silent progress
                        None)  # Progress callback user data
    dst_ds = None

input_path = Path("x.tif")
output_path = Path("y.tif")
# GCP input
xyz = [...]
row_col = [...]

gcps = []
for (x, y, z), (row, col) in zip(xyz, row_col):
    gcps.append(gdal.GCP(x, y, z, col, row))

warp_with_gcps(input_path, output_path, gcps, gcp_epsg=3301, output_epsg=3301)

Hello Valgur,
I have tested the script and receive the error

ERROR 6: SetGCPs() is only supported on newly created GeoTIFF files

I had defined GCP input as follows

xyz = [(1,1,5)] row_col = [(1,1)]

Do you have any idea how I could resolve it?

Best
Sebastian

@schoeller
looks like you've opened src_ds in readonly mode

xyz=[(74.2179978544028,19.9340873467462),(85.1321831541614,19.8968271260537),(92.39172194296,29.9992439210911),(80.5381230425356,30.0384277080736)]
row_col=[(11149,0),(1149,6234),(0,6234),(0,0)]
#Assign GCPs
gcps=[]
zp=zip(xyz,row_col)
for (x,y) , (row,col) in zp:
    gcps.append(gdal.GCP(x,y,col,row))
    print(gdal.GCP(x,y,col,row))

#It doesn't seem to be working , I do not understand what the row and col represent