Beijing 1954 (TIFF) ➜ WGS 84 (EPSG:4326) — Raster Re-projection

import os
from osgeo import gdal, osr
import pandas as pd
from typing import List

# ---- User inputs -------------------------------------------------------
input_tiff   = "beijing54.tif"         # Source raster (Beijing 54)
gcps_csv     = "54_to_84.csv"          # Control-point table (CSV)
output_tiff  = "beijing54_to_wgs84.tif"

# Transform method:  'tps' (thin-plate spline) or polynomial order (int)
transform_method = "tps"   # choose "tps" or 1/2/3 for polynomial order

# The CSV must contain:
#   x54, y54  ▶ coordinates in Beijing 54 (same units as the raster)
#   lon84, lat84 ▶ WGS 84 geographic coordinates (decimal degrees)
gdf = pd.read_csv(gcps_csv)

# Open source raster
src_ds = gdal.Open(input_tiff, gdal.GA_ReadOnly)
if src_ds is None:
    raise FileNotFoundError(f"Cannot open {input_tiff}")

gt = src_ds.GetGeoTransform()
inv_gt = gdal.InvGeoTransform(gt)

# Build GCP list
gdal_gcps: List[gdal.GCP] = []
for _, row in gdf.iterrows():
    x54, y54 = float(row['x54']), float(row['y54'])
    lon84, lat84 = float(row['lon84']), float(row['lat84'])
    # Convert world → pixel
    px, py = gdal.ApplyGeoTransform(inv_gt, x54, y54)
    gcp = gdal.GCP(lon84, lat84, 0, px, py)  # (lon, lat, z, pixel, line)
    gdal_gcps.append(gcp)

print(f"Loaded {len(gdal_gcps)} control points.")

# Create an in-memory VRT copy with GCPs
vrt_path = "/vsimem/tmp_with_gcps.vrt"
gdal.Translate(
    vrt_path,
    src_ds,
    GCPs=gdal_gcps,
    outputSRS="EPSG:4326"  # Target GCP spatial ref (WGS 84)
)
src_ds = None  # close

warp_kwargs = dict(
    dstSRS="EPSG:4326",
    format="GTiff",
    xRes=0.0002695,   # ≈ 30 m at the equator — adjust as needed
    yRes=0.0002695,
    resampleAlg="bilinear",
    multithread=True
)

if transform_method == "tps":
    warp_kwargs["tps"] = True
else:
    warp_kwargs["polynomialOrder"] = int(transform_method)

gdal.Warp(
    destNameOrDestDS=output_tiff,
    srcDSOrSrcDSTab=vrt_path,
    **warp_kwargs
)

print(f"Warp complete ➜ {output_tiff}")

ds = gdal.Open(output_tiff)
print("Target CRS:", ds.GetSpatialRef().ExportToWkt().split(',')[0])
print("Raster size :", ds.RasterXSize, "×", ds.RasterYSize)
print("GeoTransform:", ds.GetGeoTransform())
ds = None