gdal 从 TerraSAR-X Tiff 读取的地理变换不正确答案

【问题标题】：The geotransform read by gdal from a TerraSAR-X Tiff is incorrectgdal 从 TerraSAR-X Tiff 读取的地理变换不正确
【发布时间】：2017-06-14 18:02:33
【问题描述】：

我正在使用 python 的 gdal 模块来分析卫星图像 - RADARSAT-2 和 TerraSAR-X - 保存为 .tif 文件。我需要在从 shapefile 读取的坐标处获取像素值。虽然代码适用于 RS2 图像，但我在使用 TSX 图像时遇到了问题。

gdal 读取的地理变换对于 TSX 产品已关闭，这会为图像上 shapefile 特征的位置产生负像素索引。同一段代码适用于 RS2 产品。

知道发生了什么以及如何解决它吗？

来自 RS2 产品的正确地理变换示例：

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

我为 TSX 产品获得的地理变换示例：

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

代码sn-p：

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

【问题讨论】：

标签： python gdal geotiff

【解决方案1】：

TSX 产品是经过投影的，而 RS2 产品只是简单的地理参考。解决方案：通过将 GeoTIFF 文件重新投影到 WGS84 上，返回到纬度/经度坐标。

我使用 gdal 命令行工具进行重投影，如下所示：

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

【讨论】：