The geotransform read by gdal from a TerraSAR-X Tiff is incorrect

Question

I am using python's gdal module to analyse satellite images - RADARSAT-2 and TerraSAR-X - saved as .tif files. I need to fetch pixel values at coordinates read from a shapefile. While the code works fine for the RS2 images, I'm having trouble with the TSX images.

The geotransform read by gdal is off for the TSX products, which yields negative pixel indices for the location of shapefile features on the image. The same piece of code works fine for RS2 products.

Any idea of what's going on and how to fix it ?

Example of a correct geotransform from a RS2 product:

(-74.98992283355103, 7.186522272956171e-05, 0.0, 62.273587708987776, 0.0, -7.186522272956171e-05)

Example of the geotransforms I get for TSX products :

(506998.75, 2.5, 0.0, 6919001.25, 0.0, -2.5)

Code snippet :

import gdal

gdal.UseExceptions()

# Read image, band, geotransform
dataset = gdal.Open(paths['TSXtiff'])
band = dataset.GetRasterBand(band_index)        
gt = dataset.GetGeoTransform()

# Read shapefile
shapefile = ogr.Open(paths["Shapefile"])    
layer = shapefile.GetLayer()

# Add pixel_value for each feature to associated list
pixels_at_shp = []
for feature in layer :

    geometry = feature.GetGeometryRef()             

    # Coordinates in map units
    # In GDAL, mx = px*gt[1] + gt[0], my = py*gt[5] + gt[3]
    mx,my = geometry.GetX(), geometry.GetY()

    # Convert to pixel coordinates
    px = int((mx-gt[0])/gt[1])
    py = int((my-gt[3])/gt[5])

    band_values = band.ReadAsArray(px,py,1,1)

    pixels_at_shp.append(band_values)

shapefile = None

return pixels_at_shp

Answer 1

The TSX product was projected while the RS2 product was simply georeferenced. Solution : go back to lat/long degrees coordinates by reprojecting the GeoTIFF files onto WGS84.

I used the gdal command line tool to do the reprojection, like so :

for f in *.tif
do 
    gdalwarp "$f" "${f%%.*}_reproj.tif" -t_srs "+proj=longlat +ellps=WGS84"
done