Ich benutze Pythons Gdal-Modul zur Analyse von Satellitenbildern - RADARSAT-2 und TerraSAR-X - gespeichert als .tif-Dateien. Ich muss Pixelwerte an Koordinaten aus einer Shapefile abrufen. Während der Code für die RS2-Bilder gut funktioniert, habe ich Probleme mit den TSX-Bildern.Die von Gdal von einem TerraSAR-X Tiff gelesene Geotransform ist falsch
Die von gdal gelesene Geotransformation ist für die TSX-Produkte deaktiviert, was zu negativen Pixelindizes für die Position der Shapefile-Features auf dem Bild führt. Das gleiche Stück Code funktioniert gut für RS2-Produkte.
Irgendeine Idee von dem, was vor sich geht und wie man es repariert?
Beispiel eines korrekten Geotransformation von einem RS2 Produkt:
(-74.98992283355103, 7.186522272956171e-05, 0.0, 62.273587708987776, 0.0, -7.186522272956171e-05)
Beispiel der geotransforms ich für TSX Produkte erhalten:
(506998.75, 2.5, 0.0, 6919001.25, 0.0, -2.5)
-Code-Schnipsel:
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