2017-06-14 2 views
0

私は衛星画像(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 

答えて

0

TSX RS2製品は単純にジオリファレンスされていました。解決策:GeoTIFFファイルをWGS84に再投影して緯度/経度座標に戻ります。

私はそうのように、再投影を行うためのGDALコマンドラインツールを使用し:

for f in *.tif 
do 
    gdalwarp "$f" "${f%%.*}_reproj.tif" -t_srs "+proj=longlat +ellps=WGS84" 
done 
関連する問題