2017-06-14 64 views
0

我使用Python的GDAL模塊來分析衛星圖像 - RADARSAT-2和TerraSAR-X衛星 - 另存爲.TIF文件。我需要從shapefile中讀取座標處的像素值。雖然代碼對於RS2圖像正常工作,但我在TSX圖像上遇到問題。通過從TerraSAR-X衛星的Tiff @中讀取的地理座標是不正確

通過GDAL讀取的地理座標轉換爲關閉用於TSX產品,其產生用於shape文件的特徵的圖像上的位置的負像素索引。對於RS2產品,相同的代碼工作正常。

有什麼想法以及如何解決?從RS2產品正確的地理座標的

例子:

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

代碼片段:

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

的geotransforms我獲得TSX產品的實例

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 
相關問題