def geoGCP(dataFilename, outGeoName, rpcFileName):
    from osgeo import  gdal
    from osgeo import  osr
    import geo_rpc_copy

    dataset = gdal.Open(dataFilename, gdal.GA_Update)
    
    # 获取输入数据集的数据类型
    inputDataType = dataset.GetRasterBand(1).DataType

    rpc = geo_rpc_copy.read_rpc_file(rpcFileName)

    rows = dataset.RasterXSize
    cols = dataset.RasterYSize
    TLlon, TLlat = rpc.localization(0, 0, 0) # 左上角
    TRlon, TRlat = rpc.localization(rows-1, 0, 0) # 右上角
    BLlon, BLlat = rpc.localization(0, cols-1, 0) # 左下角
    BRlon, BRlat = rpc.localization(rows-1, cols-1, 0) # 右下角
    Clon, Clat = rpc.localization(int(round(rows/2)-1), int(round(cols/2)-1), 0)
    Qlon41, Qlat41 = rpc.localization(int(round(rows/4)-1), int(round(cols/4)-1), 0)
    Qlon43, Qlat43 = rpc.localization(int(round(rows*3/4)-1), int(round(cols*3/4)-1), 0)
    Qlon41_43, Qlat41_43 = rpc.localization(int(round(rows/4)-1), int(round(cols*3/4)-1), 0)
    Qlon43_41, Qlat43_41 = rpc.localization(int(round(rows*3/4)-1), int(round(cols/4)-1), 0)
    Qlon21_41, Qlat21_41 = rpc.localization(int(round(rows/2)-1), int(round(cols/4)-1), 0)
    Qlon21_43, Qlat21_43 = rpc.localization(int(round(rows/2)-1), int(round(cols*3/4)-1), 0)
    Qlon41_21, Qlat41_21 = rpc.localization(int(round(rows/4)-1), int(round(cols/2)-1), 0)
    Qlon43_21, Qlat43_21 = rpc.localization(int(round(rows*3/4)-1), int(round(cols/2)-1), 0)

    gcps = [gdal.GCP(TLlon, TLlat, 0, 0, 0),
            gdal.GCP(TRlon, TRlat, 0, rows-1, 0),
            gdal.GCP(BLlon, BLlat, 0, 0, cols-1),
            gdal.GCP(BRlon, BRlat, 0, rows-1, cols-1),
            gdal.GCP(Qlon41, Qlat41, 0, int(round(rows/4)-1), int(round(cols/4))-1),
            gdal.GCP(Qlon43, Qlat43, 0, int(round(rows*3/4)-1), int(round(cols*3/4))-1),
            gdal.GCP(Clon, Clat, 0, int(round(rows/2)-1), int(round(cols/2))-1),
            gdal.GCP(Qlon41_43, Qlat41_43, 0, int(round(rows/4)-1), int(round(cols*3/4))-1),
            gdal.GCP(Qlon43_41, Qlat43_41, 0, int(round(rows*3/4)-1), int(round(cols/4))-1),
            gdal.GCP(Qlon21_41, Qlat21_41, 0, int(round(rows/2)-1), int(round(cols/4))-1),
            gdal.GCP(Qlon21_43, Qlat21_43, 0, int(round(rows/2)-1), int(round(cols*3/4))-1),
            gdal.GCP(Qlon41_21, Qlat41_21, 0, int(round(rows/4)-1), int(round(cols/2))-1),
            gdal.GCP(Qlon43_21, Qlat43_21, 0, int(round(rows*3/4)-1), int(round(cols/2))-1),
            ]
    sr = osr.SpatialReference()
    sr.SetWellKnownGeogCS('WGS84')
    # 添加控制点
    dataset.SetGCPs(gcps, sr.ExportToWkt())
    # tps校正 重采样
    dst_ds = gdal.Warp(outGeoName, dataset, format='GTiff', tps=True,
                       resampleAlg=gdal.GRIORA_NearestNeighbour, outputType=inputDataType)
    # 备选重采样方法： 最近邻 gdal.GRIORA_NearestNeighbour
    # gdal.GRIORA_Bilinear