How to use the gdal.GCP function in GDAL

start_dir=os.getcwd()
            if dst_dir:
                os.chdir(dst_dir)

            dst_file = os.path.abspath(os.path.basename(base+ext)) # output file
            dst_drv = gdal.GetDriverByName(out_format)
            dst_ds = dst_drv.CreateCopy(
                dst_file.encode(locale.getpreferredencoding()),
                self.raster_ds,
                0
            )
            dst_ds.SetProjection(self.srs)

            #double x = 0.0, double y = 0.0, double z = 0.0, double pixel = 0.0,
            #double line = 0.0, char info = "", char id = ""
            gcps=[gdal.GCP(c[0],c[1],0,p[0],p[1],'',i) for i,p,c in self.refs]
            dst_ds.SetGCPs(gcps,self.refs.srs())
            dst_geotr=gdal.GCPsToGeoTransform(gcps) # if len(gcps) < 5 else (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
            dst_ds.SetGeoTransform(dst_geotr)
            poly,gmt_data=self.cut_poly(dst_ds)
            if poly:
                dst_ds.SetMetadataItem('CUTLINE',poly)
            if self.name:
                dst_ds.SetMetadataItem('DESCRIPTION',self.name.encode('utf-8'))

            dst_ds = None # close dataset
#            re_sub_file(dst_file, [
#                    ('^.*.*\n',''),
#                    ('^.*.*\n','')
#                    ])
        finally:
            self.raster_ds = None

original data (e.g., first line starts in upper left corner or in lower left corner). For
        Sentinel-1, the raster data is flipped in relation to the GCPs, so we need to flip the GCP
        line vector as well.

        """
        lon = self.ds.variables['GCP_longitude_'+self.ds.polarisation[:2]][:]
        lat = self.ds.variables['GCP_latitude_'+self.ds.polarisation[:2]][:]
        line = self.ds.variables['GCP_line_'+self.ds.polarisation[:2]][:]
        if flip_gcp_line:
            # Flip line vector
            line = self.ds.dimensions['y'].size - line
        pixel = self.ds.variables['GCP_pixel_'+self.ds.polarisation[:2]][:]

        gcps = []
        for i0 in range(0, self.ds.dimensions['gcp_index'].size):
            gcp = gdal.GCP(float(lon[i0]), float(lat[i0]), 0, float(pixel[i0]),
                    float(line[i0]))
            gcps.append(gcp)

        return gcps

i += 1
                dst_gcps.append(gdal.GCP(igcp.GCPX, igcp.GCPY, 0,
                                         igcp.GCPPixel - x_offset,
                                         igcp.GCPLine - y_offset, str(''), str(i)))
        n_gcps = i
        if n_gcps < 100:
            # create new 100 GPCs (10 x 10 regular matrix)
            pix_array, lin_array = np.mgrid[0:x_size:10j, 0:y_size:10j]
            pix_array = pix_array.flatten() + x_offset
            lin_array = lin_array.flatten() + y_offset

            lon_array, lat_array = self.vrt.transform_points(pix_array, lin_array,
                                                dst_srs=NSR(self.vrt.dataset.GetGCPProjection()))

            for i in range(len(lon_array)):
                dst_gcps.append(gdal.GCP(lon_array[i], lat_array[i], 0,
                                         pix_array[i] - x_offset,
                                         lin_array[i] - y_offset,
                                         str(''), str(n_gcps+i+1)))

        # set new GCPss
        self.dataset.SetGCPs(dst_gcps, self.vrt.dataset.GetGCPProjection())
        # remove geotranform which was automatically added
        self._remove_geotransform()

ngt = [gt[0],gt[1],gt[2],gt[3],gt[4],gt[5]]
            ngt[0] = gt[0] + xoff*gt[1] + yoff*gt[2];
            ngt[3] = gt[3] + xoff*gt[4] + yoff*gt[5];
            dst.SetGeoTransform( ( ngt[0], ngt[1], ngt[2], ngt[3], ngt[4], ngt[5] ) )
            
    #Check for GCPs
    elif src.GetGCPCount() > 0:
        
        if (xoff == 0) and (yoff == 0):
            dst.SetGCPs( src.GetGCPs(), src.GetGCPProjection() )
        else:
            gcps = src.GetGCPs()
            #Shift gcps
            new_gcps = []
            for gcp in gcps:
                ngcp = gdal.GCP()
                ngcp.GCPX = gcp.GCPX 
                ngcp.GCPY = gcp.GCPY
                ngcp.GCPZ = gcp.GCPZ
                ngcp.GCPPixel = gcp.GCPPixel - xoff
                ngcp.GCPLine = gcp.GCPLine - yoff
                ngcp.Info = gcp.Info
                ngcp.Id = gcp.Id
                new_gcps.append(ngcp)

            try:
                dst.SetGCPs( new_gcps , src.GetGCPProjection() )
            except:
                print "Failed to set GCPs"
                return

    return

px_height = int(src_image.RasterYSize - ul_y - 1)

            # We've constrained x & y so they are within the image.
            # If the width or height ends up negative at this point,
            # the AOI is completely outside the image
            if px_width < 0 or px_height < 0:
                print('AOI out of range, skipping\n')
                continue

            corners = [[0, 0], [px_width, 0], [px_width, px_height], [0, px_height]]
            corner_names = ['UpperLeft', 'UpperRight', 'LowerRight', 'LowerLeft']
            world_corners = model.back_project(corners, elevation)

            corner_gcps = []
            for (p, l), (x, y, h), n in zip(corners, world_corners, corner_names):
                corner_gcps.append(gdal.GCP(x, y, h, p, l, "", n))

            # Load the source data as a gdalnumeric array
            clip = src_image.ReadAsArray(ul_x, ul_y, px_width, px_height)

            # create output raster
            raster_band = src_image.GetRasterBand(1)
            output_driver = gdal.GetDriverByName('MEM')

            # In the event we have multispectral images,
            # shift the shape dimesions we are after,
            # since position 0 will be the number of bands
            try:
                clip_shp_0 = clip.shape[0]
                clip_shp_1 = clip.shape[1]
                if clip.ndim > 2:
                    clip_shp_0 = clip.shape[1]

--------
        gcsp : List with GDAL GCPs

        """
        # estimate step of GCPs
        gcp_size = np.sqrt(n_gcps)
        step0 = max(1, int(float(lat.shape[0]) / gcp_size))
        step1 = max(1, int(float(lat.shape[1]) / gcp_size))

        # generate list of GCPs
        gcps = []
        k = 0
        for i0 in range(0, lat.shape[0], step0):
            for i1 in range(0, lat.shape[1], step1):
                # create GCP with X,Y,pixel,line from lat/lon matrices
                gcp = gdal.GCP(float(lon[i0, i1]), float(lat[i0, i1]), 0, i1, i0)
                gcps.append(gcp)
                k += 1

        return gcps

