遙感影像重投影

一、簡介

　　柵格數據進行重新投影比矢量數據更復雜，對於矢量，你只需要每個頂點的新坐標就可以輕松實現投影轉換，但對於柵格，你需要處理像元發生形變和偏移的情況，以及從舊單元格位置到新單元格位置的一對一映射。新單元格位置不存在（如下圖）。確定新單元格像元值的最簡單方法是使用最接近輸出單元格映射的輸入單元格中的值，這稱為最近鄰，是最快的方法，通常是分類數據所需的方法。除此之外的所有其他采樣類型都將更改像元值，對於數據分類而言，這是我們不希望看到的結果。但是，如果使用最近鄰方法，連續數據的柵格通常看起來不太美觀，也就是圖像上的物體可能出現“斷裂”。對於這種情況，我們通常使用雙線性插值或三次卷積，它采用的是周圍像元的平均值。

二、重投影方法

　　GDAL是處理遙感影像強大的軟件包，對遙感圖像重投影，下面列舉兩種方法，

1 gdal.warp()

from osgeo import gdal
 
in_ds = gdal.Open(r'F:\algorithm\MOD09A1.A2017361.h28v06.006.2018005034659.hdf')
# 返回結果是一個list，list中的每個元素是一個tuple，每個tuple中包含了對數據集的路徑，元數據等的描述信息
# tuple中的第一個元素描述的是數據子集的全路徑
datasets = in_ds.GetSubDatasets()
 
# 取出第1個數據子集（MODIS反射率產品的第一個波段）進行轉換
# 第一個參數是輸出數據，第二個參數是輸入數據，后面可以跟多個可選項
gdal.Warp('D:/reprojection01.tif', datasets[0][0], dstSRS='EPSG:32649')  # UTM投影
gdal.Warp('D:/reprojection02.tif', datasets[0][0], dstSRS='EPSG:4326')   # 等經緯度投影
 
# 關閉數據集
root_ds = None

2 gdal.AutoCreateWarpedVRT()

除了使用GDAL附帶的gdalwarp實用程序之外，重新投影圖像的最簡單方法是使用VRT。有一個方便的功能，當你提供空間參考信息時，會為你創建重新投影的VRT數據集。具體參數如下：

• src_ds是要重新投影的數據集。
• src_wkt是源空間參考系的WKT表示。默認值為None，在這種情況下，它將使用源柵格中的SRS信息。如果此柵格沒有SRS信息，則需要在此處提供。如果使用None，你也可以在這里提供。
• dst_wkt是所需空間參照系的WKT表示。默認值為None，在這種情況下不會執行重新投影操作。
• eRasampleAlg是表1中的重采樣方法之一。默認值為GRA_NearestNeighbour。
• maxerror是你要允許的最大錯誤量（以像元為單位）。默認值為0，用於精確計算。

重采樣方法：

GRA_NearestNeighbour　　　　選取最鄰近的像元

GRA_Bilinear　　　　鄰近4個像元加權平均

GRA_Cubic　　　　鄰近的16個像元平均

GRA_CubicSpline　　　　16個像元的三次B樣條

GRA_Lanczos　　　　36個像元Lanczos窗口

GRA_Average　　　　求均值

GRA_Mode　　　　出現頻率最多的像元值

　　AutoCreateWarpedVRT函數不會在磁盤上創建VRT文件，但會返回一個數據集對象，然后可以使用CreateCopy將其保存為其他格式。以下示例采用使用UTM空間參考的自然顏色Landsat圖像，創建具有等經緯度投影WGS84的目標SRS的扭曲VRT，並將VRT復制到GeoTIFF：

from osgeo import gdal, osr
 
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
os.chdir(r'F:\algorithm')
old_ds = gdal.Open('lansat8_test.tif')
vrt_ds = gdal.AutoCreateWarpedVRT(old_ds, None, srs.ExportToWkt(),gdal.GRA_Bilinear)
gdal.GetDriverByName('gtiff').CreateCopy('landsat8_test_wgs84.tif', vrt_ds)

　　OGRSpatialReference類和OGRCoordinateTransformation類主要用來提供定義坐標系統（投影和水准面）和轉換坐標。這兩個類都基於OpenGIS的坐標轉換說明，並且使用Well Known Text格式來進行表述坐標系統。

　　一些關於OpenGIS坐標系統的資料，以及空間參考坐標抽象模型文件可以在OGC（Open Geospatial Consortium）的網站上找到。GeoTIFF投影轉換列表（GeoTIFF Projections Transform List）網頁可以更好的幫助你理解WKT的規則，同時EPSG的網站也是很有用的資料。

三、仿射變換

1 GetGeoTransform

#GetGeoTransform，在真實坐標和柵格數據坐標具有相同srs情況下，計算坐標偏移。
#作用：圖像坐標（行列號）和現實世界坐標（投影坐標或地理坐標）之間的轉換。是仿射變換，不是投影轉換，和上面不同。
#0、3 x\y坐標 起始點現實世界坐標  1、5 像素寬度和高度  2、4 x\y方向旋轉角
gt = ds.GetGeoTransform() #正變換：圖像坐標到現實世界坐標。正變換時輸入行列號，輸出的現實世界坐標是柵格圖像srs下的坐標.
inv_gt = gdal.InvGeoTransform(gt) #逆變換：現實世界坐標到圖像坐標
offsets = gdal.ApplyGeoTransform(inv_gt, 465200, 5296000) #逆變換：輸入的投影坐標具有和柵格圖像的相同的srs
xoff, yoff = map(int, offsets)  #取整

2 gdal.Transformer

#gdal.Transformer，可計算相同srs下的坐標偏移；不能用於不同srs投影轉換
#　　作用：現實世界坐標（投影坐標或地理坐標）與圖像坐標（行列號）之間的轉換、兩個柵格圖像之間像素坐標偏移（行列號），
如鑲嵌#原理就是在相同srs情況下，計算圖1的像素坐標到現實世界坐標的偏移，再從現實世界坐標偏移到圖2的像素坐標。
其實就是兩次仿射變換（正、逆），從而把圖1的像素坐標偏移到圖2的像素坐標。#所以不能用於不同srs情況，
因為該函數沒有內置不同srs的投影轉換公式。只能用於相同srs下，兩個柵格數據集坐標的偏移。
#這里in_ds和out_ds具有相同srs。轉換目的是為了把不同柵格數據的圖像坐標（行列號）進行偏移，方便鑲嵌
#in_ds是in_ds = gdal.Open()。
trans = gdal.Transformer(in_ds, out_ds, [])  #空白用於設置轉換器選項
success, xyz = trans.TransformPoint(False, 0, 0) #False基於源數據計算目標柵格，反之為True。起始坐標為左上角 0，0
x,y,z = map(int, xyz)

#圖像坐標和現實世界坐標之間轉換
trans = gdal.Transformer(out_ds, None, [])
success, xyz = trans.TransformPoint(0, 1078, 648)