私は衛星画像(RADARSAT-2とTerraSAR-X)を.tifファイルとして保存するためにpythonのgdalモジュールを使用しています。シェイプファイルから読み取った座標でピクセル値を取得する必要があります。このコードはRS2イメージで正常に機能しますが、TSXイメージに問題があります。TerraSAR-X Tiffからgdalによって読み込まれたジオトランスフォームが正しくありません
gdalで読み取られたgeotransformは、TSX製品ではオフになっており、画像上のシェイプファイルのフィーチャの位置に対して負のピクセルインデックスが生成されます。同じコードがRS2製品でうまく動作します。
何が起こっているのか、それを修正する方法はありますか? RS2製品から正しいgeotransformの
例:
(506998.75, 2.5, 0.0, 6919001.25, 0.0, -2.5)
コードスニペット:私はTSX製品の取得geotransformsの
(-74.98992283355103, 7.186522272956171e-05, 0.0, 62.273587708987776, 0.0, -7.186522272956171e-05)
例
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