gdalwarp -tps producing artifacts at image edge

[dguenzel]

Hello,

I am trying to gdalwarp some Sentinel-1 EW GRD products to EPSG:3413. The acquisitions cover the arctic region north of Greenland. If I set the -tps option I see some strange artifacts(?) in some of the produced GeoTiffs as shown below. Could this be due to a bug in the Thin Plate Spline algorithm? Thank you.

Expected behavior and actual behavior.

Warped GeoTiffs contain artifacts where there should be nodata, for example:

S1B_EW_GRDM_1SDH_20211208T160438_20211208T160538_029936_0392E4_16C7

S1B_EW_GRDM_1SDH_20211210T154815_20211210T154915_029965_0393CB_56B1

EDIT: There is also visible distortion at the image edge (top right).

Steps to reproduce the problem.

gdalwarp -t_srs EPSG:3413 -r near -of GTiff -tps s1b-ew-grd-hh-20211208t160438-20211208t160538-029936-0392e4-001.tiff s1b-ew-grd-hh-20211208t160438-20211208t160538-029936-0392e4-001_proj_tps_ver342.tiff

gdalwarp -t_srs EPSG:3413 -r near -of GTiff -tps s1b-ew-grd-hh-20211210t154815-20211210t154915-029965-0393cb-001.tiff s1b-ew-grd-hh-20211210t154815-20211210t154915-029965-0393cb-001_proj_tps_ver342.tiff

Here is the gdalinfo of the created files: gdalinfo.txt

If the -tps option is not given, the problem does not occur. Warping to a different CRS, EPSG:3995 or 5938, showed the same behaviour.

Operating system

Ubuntu 20.04.4 LTS

GDAL version and provenance

3.4.2 version from conda-forge

Tested other versions showing the same problem: 3.3.3 version from conda-forge 3.3.1 included with QGIS (Windows x64)

[jratike80]

Add also gdalinfo of the source image.

[dguenzel]

Sorry, here it is: gdalinfo_source.txt

[rouault]

it would be more convenient if you could attach the zipped output of:

gdal_translate s1b-ew-grd-hh-20211208t160438-20211208t160538-029936-0392e4-001.tiff out.tif -scale 0 1 255 255 -ot Byte -co COMPRESS=LZW -co TILED=YES

which will produce an image with the GCPs but with all pixels set to 255, which should compress very well

[dguenzel]

Sure, good idea: out.zip

[scottstanie]

I think the algorithm for warping with tps around the poles has the issue. @taliboliver ran into a similar issue warping a Sentinel-1 GRD to EPSG:3413

I had the idea to transform all the GCPs to EPSG:3413, then gdalwarp like normal and it seems to work. Here's the warped file in QGIS (shown in the EPSG:3413 projection)

I did this with the following short script

import shutil
import subprocess
import sys
from pathlib import Path

import rasterio as rio
from pyproj import CRS, Transformer

def warp_gcps(
    gcps: list[rio.control.GroundControlPoint],
    dst_epsg: int = 3413,
    src_epsg: int = 4326,
):
    """Warp the GCPs from a Rasterio raster to `dst_epsg`."""
    t = Transformer.from_crs(src_epsg, dst_epsg, always_xy=True)
    xs = [g.x for g in gcps]
    ys = [g.y for g in gcps]
    out_xs, out_ys = t.transform(xs, ys)
    for x, y, g in zip(out_xs, out_ys, gcps):
        g.x = x
        g.y = y
    return gcps

if __name__ == "__main__":
    if len(sys.argv) < 3:
        print(f"Usage: {sys.argv[0]} input_file output_file")
        sys.exit(1)

    infile = sys.argv[1]
    outfile = sys.argv[2]
    # We're altering the GCPs, so make a copy of the file
    tmp_infile = "temp.tif"
    shutil.copy(infile, tmp_infile)

    dst_epsg = 3413
    with rio.open(tmp_infile, "r+") as src:
        src_gcps, src_crs = src.get_gcps()
        warped_gcps = warp_gcps(src_gcps, src_epsg=src_crs.to_epsg(), dst_epsg=dst_epsg)
        warped_crs = CRS.from_epsg(dst_epsg)
        src.gcps = warped_gcps, warped_crs

    cmd = f"gdalwarp -tps -r bilinear -tr 100 100 -srcnodata 0 -dstnodata 0 -t_srs EPSG:{dst_epsg} {tmp_infile} {outfile}"
    print(cmd)
    subprocess.run(cmd, check=True, shell=True)
    Path(tmp_infile).unlink()