tp_lon_name = 'TP_longitude'
    tp_lat_name = 'TP_latitude'
    if tp_x_name in dataset.sizes and tp_y_name in dataset.sizes \
            and tp_lon_name in dataset and tp_lat_name in dataset:
        tp_gcps = []
        tp_width = dataset.sizes[tp_x_name]
        tp_height = dataset.sizes[tp_y_name]
        tp_lon_var = dataset[tp_lon_name]
        tp_lat_var = dataset[tp_lat_name]
        tp_lon_values = tp_lon_var.values
        tp_lat_values = tp_lat_var.values
        gcp_index = 0
        for y in np.linspace(0, tp_height - 1, tp_gcp_step, dtype=np.int32):
            for x in np.linspace(0, tp_width - 1, tp_gcp_step, dtype=np.int32):
                lon, lat = float(tp_lon_values[y, x]), float(tp_lat_values[y, x])
                tp_gcps.append(gdal.GCP(lon, lat, 0.0, x + 0.5, y + 0.5, '%s,%s' % (x, y), str(gcp_index)))
                gcp_index += 1
    else:
        tp_gcps = None
        tp_width = None
        tp_height = None
        tp_lon_var = None

    mem_driver = gdal.GetDriverByName("MEM")

    lon1 = lon_min
    lon2 = lon_min + res * dst_width
    lat1 = lat_max - res * dst_height
    lat2 = lat_max

    dst_dataset = xr.Dataset()
    dst_dataset['lon_bnds'] = (['lon', 'bnds'],

def shift_cropped_gcps(self, x_offset, x_size, y_offset, y_size):
        """Modify GCPs to fit the size/offset of cropped image"""
        gcps = self.dataset.GetGCPs()
        if len(gcps) == 0:
            return

        dst_gcps = []
        i = 0
        # keep current GCPs
        for igcp in gcps:
            if (0 < igcp.GCPPixel - x_offset < x_size and 0 < igcp.GCPLine - y_offset < y_size):
                i += 1
                dst_gcps.append(gdal.GCP(igcp.GCPX, igcp.GCPY, 0,
                                         igcp.GCPPixel - x_offset,
                                         igcp.GCPLine - y_offset, str(''), str(i)))
        n_gcps = i
        if n_gcps < 100:
            # create new 100 GPCs (10 x 10 regular matrix)
            pix_array, lin_array = np.mgrid[0:x_size:10j, 0:y_size:10j]
            pix_array = pix_array.flatten() + x_offset
            lin_array = lin_array.flatten() + y_offset

            lon_array, lat_array = self.vrt.transform_points(pix_array, lin_array,
                                                dst_srs=NSR(self.vrt.dataset.GetGCPProjection()))

            for i in range(len(lon_array)):
                dst_gcps.append(gdal.GCP(lon_array[i], lat_array[i], 0,
                                         pix_array[i] - x_offset,
                                         lin_array[i] - y_offset,